10    REM/***************************************************************************
20    REM/THIS IS A BENEAL PURPOSE PROGRAM TO SOLVE 1-DIMENSIONAL
30    REM/DIFFUSION PROBLEM IN THE FORM OF:
40    REM/
50    REM/         "RCdT/dt=1/A(X)d/dX(A(X)KdT/dX)+S"
60    REM/
70    REM/-------------------MAIN PROGRAM--------------------------------------------
80    DIM X(130), XF(130), XM(130), XP(130), R(130), RF(130), AP(130)
90    DIM AE(130), AW(130), CN(130), P(130), Q(130), T(130), TA(130)
100   DIM TG(130), TR(130), TS(130), GM(130), RC(130), H#(30)
110   K = 1: KM = 1: KP = 1: OM = 1
120   NF = 1
130   T = 0
140   KN = 1
150   KT = 1
160   REM/-----------------FIRST TO SPECIFY THE PROBLEM------------------------------
170   GOSUB 2320
180   REM/--------------COME HERE TO SET UP GRID POINTS-----------------------------
190   GOSUB 2560
200   REM/---HERE TO SPECIFY DIFF-COEFFICIENT AND SOUTCE TERM-----------------------
210   GOSUB 2870
220   REM/--------------TO SET UP AE,AW,AP AND DN-----------------------------------
230   GOSUB 520
240   REM/--------NOW TO SOLVE THE ALGEBRAIC EQUATIONS BY TDMA----------------------
250   GOSUB 910
260   ON LS GOTO 350, 270, 350, 270
270   IF DF <= EP THEN GOTO 350
280   PRINT TAB(10); "The Maximum deviation="; DF; " at IT="; KT
290   FOR I = 1 TO M1
300   TA(I) = TA(I) + OM * (T(I) - TA(I))
310   NEXT I
320   DF = 0
330   KT = KT + 1
340   GOTO 200
350   REM/----------TO PRINT OUT GENERAL RESULTS------------------------------------
360   GOSUB 1140
370   ON LS GOTO 440, 440, 380, 380
380   IF T > TM THEN 440
390   FOR I = 1 TO M1
400   TG(I) = T(I)
410   NEXT I
420   KT = 1
430   IF LS = 3 THEN GOTO 220 ELSE 200
440   REM/----------- TO PRINT OUT SPECIAL RESULTS.IF NOT,LEAVE IT OPEN---------------
450    GOSUB 3200
460    IF NF = KM THEN 490
470    NF = NF + 1
480    GOTO 130
490    BEEP
500    END
510    REM/------------SUBROUT INE TO SET UP AE，AW，AP AND CN----------------------
520    REM/-----COEFFICIENTS OF BOUNDARY POINTS-------
530    IF KI > 1 THEN GOTO 590
540    AP(1) = 1
550    AE(1) = 0
560    AW(1) = 0
570    CN(1) = TI
580    GOTO 630
590    AE(1) = GM(1) / XM(2)
600    AP(1) = AE(1) + BI
610    AW(1) = 0
620    CN(1) = AI
630    IF KE > 1 THEN 690
640    AP(M1) = 1
650    AE(M1) = 0
660    AW(M1) = 0
670    CN(M1) = TE
680    GOTO 730
690    AW(M1) = GM(M1) / XP(M2)
700    AP(M1) = AW(M1) + BE
710    AE(M1) = 0
720    CN(M1) = AE
730    REM/-----COEFFICIENTS OF INTERNAL POINTS---------
740    IF LS = 3 AND T > .5 * DT THEN 830
750    EX = 1
760    IF MD = 3 THEN EX = 2
770    AW(2) = GM(2) / XM(2) * RF(2) ^ EX
780    AE(M2) = GM(M2) / XP(M2) * RF(M1) ^ EX
790    FOR I = 2 TO (M2 - 1)
800    AE(I) = RF(I + 1) ^ EX / (XP(I) / GM(I) + XM(I + 1) / GM(I + 1))
810    AW(I + 1) = AE(I)
820    NEXT I
830    FOR I = 2 TO M2
840    CN(I) = CN(I) * (XF(I + 1) - XF(I)) * R(I) ^ EX
850    AP(I) = AE(I) + AW(I) - AP(I) * (XF(I + 1) - XF(I)) * R(I) ^ EX
860    ON LS GOTO 890, 890, 870, 870
870    CN(I) = CN(I) + RC(I) * (XF(I + 1) - XF(I)) * R(I) ^ EX / DT * TG(I)
880    AP(I) = AP(I) + RC(I) * (XF(I + 1) - XF(I)) * R(I) ^ EX / DT
890    NEXT I
900    RETURN
910    REM/------------SUBROUTINE OF TDMA------------------------------------------
920    REM/--------ELIMINATION--------
930    P(1) = AE(1) / AP(1)
940    Q(1) = CN(1) / AP(1)
950    FOR I = 2 TO M1
960    P(I) = AE(I) / (AP(I) - AW(I) * P(I - 1))
970    Q(I) = (CN(I) + AW(I) * Q(I - 1)) / (AP(I) - AW(I) * P(I - 1))
980    NEXT I
990    REM/-------BACK  SUBSTITUTION-------
1000    T(M1) = Q(M1)
1010    FOR J = 1 TO (M1 - 1)
1020    I = M1 - J
1030    T(I) = P(I) * T(I + 1) + Q(I)
1040    NEXT J
1050    IF LS = 2 OR LS = 4 THEN 1060 ELSE 1130
1060    FOR I = 2 TO M1
1070    IF T(I) <= 9.999999E-21 THEN 1100
1080    DS = ABS(T(I) - TA(I)) / T(I)
1090    GOTO 1110
1100    DS = ABS(T(I) - TA(I))
1110    IF DF < DS THEN DF = DS
1120    NEXT I
1130    RETURN
1140    REM/-------------SUBROUTINE TO PRINTING OUT---------------------------------
1150    ON LS GOTO 1190, 1190, 1160, 1160
1160    M = INT((T + .5 * DT) / (K * DT))
1170    IF M = KN THEN 1180 ELSE 1630
1180    KN = KN + 1
1190    IF KP = 2 THEN 1630
1200    LPRINT
1210    LPRINT
1220    LPRINT TAB(20); HO
1230    LPRINT
1240    ON MD GOTO 1250, 1270, 1290, 1310
1250    LPRINT TAB(25); "Cartisian Coorsinates"
1260    GOTO 1320
1270    LPRINT TAB(24); "Cylindrical Coordinates"
1280    GOTO 1320
1290    LPRINT TAB(27); "Polar Coordinates"
1300    GOTO 1320
1310    LPRINT TAB(23); "Nonuniform Cross Section"
1320    ON LS GOTO 1330, 1350, 1380, 1400
1330    LPRINT TAB(24); "Linear Steady Problem"
1340    GOTO 1430
1350    LPRINT TAB(23); "Nonl i near Steady Problem"
1360    LPRINT TAB(25); "Iterative Times="; KT
1370    GOTO 1430
1380    LPRINT TAB(23); "Linear Unsteady Problem"
1390    GOTO 1420
1400    LPRINT TAB(20); "Nonl i near Unsteady Problem"
1410    LPRINT TAB(25); "Iterative Times="; KT
1420    LPRINT TAB(29); "At Time="; T
1430    JE = 0
1440    JB = JE + 1
1450    JE = JE + 4
1460    IF JE > H1 THEN JF = M1
1470    LPRINT
1480    LPRINT TAB(2); "J",
1490    FOR J = JB TO JE
1500    LPRINT J,
1510    NEXT J
1520    IF MD = 2 OR MD = 3 THEN A$ = "R" ELSE A$ = "X"
1530    LPRINT TAB(2); A#,
1540    FOR J = JB TO JE
1550    'PRINT X(J),
1560    NEXT J
1570    LPRINT TAB(2); B#,
1580    FOR J = JB TO JE
1590    LPRINT T(J),
1600    NEXT J
1610    LPRINT
1620    IF J < M1 THEN 1440
1630    T = T + DT
1640    RETURN
1650   REM/************************************************************************
2300   REM/************************************************************************
2310   REM/---Fully Developed Turbulent Heat Transfer in a Tvbe--------------------
2320   REM/--------------------SUBROUTINE SPECI------------------------------------
2330   IF NF = 2 THEN GOTO 2470
2340   LS = 2
2350   MD = 2
2360   KE = 1
2370   AI = 0
2380   BI = 0
2390   TE = 0
2400   KI = 2
2410   EP = .0001
2420   OM = .5
2425   B$ = "W"
2430   H$ = "----- Dimensionless Velocity----------------------"
2440   PR = .7
2445   PT = 1!
2450   RE = 120000!
2455   KR = 2
2460   RETURN
2470   LS = 1
2480   H$ = "-----Dimensionless Temperature--------------------"
2485   B$ = "THETA"
2490   RETURN
2560  REM/-----------------------SUBROUTINE GRID-----------------------------------
2570  XI = .1
2580  XE = 11
2590  PW = .25
2600  M1 = 40
2610  M2 = M1 - 1
2620  FOR I = 2 TO M1
2630  XF(I) = ((I - 2) / (M2 - 1)) ^ PW * (XE - XI) + XI
2640  RF(I) = XF(I)
2650  NEXT I
2660  X(1) = XI
2670  X(M1) = XE
2680  FOR I = 2 TO M2
2690  X(I) = .5 * (XF(I + 1) + XF(I))
2700  XP(I) = XF(I + 1) - X(I)
2710  XM(I) = X(I) - XF(I)
2720  R(I) = X(I)
2730  NEXT I
2740  R(1) = X(1)
2750  REM/----TO SPECIFY THE INITIAL FIELD------------
2755  IF KR > 1 THEN 2770
2760  FOR I = 1 TO M1
2762  TA(I) = 1 - X(I) ^ 2
2764  NEXT I
2766  RETURN
2770  OPEN "I", #1, "DATA"
2774  FOR I = 1 TO M1
2776  INPUT #1, TA(I)
2780  NEXT I
2785  CLOSE #1
2790  RETURN
2870  REM/----------------------SUBROUTINE GAMSSR----------------------------------
2880  ON NF GOTO 2890, 3070
2890  WM = 0!
2900  FOR I = 2 TO M2
2910  WM = WM + TA(I) * (XF(I + 1) ^ 2 - XF(I) ^ 2)
2920  NEXT I
2930  DW = (TA(M2) - TA(M1)) / (X(M1) - X(M2)) / WM
2935  SW = SQR(RE / 2) * SQR(DW) / 26
2940  FOR I = 2 TO M2
2950  PN = (1 - X(I)) * SW
2960  IF PN > 50 THEN EN = 0
2970  IF PN < 50 THEN EN = 1 / EXP(PN)
2980  LR = (.14 - .03 * X(I) ^ 2 - .06 * X(I) ^ 4) * (1 - EN)
2990  RX = (X(I + 1) - X(I)) / (X(I) - X(I - 1))
3000  DW = (TA(I + 1) + TA(I) * (RX ^ 2 - 1) - RX ^ 2 * TA(I - 1))
3005  DW = DW / ((RX + 1) * (X(I + 1) - X(I)))
3010  GM(I) = 1 + .5 * RE * LR ^ 2 * ABS(DW) / WM
3020  CN(I) = 1
3025  AP(I) = 0
3030  NEXT I
3040  GM(I) = 1
3050  GM(M1) = 1
3060  RETURN
3070  FOR I = 2 TO M2
3075  GM(I) = 1 + PR / PT * TS(I)
3080  CN(I) = TR(I) / (3.14159 * WM)
3085  AP(I) = 0
3090  NEXT I
3100  GM(1) = 1
3110  GM(M1) = 1
3120  RETURN
3200  REM/---------------------SUBROUTINE SPRINT------------------------------------
3210  ON NF GOTO 3220, 3250
3220  FOR I = 1 TO M1
3225  TR(I) = T(I)
3230  TS(I) = GM(I) - 1
3235  NEXT I
3240  FR = 8 / WM
3245  RETURN
3250  LPRINT
3255  LPRINT TAB(17); "------Turbulent Diffusion Coefficients------------------"
3260  LPRINT
3262  JE = 0
3264  JB = JE + 1
3266   JE = JE + 4
3268   IF JE > M1 THEN JE = M1
3270   LPRINT TAB(2); "J",
3272   FOR J = JB TO JE
3274   LPRINT J,
3278   NEXT J
3280   LPRINT TAB(2); "Mt/M1",
3282   FOR J = JB TO JE
3284   LPRINT TS(J),
3286   NEXT J
3290   LPRINT TAB(2); "Kt/K1",
3292   FOR J = JB TO JE
3294   LPRINT GM(J) - 1,
3296   NEXT J
3298   LPRINT
3300   LPRINT
3302   IF JE < M1 THEN GOTO 3264
3330   TA = 0
3340   FOR I = 2 TO M2
3342   TA = TA + T(I) * (XF(I + 1) ^ 2 - XF(I) ^ 2) * TR(I)
3344   NEXT I
3346   TA = TA / WM
3348   NU = 1 / (3.14159 * TA)
3350   FF = 1 / (1.82 * LOG(RE) / 2.30259 - 1.64) ^ 2
3352   NP = (FF / 8) * RE * PR / (1.07 + 12.7 * SQR(FF / 8) * (PR ^ .6666 - 1))
3354   DF = (FR / RE - FF) / FF * 100
3356   DN = (NU - NP) / NP * 100
3358   LPRINT
3360   LPRINT TAB(20); "F="; FR / RE; ", for Re="; RE; ", from computation"
3362   LPRINT TAB(20); "Nu="; NU; ", for Re="; RE; ", Pr="; PR; " and Prt="; PT
3364   LPRINT TAB(20); "F="; FF; " from Filonenko  equation"
3366   LPRINT TAB(20); "Nu="; NP; " from Petukhov equation"
3368   LPRINT TAB(20); "Deviation of F="; DF; "%"
3370   LPRINT TAB(20); "Deviation of Nu="; DN; "%"
3375   LPRINT TAB(20); "Epsilon="; EP; ", Omeqa="; QM; ",PW="; PW
3380   OPEN "O", #1, "DATA"
3382   FOR I = 1 TO M1
3384   WRITE #1, TR(I)
3386   NEXT I
3388   RETURN
3390   REM/**************************************************************************/