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