10   OPTION BASE 1
20   PRINT CHR$(12)
40   PRINT "Please input ROOT CHORD and TIP CHORD as fractions of the span"
50   INPUT cr, ct
60   PRINT cr, ct
70   PRINT "Please input the SWEEPBACK ANGLE in degrees."
80   INPUT l
90   PRINT l
95   pi = 3.141593
100  m = 2
110  n = 6
120  tot = m * n
130  n1 = tot / 2
140  DIM xl(12), yl(12), xr(12), yr(12), xc(12), yc(12), xroot(12), yroot(12)
150  FOR i = 1 TO m
160      xtip(i) = .5 * TAN(l * pi / 180) + ct / m * (i - .75)
170      ytip(i) = -.5
180      xroot(i) = cr / m * (i - .75)
190      yroot(i) = 0
200   FOR j = 1 TO n / 2
210   k = (i - 1) * n / 2 + j
220   k1 = tot / 2 + n / 2 * i + 1 - j
230   xl(k) = xtip(i) - (xtip(i) - xroot(i)) / (n / 2) * (j - 1)
240   yl(k) = ytip(i) + (yroot(i) - ytip(i)) / (n / 2) * (j - 1)
250   xr(k) = xl(k) - (xtip(i) - xroot(i)) / (n / 2)
260   yr(k) = yl(k) + (yroot(i) - ytip(i)) / (n / 2)
270   xr(k1) = xl(k)
280   yr(k1) = -yl(k)
290   xl(k1) = xr(k)
300   yl(k1) = -yr(k)
310   NEXT j
320   NEXT i
330   PRINT "Panel Numbers and Coordinates of the bound vortex end points."
350   FOR k = 1 TO tot
360   PRINT USING "##     ###.###     ###.###     ###.###     ###.###"; k; xl(k); yl(k); xr(k); yr(k)
370   NEXT k
380   GOSUB 1710
390   FOR i = 1 TO m
400   xtip(i) = xtip(i) + ct / (2 * m)
410   xroot(i) = xroot(i) + cr / (2 * m)
420   FOR j = 1 TO n / 2
430   k = (i - 1) * n / 2 + j
440   k1 = tot / 2 + i * n / 2 + 1 - j
450   xc(k) = xtip(i) - (xtip(i) - xroot(i)) / (n / 2) * (j - .5)
460   yc(k) = (yl(k) + yr(k)) / 2
470   xc(k1) = xc(k)
480   yc(k1) = -yc(k)
490   NEXT j
500   NEXT i
510   PRINT "      Control Point Coordinates"
530   FOR k = 1 TO tot
540   PRINT USING "##     ###.###     ###.###"; k; xc(k); yc(k)
550   NEXT k
560   GOSUB 1710
570   PRINT "Calculations from this point on are for the left wing only!!!"
580   PRINT "Values for the right wing may be obtained using wing symmetry."
590   GOSUB 1710
600   PRINT "The matrix of computed velocities induced at each control point"
610   PRINT "on the left wing by the horseshoe vorticies covering the whole"
620   PRINT "wing is proportional to:"
630   PRINT
640   DIM ind(6, 6)
650   FOR k = 1 TO tot / 2
660   FOR kk = 1 TO tot / 2
670   ind(k, kk) = 0
680   NEXT kk
690   NEXT k
700   FOR k = 1 TO tot / 2
710   FOR kk = 1 TO tot
720   d1 = SQR((xl(kk) - xc(k)) ^ 2 + (yl(kk) - yc(k)) ^ 2)
730   d2 = SQR((xr(kk) - xc(k)) ^ 2 + (yr(kk) - yc(k)) ^ 2)
740   d3 = (xr(kk) - xc(k)) * (yl(kk) - yc(k)) - (xl(kk) - xc(k)) * (yr(kk) - yc(k))
750   d4 = yl(kk) - yc(k)
760   d5 = yr(kk) - yc(k)
770   vbc = ((xr(kk) - xl(kk)) * (xl(kk) - xc(k)) + (yr(kk) - yl(kk)) * (yl(kk) - yc(k))) / d1
780   vbc = vbc - ((xr(kk) - xl(kk)) * (xr(kk) - xc(k)) + (yr(kk) - yl(kk)) * (yr(kk) - yc(k))) / d2
790   vbc = vbc / d3
800   vab = 1 - (xl(kk) - xc(k)) / d1
810   vab = vab / d4
820   vcd = 1 - (xr(kk) - xc(k)) / d2
830   vcd = -vcd / d5
840   w = vab + vbc + vcd
850   IF kk > tot / 2 THEN 880
860   ind(k, kk) = w
870   GOTO 910
880   j = INT((kk - tot / 2 - .5) / (n / 2)) + 1
890   kkk = (2 * j - 1) * (n / 2) + 1 - (kk - tot / 2)
900   ind(k, kkk) = ind(k, kkk) + w
910   NEXT kk
920   PRINT USING "#####.#####"; ind(k, 1); ind(k, 2); ind(k, 3); ind(k, 4); ind(k, 5); ind(k, 6)
940   NEXT k
970   DIM a(6, 12), c(6)
980   d2 = 1
990   FOR i = 1 TO n1
1000  FOR j = 1 TO n1
1010  a(i, j) = ind(i, j)
1020  IF i = j THEN 1050
1030  a(i, j + n1) = 0
1040  GOTO 1060
1050  a(i, j + n1) = 1
1060  NEXT j
1070  c(i) = i
1080  NEXT i
1090  FOR i = 1 TO n1
1100  p = 0
1110  FOR j = 1 TO n1
1120  IF p > ABS(a(i, j)) THEN 1150
1130  p = ABS(a(i, j))
1140  r = j
1150  NEXT j
1160  IF i = n1 THEN 1270
1170  IF r = i1 THEN 1250
1180  d2 = -d2
1190  FOR j = 1 TO 2 * n1
1200  t = a(i, j)
1210  a(i, j) = a(r, j)
1220  a(r, j) = t
1230  NEXT j
1240  t = c(r)
1250  c(r) = c(i)
1260  c(i) = t
1270  d = a(i, i)
1280  d2 = d2 * d
1290  FOR j = i TO 2 * n1
1300  a(i, j) = a(i, j) / d
1310  NEXT j
1320  FOR ii = 1 TO n1
1330  IF ii = i THEN 1380
1340  d = a(ii, i)
1350  FOR j = i TO 2 * n1
1360  a(ii, j) = a(ii, j) - d * a(i, j)
1370  NEXT j
1380  NEXT ii
1390  NEXT i
1400  GOSUB 1710
1410  PRINT "The inverse matrix is:"
1420  FOR i = 1 TO n1
1430  FOR j = 1 TO n1
1440  a(i, j) = a(i, j + n1)
1445  NEXT j
1450  PRINT USING "#####.#####"; a(i, 1); a(i, 2); a(i, 3); a(i, 4); a(i, 5); a(i, 6)
1480  NEXT i
1490  GOSUB 1710
1500  gg = 0
1510  PRINT "Vortex Strengths:"
1520  FOR i = 1 TO n1
1530  g = 0
1540  FOR j = 1 TO n1
1550  g = g + a(i, j)
1560  NEXT j
1570  PRINT USING "g(#)/4.pi.b.u.alpha)   #.#####"; i; -g
1580  gg = gg - g
1590  NEXT i
1600  GOSUB 1710
1610  lcs = 32 * pi * gg / ((cr + ct) * n)
1630  PRINT USING "Wing Lift Curve Slope               ##.###"; lcs
1640  PRINT USING "Wind Area                           ##.###"; (cr + ct) / 2
1650  PRINT USING "Aspect Ratio                        ##.###"; 2 / (cr + ct)
1660  PRINT USING "Taper Ratio                         ##.###"; ct / cr
1670  PRINT USING "Leading Edge Sweepback Angle        ##.###"; l
1680  PRINT USING "Trailing Edge Sweepback Angle       ##.###"; (ATN((.5 * TAN(l * pi / 180) + ct - cr) / .5)) * 180 / pi
1690  PRINT USING "Quarter Chord Line Sweepback Angle  ##.###"; (ATN((.5 * TAN(l * pi / 180) + (ct - cr) / 4) / .5)) * 180 / pi
1700  GOTO 1750
1710  PRINT
1720  PRINT "                                                      Press RETURN"
1730  INPUT a$
1740  RETURN
1750  END