10  OPTION BASE 1
20  DIM an(4), ma(4, 4), main(4, 4), bn(4), u(4, 8), v(4)
30  REM             data input
40  PRINT "Please input all angles in degrees"
50  PRINT "Please input Wing Span, Root Chord, and Tip Chord"
60  INPUT b, cr, ct
70  PRINT b, cr, ct
80  s = (cr + ct) * b / 2
90  ar = b * b / s
100 PRINT "Please input tip twist relative to root (nose up positive)"
110 INPUT tw
120 PRINT "Please input the aerofoil section lift curve slope (per rad)"
130 INPUT a
140 PRINT "Please input the aerofoil section zero lift angle"
150 INPUT a10
160 PRINT "Please input the root section incidence"
170 INPUT alr
180 REM    Calculates the coefficients of the fourier exspansions
190 REM    at phi=22.5,45,67.5 and 90 degrees, and prints them
200 PRINT
210 PRINT "The array of coefficients is:-"
220 e = 22.5: pi = 3.141527
230 FOR i = 1 TO 4
240     fi = i * e * pi / 180
250     al = (alr + tw * COS(fi) - a10) * pi / 180
260     ch = cr - (cr - ct) * COS(fi)
270     bn(i) = ch * a * al / (4 * b) * SIN(fi)
280         FOR j = 1 TO 4
290         ma(i, j) = SIN((2 * j - 1) * fi) * (SIN(fi) + (2 * j - 1) * ch * a / (4 * b))
310         NEXT j
315         PRINT USING "##.####  ##.####  ##.####  ##.####"; ma(i, 1); ma(i, 2); ma(i, 3); ma(i, 4)
330     NEXT i
335     PRINT
340 REM Finds main, the inverse matrix of ma, and prints it
350 d2 = 1
360 FOR i = 1 TO 4
370  FOR j = 1 TO 4
380  u(i, j) = ma(i, j)
390  IF i = j THEN 420
400  u(i, j + 4) = 0
410  GOTO 430
420  u(i, j + 4) = 1
430  NEXT j
440  v(i) = i
450 NEXT i
470 FOR i = 1 TO 4
480 m = 0
490  FOR j = 1 TO 4
500  IF m > ABS(u(i, j)) THEN 530
510  m = ABS(u(i, j))
520  r = j
530  NEXT j
540  IF i = 4 THEN 650
550  IF r = i THEN 650
560  d2 = d2 * (-1)
570 FOR j = 1 TO 8
580 t = u(i, j)
590 u(i, j) = u(r, j)
600 u(r, j) = t
610 NEXT j
620 t = v(r)
630 v(r) = v(i)
640 c(i) = t
650 d = u(i, i)
660 d2 = d2 * d
670   FOR j = 1 TO 8
680   u(i, j) = u(i, j) / d
690   NEXT j
700 FOR i1 = 1 TO 4
710 IF i1 = i THEN 760
720 d = u(i1, i)
730    FOR j = i TO 8
740    u(i1, j) = u(i1, j) - d * u(i, j)
750    NEXT j
760     NEXT i1
770      NEXT i
780 PRINT
790 PRINT "The inverse matrix is:"
800 FOR i = 1 TO 4
810  FOR j = 1 TO 4
820  main(i, j) = u(i, j + 4)
840  NEXT j
850  PRINT USING "##.####  ##.####  ##.####  ##.####"; main(i, 1); main(i, 2); main(i, 3); main(i, 4)
860  NEXT i
870 REM Multiply the inverse matrix by the vector bn
880 REM and print the result
890 PRINT
900 PRINT "The fourier coefficients are:-"
910 FOR i = 1 TO 4
920 an(i) = 0
930  FOR j = 1 TO 4
940  an(i) = an(i) + main(i, j) * bn(j)
950  NEXT j
960 PRINT "a"; 2 * i - 1; " = "; an(i)
970 NEXT i
980 PRINT "                      Press return when ready"
990 INPUT a$
1000 REM     Prints out the results for the wing
1010 PRINT "aspect ratio = ", ar
1020 PRINT "taper ratio = ", ct / cr
1030 PRINT "root incidence = ", alr
1040 PRINT "washout =", -tw
1050 PRINT
1060 cl = pi * ar * an(1)
1070 PRINT USING "cl = #.####"; cl
1080 cdi = 0
1090 FOR i = 1 TO 4
1100 cdi = cdi + (2 * i - 1) * an(i) * an(i)
1110 NEXT i
1120 cdi = cdi * pi * ar
1130 PRINT USING "cdi = #.#####"; cdi
1140 PRINT USING "Elliptic loading cdi = #.#####"; cl * cl / (pi * ar)
1150 REM             Outputs local cl across span
1160 PRINT
1170 PRINT "The spanwise cl distribution is:-"
1180 PRINT "  phi        local cl      cos(phi)    local cl/wing cl"
1190  FOR phi = 0 TO 90 STEP 10
1200  fi = phi * pi / 180
1210  lc = 0
1220  FOR i = 1 TO 4
1230  lc = lc + an(i) * SIN((2 * i - 1) * fi)
1240  NEXT i
1250  IF ct = 0 AND phi = 0 THEN 1280
1260  lc = 4 * b * lc / (cr - (cr - ct) * COS(fi))
1270 GOTO 1290
1280 lc = 99.9999
1290 PRINT USING "  ##         ##.####         #.##        ###.####"; phi; lc; COS(fi); lc / cl
1300 NEXT phi
1310 END