rem Note it requires input data file
rem Yes I like line numbers!
5 DEFINT I-N
IM = 64
10 DIM AH(IM, 7) AS INTEGER, AX(IM, IM) AS INTEGER
20 DIM TP(1024), F(1024), X(IM, IM), E(IM, IM)
30 DIM aj(21, 21), VH(7), XX(1024), Y(1024), YD(30), XD(40)
40 OPEN "I", #1, "c:\NMR2.DAT"
LINE INPUT #1, B$
rem reading number of protons
50 INPUT #1, N
rem reading chem. shifts
60 LINE INPUT #1, Z$
70 GOSUB 5000
80 REM REM PRINT "INPUT SPIN-SPIN COUPLING J12,J13...J23,J24..etc"
90 LINE INPUT #1, Z$
100 GOSUB 6000
110 INPUT #1, FMIN, FMAX, t2
185 TK = TIMER
188 REM CALCULATING ELEMENTS IN HAMILTONIAN MATRIX
189 REM VH(2) = 10362!
REM INPUT t2
xsca = 570
ysca = 200
190 N2 = 2 ^ N
210 FOR I = 0 TO N2 - 1
220 II = I
230 FOR IN = N - 1 TO 0 STEP -1
240 JJ = 2 ^ IN
250 IF (II - JJ) < 0 THEN TJ = -1 ELSE TJ = 1
270 IF TJ = 1 THEN II = II - JJ
280 AH(I + 1, IN + 1) = TJ
300 NEXT IN
310 NEXT I
600 FOR I = 1 TO N2
605 E(I, I) = 0: AX(I, I) = 0
610 FOR L = 1 TO N
620 E(I, I) = E(I, I) + VH(L) * .5 * AH(I, L)
625 TD = AH(I, L) * .25
630 FOR K = L + 1 TO N
640 E(I, I) = E(I, I) + aj(L, K) * TD * AH(I, K)
650 NEXT K
660 NEXT L
670 NEXT I
680 FOR I = 1 TO N2
685 X(I, I) = 1
690 FOR L = I + 1 TO N2
700 MA = 0: MB = 0: K1 = 0: K2 = 0
710 X(I, L) = 0: X(L, I) = 0
740 FOR K = 1 TO N
745 TG = AH(I, K) * AH(L, K) - 1
750 MA = MA + TG
760 MB = MB + TG * AH(L, K)
770 IF TG = 0 GOTO 820
780 IF K1 > 0 THEN K2 = K ELSE K1 = K
820 NEXT K
830 IF MA <> -4 GOTO 870
840 IF MB <> 0 GOTO 870
850 E(I, L) = .5 * aj(K1, K2)
860 E(L, I) = E(I, L)
870 NEXT L
880 NEXT I
GOTO 890
FOR LQ = 1 TO N2
FOR KQ = 1 TO N2
REM PRINT E(LQ, KQ);
NEXT KQ
REM PRINT ""
NEXT LQ
INPUT SA$
890 GOSUB 4000
895 REM PRINT " "; L; " Iterations, time required="; TIMER - TK; " Sec"
960 REM CALCULATING FREQUENCIES & TRANSISTION PROBABILITIES
970 IAY = 0
980 FOR I = 1 TO N2
990 FOR J = I + 1 TO N2
1000 FOR K = 1 TO N
1010 IAY = IAY + AH(I, K) * AH(J, K)
1020 NEXT K
1040 IF IAY = N - 2 THEN AX(I, J) = 1 ELSE AX(I, J) = 0
1050 AX(J, I) = AX(I, J)
1100 IAY = 0
1110 NEXT J
1120 NEXT I
1130 M = 0
1170 FOR L = 1 TO N2 - 1
1180 FOR K = L + 1 TO N2
1185 FF = ABS(E(K, K) - E(L, L))
1186 IF ABS(10000 - FF) > 50000 GOTO 1320
1200 TPP = 0
1210 FOR J = 1 TO N2
1218 IF ABS(X(L, J)) < .001 GOTO 1260
1220 FOR I = 1 TO N2
1230 IF AX(J, I) = 0 GOTO 1250
1245 TPP = TPP + X(K, I) * X(L, J)
1250 NEXT I
1260 NEXT J
1270 IF ABS(TPP) < .1 GOTO 1320
1280 M = M + 1
1300 F(M) = FF - 10000
1310 TP(M) = TPP * TPP
1320 NEXT K
1330 NEXT L
1350 GOSUB 1500
1380 GOSUB 2700
1390 GOSUB 3140
1395 CLOSE #1
1400 END
1500 FOR I = 1 TO M - 1
1505 FM = F(I): JM = I
1510 FOR J = I + 1 TO M
1520 IF F(J) < FM GOTO 1540
1530 FM = F(J): JM = J
1540 NEXT J
1550 F(JM) = F(I): F(I) = FM
1560 TPT = TP(JM)
1570 TP(JM) = TP(I): TP(I) = TPT
1580 NEXT I
1582 REM PRINT "Time taken after TP calculation="; TIMER - TK; " Sec"
1585 IF M > 20 GOTO 1620
1586 REM PRINT : REM PRINT " Freq. Intensity"
1590 FOR I = 1 TO M
1595 REM PRINT USING "##"; I;
1596 REM PRINT USING "#####.###"; F(I); TP(I)
IF F(I) < 0 THEN F(I) = -F(I)
1600 NEXT I
1610 GOTO 2000
1620 MMM = INT(M / 2)
1621 REM PRINT
1622 REM PRINT " Freq. Intensity Freq. Intensity"
1630 FOR I = 1 TO MMM
1640 REM PRINT USING "###"; I;
1642 REM PRINT USING "###.###"; F(I); TP(I);
1644 REM PRINT " ";
1645 REM PRINT USING "###"; I + MMM;
1646 REM PRINT USING "###.###"; F(I + MMM); TP(I + MMM)
1650 NEXT I
1660 IF MMM = M / 2 GOTO 2000
1662 REM PRINT USING "###"; M;
1664 REM PRINT USING "###.###"; TP(M); F(M)
2000 REM PRINT :
RETURN
2700 REM CONSTRUCTING LORENTZIANS
2930 REM INPUT #1, FMIN, FMAX, T2
2940 SW = FMAX - FMIN
2950 MM = 639
2960 AI = SW / MM
2970 FOR N5 = 1 TO MM + 1
2980 XX(N5) = FMIN + AI * (N5 - 1)
2990 Y(N5) = 0
3000 NEXT N5
3012 WID = 4 / t2
3015 TA = 4 * 3.14159 * 3.14159 * t2 * t2
3017 TB = WID / AI
3020 FOR N5 = 1 TO M
3025 C2 = (F(N5) - FMIN) / AI
IF ABS(F(N5)) > 4000 GOTO 3130
IF ABS(C2 - TB) > 30000 THEN IJC = 30000
IF IJC = 30000 GOTO 3040
3030 IJC = INT(C2 - TB)
3040 IF IJC < 1 THEN IJC = 1
3050 IF IJC > MM GOTO 3130
3060 IJD = INT(C2 + TB)
3070 IF IJD > MM THEN IJD = MM
3075 C1 = 2 * t2 * TP(N5)
3080 FOR MI = IJC TO IJD
3090 C3 = F(N5) - XX(MI)
3110 Y(MI) = Y(MI) + C1 / (1 + TA * C3 * C3)
3120 NEXT MI
3130 NEXT N5
3135 RETURN
3140 REM PLOTTING OUT ON SCREEN
GOSUB 7000
GOTO 3290
3150 YMIN = Y(5)
3160 YMAX = Y(5)
3170 FOR I = 1 TO MM
3180 IF Y(I) > YMAX THEN YMAX = Y(I)
3190 IF Y(I) < YMIN THEN YMIN = Y(I)
3200 NEXT I
3210 SCREEN 2
3220 CLS
REM y(1); ymax: INPUT sa$
3230 YB = 199 - (Y(1) / YMAX) * 199
3235 IB = 1
3240 FOR I = 2 TO MM
3250 YA = 199 - (Y(I) / YMAX) * 199
3255 REM IF ABS(YA - YB) < .005 GOTO 3280
3260 LINE (640 - IB, YB)-(640 - I, YA)
3270 YB = YA
3275 IB = I
3280 NEXT I
3282 LINE (640 - IB, YB)-(640 - MM, YA)
3283 REM PRINT TIMER - TK
LOCATE 1, 1
REM PRINT aj(2, 3)
3284 LINE INPUT SA$
3286 SCREEN 0
3288 CLS
3290 RETURN
4000 REM DIAGONALISATION USING JACOBI METHOD (D.S. 1993)
4005 M = (N2 * N2 - N2) / 2
4010 TRA = M * .001
4090 FOR L = 1 TO M * 2
4100 Q = 0: AM = -1
4110 FOR J = 1 TO N2 - 1
4120 FOR I = J + 1 TO N2
4130 Z = ABS(E(I, J))
4140 Q = Q + Z
4150 IF Z < AM GOTO 4170
4160 AM = Z: IX = I: JX = J
4170 NEXT I
4180 NEXT J
4190 IF Q < TRA THEN RETURN
4200 A1 = E(IX, IX) - E(JX, JX)
4210 A5 = E(IX, JX)
4220 IF ABS(A5) * 50 > ABS(A1) GOTO 4270
4230 A6 = 2 * A1 * A1 + 3 * A5 * A5
4240 C = -2 * A1 * A5 / A6
4250 S = 1 - A5 * A5 / A6
4260 GOTO 4320
4270 A2 = SQR(A1 * A1 + 4 * A5 * A5)
4280 IF A1 > 0 THEN A2 = -A2
4290 A3 = (A1 + A2) / (2 * A5)
4300 S = 1 / SQR(A3 * A3 + 1)
4310 C = A3 * S
4320 FOR J = 1 TO N2
4330 TM = E(JX, J) * S + E(IX, J) * C
4340 E(IX, J) = E(IX, J) * S - E(JX, J) * C
4350 E(JX, J) = TM
4360 TM = X(JX, J) * S + X(IX, J) * C
4370 X(IX, J) = X(IX, J) * S - X(JX, J) * C
4380 X(JX, J) = TM
4390 NEXT J
4400 FOR J = 1 TO N2
4410 TM = E(J, JX) * S + E(J, IX) * C
4420 E(J, IX) = E(J, IX) * S - E(J, JX) * C
4430 E(J, JX) = TM
4440 NEXT J
4450 NEXT L
4460 PRINT "CONVERGENCE FAILURE": CLOSE #1
4490 END
5000 REM SEPARATING INPUT CHARACTER INTO NUMBERS
5010 FOR L = 1 TO N - 1
5020 II = LEN(Z$)
5030 JJ = INSTR(Z$, ",")
5040 VH(L) = VAL(LEFT$(Z$, JJ - 1)) + 10000
5050 Z$ = RIGHT$(Z$, II - JJ)
5060 NEXT L
5070 VH(N) = VAL(Z$) + 10000
5080 RETURN
5090 END
6000 REM SEPARATING INPUT CHARACTER INTO NUMBERS
6010 FOR I = 1 TO N - 2
6015 FOR J = I + 1 TO N
6020 II = LEN(Z$)
6030 JJ = INSTR(Z$, ",")
6040 aj(I, J) = VAL(LEFT$(Z$, JJ - 1))
6050 Z$ = RIGHT$(Z$, II - JJ)
6060 NEXT J
6065 NEXT I
6070 aj(N - 1, N) = VAL(Z$)
6080 RETURN
6090 END
7000 SCREEN 12
WIDTH 80, 60
LINE (0, 0)-(640, 480), 0, BF
B$ = B$ + "(100 MHz)"
7025 YMIN = 5000!: YMAX = -1!: XMIN = 5000!: XMAX = 0!
C1 = (FMAX - FMIN) / 640: I = 640
7030 FOR K = 1 TO I
REM PRINT Y(K); K
IF YMIN > Y(K) THEN YMIN = Y(K)
IF YMAX < Y(K) THEN KMAX = K
IF YMAX < Y(K) THEN YMAX = Y(K)
XX(K) = (C1 * K + FMIN) / 100
IF XMIN > XX(K) THEN XMIN = XX(K)
IF XMAX < XX(K) THEN XMAX = XX(K)
NEXT K
REM PRINT YMAX; YMIN; X(KMAX); Y(KMAX): INPUT SA$
XO = XX(1): YO = Y(1)
XO = ((XMAX - XX(1)) / (XMAX - XMIN) * xsca) + 60
YO = ((Y(1) - YMIN) / (YMAX - YMIN) * ysca) + 30
XO = INT(XO + .5): YO = INT(YO + .5)
FOR K = 2 TO I
XX = ((XMAX - XX(K)) / (XMAX - XMIN) * xsca) + 60
YY = ((Y(K) - YMIN) / (YMAX - YMIN) * ysca) + 30
YN = (YY + YO * TC) / (1 + TC)
XX = INT(XX + .5): YY = INT(YN + .5)
LINE (XO, 480 - YO)-(XX, 480 - YY)
XO = XX: YO = YN
NEXT K
ytop = 440
LINE (60, ytop - ysca)-(60, 455)
LINE (60, 455)-(xsca + 60, 455)
LINE (60, ytop - ysca)-(xsca + 60, ytop - ysca)
LINE (xsca + 60, 455)-(xsca + 60, ytop - ysca)
L = 0
IF N = 0 THEN kk = 10 ELSE kk = 50
IF N > 1 THEN kk = 100
FOR I = -5 TO 20
L = L + 1
YD(L) = I / kk
NEXT I
KG = -1
7040 K = 0
KG = KG + 1
G(0) = .05: G(1) = .1: G(2) = .2: G(3) = .5: G(4) = 1: G(6) = 20: G(7) = 50: G(8) = 100: G(9) = 200
IF XMIN = 0 THEN XS = -G(KG)
IF XMIN <> 0 THEN XS = -G(KG)
XI = XS
7045 XI = XI + G(KG)
IF XI > XMAX GOTO 7047
K = K + 1
IF K > 20 GOTO 7040
XD(K) = XI
GOTO 7045
7047 KG = 0
7050 KY = 0
KG = KG + 1
G(1) = .1: G(2) = .2: G(3) = .5: G(4) = 1: G(6) = 2: G(7) = 50: G(8) = 100: G(9) = 500
IF YMIN = 0 THEN YS = 0
IF YMIN <> 0 THEN YS = G(KG)
FOR I = YS TO YMAX STEP G(KG)
KY = KY + 1
IF KY > 30 GOTO 7050
YD(KY) = I
NEXT I
GOTO 7091
7060 FOR I = 1 TO KY
IF YD(I) > YMAX GOTO 7090
IF YD(I) < YMIN GOTO 7090
YYD = (((YD(I) - YMIN) / (YMAX - YMIN)) * ysca) + 30
YC = INT(((YYD - 30) / ysca) * 54.5 + 4)
YYD = INT(YYD + .5)
REM PRINT YD(I); YYD; YMAX; YMIN
REM PRINT YC: INPUT SA$
IF 60 - YC > 59 GOTO 7090
IF 60 - YC < 1 GOTO 7090
LOCATE 60 - YC, 3
PRINT YD(I)
LINE (53, 480 - YYD)-(60, 480 - YYD)
7090 NEXT I
IF N < 4 THEN AA$ = A$(N + 1) ELSE AA$ = CHR$(N + 48) + " th Derivative"
FF = LEN(AA$)
FOR I = 1 TO FF
LOCATE 32 - I, 1
PRINT MID$(AA$, FF + 1 - I, 1);
NEXT I
7091 FOR I = 1 TO K
REM PRINT XD(I); XMIN: INPUT SA$
IF XD(I) > XMAX GOTO 7100
IF XD(I) < XMIN GOTO 7100
XXD = (((XMAX - XD(I)) / (XMAX - XMIN)) * xsca) + 60
XC = INT(((XXD - 20) / xsca) * xsca / 8 + 2)
XXD = INT(XXD + .5)
REM PRINT YD(I); YYD; YMAX; YMIN
REM PRINT YC: INPUT SA$
IF XC > 75 GOTO 7100
IF XC < 0 GOTO 7100
LOCATE 59, XC
PRINT INT(XD(I) * 10 + .5) / 10;
LINE (XXD, 455)-(XXD, 460)
7100 NEXT I
A$ = "Chemical Shift ë (ppm)"
FOR I = 1 TO 25
LOCATE 60, 30 + I
PRINT MID$(A$, I, 1);
NEXT I
ypos = (ytop - ysca) / 8
xpos = xsca / 12
LOCATE ypos + 2, xpos - 15
PRINT B$;
7500 CLOSE #1
7502 H$ = INKEY$
IF H$ = "" GOTO 7502
7510 RETURN
|