REM Copyright (C) 2001, John J. G. Savard

DIM im%(40, 640), ma%(16)
DIM ni(20), ns(20), si(20), ss(20)
DIM r1(5, 30), r2(5, 30)

REM DEF FNQ (D) = SGN(D) * (ABS(D) ^ .3333333)

FOR i = 0 TO 15
READ ma%(i)
NEXT i

FOR i = 1 TO 5
ii = 1
2 READ r1(i, ii), r2(i, ii)
IF r1(i, ii) < 99 THEN ii = ii + 1: GOTO 2
NEXT i

REM start of code to set pixels in image to 1

REM Draw a map of the world

pi# = 3.141592653589793#
dr# = pi# / 180#
rd# = 180# / pi#

3 PRINT "Select projection (0 to finish):"
REM PRINT "Cylindrical projections:"
PRINT " 1) Mercator              2) Gall                    3) Miller"
PRINT " 4) Plate Carree          5) Cylindrical Equal-Area"
REM PRINT "Azimuthal projections (in hemispheres unless otherwise noted):"
PRINT " 6) Gnomonic (incomplete) 7) Stereographic           8) Azimuthal Equidistant"
PRINT " 9) Lambert Azimuthal Equal-area                    10) Orthographic"
PRINT "11) Azimuthal Equidistant (single circle)"
REM PRINT "Conic projections:"
PRINT "12) Lambert conformal    13) Simple conic           14) Albers equal-area"
REM PRINT "Pseudocylindrical and pseudoconic projections:"
PRINT "15) Sinusoidal           16) Bonne's                17) Mollweide's"
PRINT "18) Eckert IV            19) MacBryde's Flat Polar Quartic"
PRINT "20) Adams' Equal-Area    21) Parabolic stereographic"
PRINT "22) Loximuthal"
REM PRINT "Polyconic projections:"
PRINT "23) Polyconic"
REM PRINT "Conventional projections:"
PRINT "24) Globular             25) Van der Grinten's      26) Van der Grinten IV"
PRINT "27) Savard Egg           28) Ster Modified VdG IV   29) Globular (whole world)"
REM PRINT "Other conformal projections:"
PRINT "30) Lagrange's conformal 31) August's conformal     32) Cosine-based"
PRINT "33) Eisenlohr            34) Guyou's doubly-periodic"
PRINT "35) World in rectangle   36) World in diamond"
REM PRINT "Other equal-area projections:"
PRINT "37) Hammer-Aitoff        38) Generalized Aitoff Equal-Area"
REM PRINT "Other projections:"
PRINT "39) Aitoff-Wagner        40) Aitoff projection    41) Winkel's Tripel"
PRINT "42) Pole-line Aitoff     43) Rectangular Polyconic"
PRINT "44) Four-square world    45) Quadfoil 1           46) New Conventional"
PRINT "47) Conventional II      48) Modified VdG IV      49) Squashed Conformal"
PRINT "50) Ginzoid              51) Larrivee             52) Laskowski Tri-Optimal"
PRINT "53) Generalized Aitoff   54) Approximate Canters"
PRINT "55) Canters' Minimum-Error Polyconic              56) Canters' Equal-Scale"
INPUT pj%

IF pj% = 0 THEN 990

px% = (pj% - 1) \ 10
py% = pj% - 10 * px%
px% = px% + 1

LINE INPUT "Output file name? "; fi$
OPEN fi$ FOR OUTPUT AS #2
PRINT "Maximum image size: 640 by 640."
INPUT "Image width: ", w%
INPUT "Image height: ", h%

wb% = w% \ 16
hw% = w% \ 2
hh% = h% \ 2

REM PRINT "Clearing image array."

FOR j% = 1 TO h%
FOR i% = 1 TO wb%
im%(i%, j%) = 0
NEXT i%
NEXT j%

PRINT "Tilt by 180 degrees for conics with SP in southern hemisphere"
PRINT "Transverse rotation not required for any orientation, but allows"
PRINT "for ease in matching projections without solving spherical triangles."
PRINT
INPUT "Meridian on globe along which to tilt pole: ", a1
INPUT "Amount of tilt: ", a2
INPUT "Transverse clockwise rotation: ", a4
INPUT "Meridian on map to assign meridian to: ", a3

ti% = 1
IF (ABS(a2) < .001) AND (ABS(a4) < .001) THEN ti% = 0: a3 = a3 - a1: a1 = 0: GOTO 4
tc = COS(dr# * a2)
ts = SIN(dr# * a2)

IF ABS(a4) < .001 THEN 4
ec = COS(dr# * a4)
es = SIN(dr# * a4)

REM initialization for projection
4 ip% = 0
pa% = (pj% - 1) \ 20
pb% = pj% - 20 * pa%
pa% = pa% + 1

ON pa% GOTO 5, 6, 7

5 ON pb% GOTO 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200
6 ON pb% GOTO 210, 220, 230, 240, 250, 260, 270, 280, 290, 300, 310, 320, 330, 340, 350, 360, 370, 380, 390, 400
7 ON pb% GOTO 410, 420, 430, 440, 450, 460, 470, 480, 490, 500, 510, 520, 530, 540, 550, 560

REM Assign topology type, border, aspect ratio (width/height)

REM Topolgies:
REM  1) Global azimuthal; check for too large a line at far pole
REM  2) Cylindrical type; check for crossing date line
REM  3) Hemispherical azimuthal; check for crossing equator
REM  4) Hemispherical conventional; check for crossing date line/prime meridian
REM  5) Polyconic gore; check for switching gores
REM  6) Conic
REM  7) Gnomonic cube
REM  8) Interrupted projection
REM  9) Guyou's
REM 10) Curved pole-line; cylindrical and too large a line at either pole
REM 11) Check for crossing 0, 90, 180, -90 in southern hemisphere

REM Borders:
REM  1) Circle/Ellipse (shape depends on aspect ratio)
REM  2) Two circles
REM  3) Rectangle (shape depends on aspect ratio)
REM  4) Conic
REM  5) Eckert IV type
REM  6) Sinusoidal, August's, or other projection-native
REM  7) Quadrifoil Conformal
REM  8) Projection-native with linear poles (i.e., flat polar quartic)
REM  9) Projection-native with curved poles (i.e., certain Aitoff)
REM 10) Guyou's, Gnomonic near-hemispheres
REM 11) Diamond
REM 12)
REM 13) Conventional II

10 tp% = 2: bo% = 3
INPUT "North latitude limit: ", ql
INPUT "South latitude limit: ", qs
qn = rd# * LOG((1 / COS(dr# * ql)) + TAN(dr# * ql))
qv = rd# * LOG((1 / COS(dr# * qs)) + TAN(dr# * qs))
ar = 360 / (qn + qv)
qv = .5 * (qv - qn)
GOTO 665
20 tp% = 2: bo% = 3
k1 = rd# * (1 + SQR(2))
ar = 180 / k1
GOTO 665
30 tp% = 2: bo% = 3: ar = 1.363886
GOTO 665
40 tp% = 2: bo% = 3: ar = 2
GOTO 665
50 tp% = 2: bo% = 3: ar = 3.141592
GOTO 665
60 tp% = 3: bo% = 10: ar = 2
INPUT "Angular limit from center: ", ga
gt = TAN(dr# * ga)
gb = 90 - rd# * ATN(1.414213 * gt)
gs = 90 / gt
GOTO 665
70 tp% = 3: bo% = 2: ar = 2
GOTO 665
80 tp% = 3: bo% = 2: ar = 2
GOTO 665
90 tp% = 3: bo% = 2: ar = 2
GOTO 665
100 tp% = 3: bo% = 2: ar = 2
GOTO 665
110 tp% = 1: bo% = 1: ar = 1
GOTO 665
120 INPUT "Standard parallel: ", sp
k1 = SIN(dr# * sp)
INPUT "Southern limit of latitude: ", qs
k2 = 180 / (TAN(.5 * dr# * (90 + qs)) ^ k1)
tp% = 2: bo% = 12: ar = 2
GOTO 665
130 INPUT "Standard parallel: ", sp
tp% = 2: bo% = 9: ar = 2
k1 = SIN(dr# * sp)
k2 = rd# / TAN(dr# * sp)
GOTO 665
140 INPUT "Standard parallel: ", sp
tp% = 2: bo% = 9: ar = 2
k1 = SIN(dr# * sp)
k2 = COS(dr# * sp)
k3 = k1 + .5 * k2 / k1
r1 = SQR(k3 - 1)
r2 = SQR(k3 + 1)
k4 = 180 / (r2 - r1)
k5 = .5 * k4 * (r2 + r1)
GOTO 665
150 sp% = 0: GOSUB 3900
tp% = 8: bo% = 6: ar = 2
GOTO 665
160 PRINT "(If integer, interruptions allowed)"
INPUT "Standard parallel: ", sp
sp% = sp
IF sp% = sp THEN 165
sp% = 0: ip% = 1
ni(0) = -180: ni(1) = 180: ns(1) = 0
si(0) = -180: si(1) = 180: ss(1) = 0
GOTO 167
165 GOSUB 3900
167 tp% = 8: bo% = 6: ar = 2
REM k1 = rd# / TAN(dr# * sp) valid, except for sp=90, when
REM cos is zero, but sin is 1; changed form allows Werner's
REM projection
k1 = rd# * COS(dr# * sp) / SIN(dr# * sp)
k2 = dr# * SIN(dr# * sp)
GOTO 665
170 sp% = 0: GOSUB 3900
tp% = 8: bo% = 6: ar = 2
REM Border was 1 when not interrupted
GOTO 665
180 sp% = 0: GOSUB 3900
tp% = 8: bo% = 8: ar = 2
REM Border was 5 when not interrupted
GOTO 665
190 sp% = 0: GOSUB 3900
tp% = 8: bo% = 8: ar = 2.221442
GOTO 665
200 sp% = 0: GOSUB 3900
tp% = 8: bo% = 6: ar = 2.221442
GOTO 665
210 sp% = 0: GOSUB 3900
tp% = 8: bo% = 6: ar = 1.570796
GOTO 665
220 INPUT "Standard parallel? ", sp
k1 = rd# * LOG((1 / COS(dr# * sp)) + TAN(dr# * sp))
tp% = 2: bo% = 6: ar = 2
GOTO 665
230 GOSUB 3800
tp% = 8: bo% = 6: ar = 2
GOTO 665
240 tp% = 4: bo% = 2: ar = 2
GOTO 665
250 tp% = 2: bo% = 1: ar = 1
GOTO 665
260 tp% = 2: bo% = 6: ar = 1.6
GOTO 665
270 tp% = 2: bo% = 1: ar = 1.5
GOTO 665
280 tp% = 2: bo% = 6: ar = 1.63
GOTO 665
290 tp% = 2: bo% = 6: ar = 1.6
GOTO 665
300 tp% = 2: bo% = 1: ar = 1
GOTO 665
310 tp% = 2: bo% = 6: ar = 1.42
GOTO 665
320 tp% = 2: bo% = 6: ar = 1.5
GOTO 665
330 tp% = 2: bo% = 6: ar = 1.5
GOTO 665
340 tp% = 9: bo% = 10: ar = 2
GOTO 665
350 INPUT "Parameter (m, k^2, (sin(alpha))^2) (.5: square): ", k1
w = .5 * pi#: m = k1: GOSUB 3100: wi = w
w = .5 * pi#: m = 1 - k1: GOSUB 3100: ht = w
tp% = 2: bo% = 3: ar = wi / ht
k2 = 180 / wi
PRINT "Aspect ratio will be: "; ar
GOTO 665
360 tp% = 2: bo% = 11: ar = 1
GOTO 665
370 tp% = 2: bo% = 1: ar = 2
GOTO 665
380 PRINT "Compression factor (or its reciprocal)"
INPUT ".9: original Aitoff, 1.5: 3:2 ellipse: ", p
IF p < 1 THEN sh = p: ex = 1 / p: GOTO 385
sh = 1 / p: ex = p
385 tp% = 10: th% = 16: bo% = 9: ar = 2 * sh
GOTO 665
390 tp% = 2: bo% = 9: ar = 1.6
GOTO 665
400 tp% = 2: bo% = 1: ar = 2
GOTO 665
410 tp% = 2: bo% = 8: ar = 1.766044
GOTO 665
420 INPUT "Compression factor (or its reciprocal): ", p
IF p < 1 THEN sh = p: ex = 1 / p: GOTO 425
sh = 1 / p: ex = p
425 tp% = 10: th% = 16: bo% = 9: ar = 2 * sh
GOTO 665
430 INPUT "Standard parallel: ", sp
GOSUB 3800
k1 = SIN(dr# * sp)
tp% = 8: bo% = 6: ar = 2
GOTO 665
440 tp% = 11: bo% = 3: ar = 1
GOTO 665
450 tp% = 11: bo% = 7: ar = 1
GOTO 665
460 tp% = 10: th% = 15: bo% = 13: ar = 1.333333
GOTO 665
470 tp% = 10: th% = 15: bo% = 13: ar = 1.333333
GOTO 665
480 INPUT "Projection distance (180 for standard VdG IV): ", D
k1 = 2 * D
k2 = D + 90
k3 = D * D
k4 = 4 * k3
tp% = 2: bo% = 6: ar = 1.6
GOTO 665
490 INPUT "Number of stages (2-5): ", k1
ON k1 - 1 GOTO 492, 493, 494, 495
492 ar = 1.311: GOTO 499
493 ar = 1.454: GOTO 499
494 ar = 1.403: GOTO 499
495 ar = 1.37: GOTO 499
499 tp% = 2: bo% = 6
GOTO 665
500 GOSUB 3800
INPUT "Overall size scaling factor: ", k1
IF k1 < 1 THEN k1 = 1 / k1
INPUT "Vertical radius scaling factor: ", k2
IF k2 < 1 THEN k2 = 1 / k2
INPUT "Longitude size scaling factor: ", k3
IF k3 < 1 THEN k3 = 1 / k3
INPUT "Constant admixture to longitude: ", k4
IF k4 > 1 THEN k4 = 1 / k4
INPUT "Vertical stretch: ", k5
tp% = 8: th% = 8: bo% = 9: ar = 1.8
GOTO 665
510 tp% = 10: th% = 8: bo% = 9: ar = 1.3
GOTO 665
520 tp% = 10: th% = 24: bo% = 9: ar = 1.8
GOTO 665
530 INPUT "Latitude compression factor (or its reciprocal): ", p
IF p < 1 THEN k1 = p: k2 = 1 / p: GOTO 533
k1 = 1 / p: k2 = p
533 INPUT "Longitude compression factor (or its reciprocal): ", p
IF p < 1 THEN k3 = p: k4 = 1 / p: GOTO 535
k3 = 1 / p: k4 = p
535 INPUT "Vertical stretch: ", p
k5 = p: IF p < 1 THEN k5 = 1 / p
tp% = 10: th% = 16: bo% = 9: ar = 2 * k1 / k5
GOTO 665
540 tp% = 10: th% = 16: bo% = 9: ar = 1.5
GOTO 665
550 tp% = 10: th% = 16: bo% = 9: ar = 1.5
GOTO 665
560 tp% = 10: th% = 16: bo% = 9: ar = 1.5
GOTO 665

665 vd = 0: hd = 0
INPUT "Zoom factor? ", zm
w2% = h% * ar
IF w% < w2% THEN w2% = w%
sf = zm * (w2% - 8) / 360

IF zm < 1.01 THEN 670
INPUT "Center (x,y): ", hd, vd
hd = -hd * zm
vd = -vd * zm

rt% = 0
INPUT "Clockwise rotation of final map: ", a
IF ABS(a) < .001 THEN 670
rt% = 1
rc = COS(dr# * a)
rs = SIN(dr# * a)

670 PRINT "1: line map 4: new line maps 8: parallels 16: meridians 32: point map"
PRINT "64: special parallels"
INPUT "Map only (1), Graticule only (2), or map and graticule(3)? ", ac%

IF (ac% AND 1) = 0 THEN 675
LINE INPUT "Input file name? "; fj$
OPEN fj$ FOR INPUT AS #1

675 IF (ac% AND 2) = 0 THEN 680

ac% = ac% OR 24

INPUT "Graticule increment: ", gi
gx = gi: gy = gi
GOTO 690

680 IF (ac% AND 16) = 0 THEN 685
INPUT "Longitude increment for meridians: ", gx

685 IF (ac% AND 8) = 0 THEN 690
INPUT "Latitude increment for parallels: ", gy

690 IF (ac% AND 64) = 0 THEN 695
INPUT "Numerical increment: "; gn
INPUT "Limit: "; gt

695 IF (ac% AND 32) = 0 THEN 700
LINE INPUT "Point data input file? "; fp$
OPEN fp$ FOR INPUT AS #3

REM draw map border

700 pv% = 0

REM PRINT "Drawing border of type "; bo%

IF bo% = 0 THEN 770
ON bo% GOTO 705, 710, 715, 720, 725, 730, 735, 740, 745, 750, 755, 760, 765

705 FOR i = 0 TO 360
ag = dr# * i
x2 = 180 * COS(ag)
y2 = 180 * SIN(ag) / ar
GOSUB 4000
NEXT i
GOTO 770

710 FOR i = 0 TO 360
ag = dr# * i
x2 = (90 * COS(ag) - 90)
y2 = 180 * SIN(ag) / ar
GOSUB 4000
NEXT i
pv% = 0
FOR i = 0 TO 360
ag = dr# * i
x2 = (90 * COS(ag) + 90)
y2 = 180 * SIN(ag) / ar
GOSUB 4000
NEXT i
GOTO 770

715 xb = 180: yb = 180 / ar
x2 = xb: y2 = yb
GOSUB 4000
x2 = xb: y2 = -yb
GOSUB 4000
x2 = -xb: y2 = -yb
GOSUB 4000
x2 = -xb: y2 = yb
GOSUB 4000
x2 = xb: y2 = yb
GOSUB 4000
GOTO 770

720 REM

725 FOR i = 90 TO 270
ag = dr# * i
x2 = 90 * COS(ag) - 90
y2 = 180 * SIN(ag) / ar
GOSUB 4000
NEXT i
FOR i = -90 TO 90
ag = dr# * i
x2 = 90 * COS(ag) + 90
y2 = 180 * SIN(ag) / ar
GOSUB 4000
NEXT i
x2 = -90
y2 = 180 / ar
GOSUB 4000
GOTO 770

730 IF ip% = 1 THEN 731
ix = -180
FOR iy = -90 TO 90
aa = iy: ao = ix
GOSUB 1080
NEXT iy
ix = 180
FOR iy = 90 TO -90 STEP -1
aa = iy: ao = ix
GOSUB 1080
NEXT iy
GOTO 770

731 ix = -180
ao = ix
vo = ix - ss(1)
dx = ss(1)
FOR iy = -90 TO sp%
aa = iy
GOSUB 1900
NEXT iy
ao = ix
vo = ix - ns(1)
dx = ns(1)
FOR iy = sp% + 1 TO 90
aa = iy
GOSUB 1900
NEXT iy

REM Northern hemisphere
ic% = 1
732 ix = ni(ic%)
ao = ix
vo = ix - ns(ic%)
dx = ns(ic%)
FOR iy = 89 TO sp% STEP -1
aa = iy
GOSUB 1900
NEXT iy

IF ni(ic%) > 179.9 THEN ic% = 1: pv% = 0: GOTO 733

ao = ix
vo = ix - ns(ic% + 1)
dx = ns(ic% + 1)
FOR iy = sp% + 1 TO 90
aa = iy
GOSUB 1900
NEXT iy

ic% = ic% + 1: GOTO 732

REM Southern hemisphere

733 ix = si(ic%)
ao = ix
vo = ix - ss(ic%)
dx = ss(ic%)
FOR iy = -90 TO sp%
aa = iy
GOSUB 1900
NEXT iy

IF si(ic%) > 179.9 THEN 770

ao = ix
vo = ix - ss(ic% + 1)
dx = ss(ic% + 1)
FOR iy = sp% - 1 TO -89 STEP -1
aa = iy
GOSUB 1900
NEXT iy

ic% = ic% + 1: GOTO 733

735 FOR iy = -100 TO 100
x = -1
y = .01 * iy
x2 = 20 * (5 * x - x * x * x * x * x + 10 * x * x * x * y * y - 5 * x * y * y * y * y)
y2 = 20 * (5 * y - 5 * x * x * x * x * y + 10 * x * x * y * y * y - y * y * y * y * y)
GOSUB 4000
NEXT iy
FOR ix = -99 TO 100
x = .01 * ix
y = 1
x2 = 20 * (5 * x - x * x * x * x * x + 10 * x * x * x * y * y - 5 * x * y * y * y * y)
y2 = 20 * (5 * y - 5 * x * x * x * x * y + 10 * x * x * y * y * y - y * y * y * y * y)
GOSUB 4000
NEXT ix
FOR iy = 99 TO -100 STEP -1
x = 1
y = .01 * iy
x2 = 20 * (5 * x - x * x * x * x * x + 10 * x * x * x * y * y - 5 * x * y * y * y * y)
y2 = 20 * (5 * y - 5 * x * x * x * x * y + 10 * x * x * y * y * y - y * y * y * y * y)
GOSUB 4000
NEXT iy
FOR ix = 99 TO -100 STEP -1
x = .01 * ix
y = -1
x2 = 20 * (5 * x - x * x * x * x * x + 10 * x * x * x * y * y - 5 * x * y * y * y * y)
y2 = 20 * (5 * y - 5 * x * x * x * x * y + 10 * x * x * y * y * y - y * y * y * y * y)
GOSUB 4000
NEXT ix
GOTO 770

740 IF ip% = 1 THEN 741
ix = -180
FOR iy = -90 TO 90
aa = iy: ao = ix
GOSUB 1080
NEXT iy
ix = 180
FOR iy = 90 TO -90 STEP -1
aa = iy: ao = ix
GOSUB 1080
NEXT iy
ao = -180: aa = -90
GOSUB 1080
GOTO 770

741 ix = -180
ao = ix
vo = ix - ss(1)
dx = ss(1)
FOR iy = -90 TO 0
aa = iy
GOSUB 1900
NEXT iy
ao = ix
vo = ix - ns(1)
dx = ns(1)
FOR iy = 1 TO 90
aa = iy
GOSUB 1900
NEXT iy

REM Northern hemisphere
ic% = 1
742 ix = ni(ic%)
ao = ix
vo = ix - ns(ic%)
dx = ns(ic%)
FOR iy = 90 TO 0 STEP -1
aa = iy
GOSUB 1900
NEXT iy

IF ni(ic%) > 179.9 THEN 743

ao = ix
vo = ix - ns(ic% + 1)
dx = ns(ic% + 1)
FOR iy = 0 TO 90
aa = iy
GOSUB 1900
NEXT iy

ic% = ic% + 1: GOTO 742

REM Southern hemisphere

743 ic% = 1
pv% = 0

ix = -180
ao = ix
vo = ix - ss(1)
dx = ss(1)
iy = -90
aa = iy
GOSUB 1900

744 ix = si(ic%)
ao = ix
vo = ix - ss(ic%)
dx = ss(ic%)
FOR iy = -90 TO 0
aa = iy
GOSUB 1900
NEXT iy

IF si(ic%) > 179.9 THEN 770

ao = ix
vo = ix - ss(ic% + 1)
dx = ss(ic% + 1)
FOR iy = -1 TO -90 STEP -1
aa = iy
GOSUB 1900
NEXT iy

ic% = ic% + 1: GOTO 744

745 ix = -180
FOR iy = -90 TO 90
aa = iy: ao = ix
GOSUB 1080
NEXT iy
iy = 90
FOR ix = -179 TO 180
aa = iy: ao = ix
GOSUB 1080
NEXT ix
ix = 180
FOR iy = 89 TO -90 STEP -1
aa = iy: ao = ix
GOSUB 1080
NEXT iy
iy = -90
FOR ix = 179 TO -180 STEP -1
aa = iy: ao = ix
GOSUB 1080
NEXT ix
GOTO 770

750 xb = 180: yb = 180 / ar
x2 = xb: y2 = yb
GOSUB 4000
x2 = xb: y2 = -yb
GOSUB 4000
x2 = -xb: y2 = -yb
GOSUB 4000
x2 = -xb: y2 = yb
GOSUB 4000
x2 = xb: y2 = yb
GOSUB 4000
pv% = 0
x2 = 0: y2 = yb
GOSUB 4000
x2 = 0: y2 = -yb
GOSUB 4000
GOTO 770

755 xb = 180: yb = 180 / ar
x2 = 0: y2 = yb
GOSUB 4000
x2 = xb: y2 = 0
GOSUB 4000
x2 = 0: y2 = -yb
GOSUB 4000
x2 = -xb: y2 = 0
GOSUB 4000
x2 = 0: y2 = yb
GOSUB 4000
GOTO 770

760 x2 = 0: y2 = 90
GOSUB 4000
FOR ix = -180 TO 180
x2 = 180 * SIN(dr# * ix * k1)
y2 = 90 - 180 * COS(dr# * ix * k1)
GOSUB 4000
NEXT ix
x2 = 0: y2 = 90
GOSUB 4000
GOTO 770

765 FOR ix = -90 TO 90
x2 = 45 * SIN(dr# * ix)
y2 = 135 - 45 * COS(dr# * ix)
GOSUB 4000
NEXT ix
FOR iy = 89 TO -90 STEP -1
x2 = 45 + 135 * COS(dr# * iy)
y2 = 135 * SIN(dr# * iy)
GOSUB 4000
NEXT iy
FOR ix = 89 TO -90 STEP -1
x2 = 45 * SIN(dr# * ix)
y2 = 45 * COS(dr# * ix) - 135
GOSUB 4000
NEXT ix
FOR iy = -89 TO 90
x2 = (-45) - 135 * COS(dr# * iy)
y2 = 135 * SIN(dr# * iy)
GOSUB 4000
NEXT iy
GOTO 770

REM draw latitude and longitude lines

770 IF (ac% AND 16) = 0 THEN 780

FOR ix = 0 TO -180 STEP -gx
pv% = 0
FOR iy = -85 TO 85
la = iy: lo = ix
GOSUB 1000
NEXT iy
NEXT ix

FOR ix = gx TO 179 STEP gx
pv% = 0
FOR iy = -85 TO 85
la = iy: lo = ix
GOSUB 1000
NEXT iy
NEXT ix

780 IF (ac% AND 8) = 0 THEN 790

FOR iy = 0 TO -89 STEP -gy
pv% = 0
FOR ix = -180 TO 180
la = iy: lo = ix
GOSUB 1000
NEXT ix
NEXT iy

FOR iy = gy TO 89 STEP gy
pv% = 0
FOR ix = -180 TO 180
la = iy: lo = ix
GOSUB 1000
NEXT ix
NEXT iy

790 IF (ac% AND 64) = 0 THEN 800

FOR iy = 0 TO -gt STEP -gn
pv% = 0
ht = 2 * rd# * ATN(EXP(iy)) - 90
FOR ix = -180 TO 180
la = ht: lo = ix
GOSUB 1000
NEXT ix
NEXT iy

FOR iy = gn TO gt STEP gn
pv% = 0
ht = 2 * rd# * ATN(EXP(iy)) - 90
FOR ix = -180 TO 180
la = ht: lo = ix
GOSUB 1000
NEXT ix
NEXT iy

REM PRINT "Beginning actual draw."

800 IF (ac% AND 1) = 0 THEN 820
pv% = 0
pp% = 0

REM read map file and draw

810 IF EOF(1) THEN 820
LINE INPUT #1, tx$
REM PRINT tx$
REM P      0    1240 1 W 79793839.900 -17.53378 190.32600 -34.83044  77.71625
REM 123456789 123456789 123456789
np% = VAL(MID$(tx$, 10, 7))
PRINT np%;
dw% = 1
lv% = VAL(MID$(tx$, 17, 2))
IF ((lv% AND 1) = 0) AND (np% < 8) THEN dw% = 0
pv% = 0

FOR ix% = 1 TO np%
INPUT #1, lo, la
IF lo > 180 THEN lo = lo - 360
IF lo < -180 THEN lo = lo + 360

IF dw% = 1 THEN GOSUB 1000

NEXT ix%

GOTO 810

REM Draw new format maps
820 IF (ac% AND 4) = 0 THEN 850
821 PRINT
LINE INPUT "New format layer file name: "; fz$
IF LEN(fz$) < 1 THEN 850
OPEN fz$ FOR INPUT AS #3

pv% = 0
pp% = 0

REM read map file and draw

825 IF EOF(3) THEN 821
LINE INPUT #3, tx$
IF LEFT$(tx$, 1) = "#" THEN 825
IF LEFT$(tx$, 1) = ">" THEN pv% = 0: PRINT "."; : GOTO 825
REM IF LEFT$(tx$, 1) = ">" THEN pv% = 0: GOTO 825

REM PRINT tx$

n1 = INSTR(tx$, " ")

IF n1 < 2 THEN n1 = INSTR(tx$, CHR$(9))
IF n1 < 2 THEN 825

REM PRINT n1

lo = VAL(MID$(tx$, 1, n1 - 1))
la = VAL(LTRIM$(MID$(tx$, n1 + 1)))

REM PRINT lo, la

IF lo > 180 THEN lo = lo - 360
IF lo < -180 THEN lo = lo + 360

GOSUB 1000

GOTO 825

REM Draw a map of isolated points (stars, cities)
850 IF (ac% AND 32) = 0 THEN 900

REM PRINT "Plotting point layer"

pv% = 0
pp% = 1
860 IF EOF(3) THEN 900
INPUT #3, la, lo, mg, st$

REM PRINT la, lo, mg, st$

REM lo and mg interpreted astronomically;
REM different code used for cities
lo = -lo
mg = INT(mg + .5)

IF mg > 5 THEN 860
IF mg < 1 THEN mg = 1

REM PRINT la, lo, mg

GOSUB 1000

GOTO 860

900 PRINT

IF (ac% AND 1) > 0 THEN CLOSE 1
IF (ac% AND 4) > 0 THEN CLOSE 3

os$ = "  "
la$ = ""
wi% = -1
PRINT #2, "#define image_width "; w%
PRINT #2, "#define image_height "; h%
PRINT #2, "static unsigned char image_bits[] = {"
FOR j% = 1 TO h%
FOR i% = 1 TO wb%
t& = im%(i%, j%)
a$ = "0x" + RIGHT$("0" + HEX$((65535 AND t&) \ 256), 2)
GOSUB 10000
a$ = "0x" + RIGHT$("0" + HEX$(255 AND t&), 2)
GOSUB 10000
NEXT i%
NEXT j%
GOSUB 10100
CLOSE 2

GOTO 3

990 END

REM map drawing routine

REM coordinate transform

1000 IF ti% = 0 THEN 1050
ca = la
co = lo - a1
IF co < -180 THEN co = co + 360
IF co > 180 THEN co = co - 360

REM Fast flip, needed since conic SP only N
IF a2 = 180 THEN aa = -ca: ao = 180 - co: GOTO 1070

REM Convert to Cartesian coordinates
cy = SIN(dr# * ca)
vv = COS(dr# * ca)
cx = COS(dr# * co) * vv
cz = SIN(dr# * co) * vv

REM This matrix multiplication performs the actual rotation
qy = cy * tc - cx * ts
qx = cx * tc + cy * ts
qz = cz

IF ti% = 1 THEN 1005

REM Secondary rotation (for ease of matching)

cy = qy * ec - qz * es
cz = qz * ec + qy * es

qy = cy: qz = cz

1005 rr = SQR(qx * qx + qz * qz)
IF rr > .005 THEN aa = rd# * ATN(qy / rr): GOTO 1010
aa = 90: IF qy < 0 THEN aa = -90
ao = 0
REM Near the poles, computing longitude is meaningless
GOTO 1040

1010 IF ABS(qz) > ABS(qx) THEN 1020
ao = rd# * ATN(qz / qx)
IF qx < 0 THEN ao = ao + 180
GOTO 1040

1020 ao = 90 - rd# * ATN(qx / qz)
IF qz < 0 THEN ao = ao - 180

1040 ao = ao + a3
GOTO 1070

1050 aa = la
REM ao = (lo - a1) + a3
REM should be the equation, but is not needed, due to
REM precomputation in this case
ao = lo + a3

1070 IF ao < -180 THEN ao = ao + 360
IF ao > 180 THEN ao = ao - 360

1080 IF pp% = 1 THEN pv% = 1: GOTO 1900
IF (tp% <> 8) AND (pv% = 0) THEN 1900
ON tp% GOTO 1110, 1120, 1130, 1140, 1150, 1160, 1170, 1180, 1190, 1200, 1210

REM Azimuthal equidistant case
1110 IF aa > -30 THEN 1900
uu = ABS(ao - po)
IF uu < 8 THEN 1900
IF uu > 352 THEN 1900
pv% = 0
GOTO 1900

REM Cylindrical case
1120 IF ABS(ao) < 90 THEN 1900
IF ABS(po) < 90 THEN 1900
IF ao * po > 0 THEN 1900
pv% = 0
GOTO 1900

REM Hemispheres (N/S) most azimuthals
1130 IF aa * pa > 0 THEN 1900
pv% = 0
GOTO 1900

REM Hemispheres (E/W) globular
1140 IF ao * po > 0 THEN 1900
pv% = 0
GOTO 1900

1150 REM
1160 REM
1170 REM

REM Interrupted projection
1180 cu% = 1

REM PRINT "Was:"; aa; ao

IF aa < sp% THEN 1185
1182 IF ao > ni(cu%) THEN cu% = cu% + 1: GOTO 1182
vo = ao - ns(cu%)
dx = ns(cu%)
IF pv% = 1 THEN IF pa > sp% THEN IF pc% <> cu% THEN pv% = 0
GOTO 1189

1185 IF ao > si(cu%) THEN cu% = cu% + 1: GOTO 1185
vo = ao - ss(cu%)
dx = ss(cu%)
IF pv% = 1 THEN IF pa < sp% THEN IF pc% <> cu% THEN pv% = 0
1189 pc% = cu%

REM PRINT "Became:"; aa; ao; vo; dx

REM also check for this, in case of no interruptions
IF pv% = 0 THEN 1900
IF ABS(ao) < 90 THEN 1900
IF ABS(po) < 90 THEN 1900
IF ao * po > 0 THEN 1900
pv% = 0
GOTO 1900

REM Guyou's projection case
1190 IF ABS(aa) > 45 THEN 1195
IF ABS(ao) < 90 THEN 1900
IF ABS(po) < 90 THEN 1900
IF ao * po > 0 THEN 1900
pv% = 0
GOTO 1900

1195 IF ao * po > 0 THEN 1900
pv% = 0
GOTO 1900

REM Aitoff case
1200 IF ABS(ao) < 90 THEN 1205
IF ABS(po) < 90 THEN 1205
IF ao * po > 0 THEN 1205
pv% = 0
GOTO 1900
1205 IF aa < -70 THEN 1207
IF aa > 70 THEN 1207
GOTO 1900
1207 uu = ABS(ao - po)
IF uu < th% THEN 1900
IF uu > 360 - th% THEN 1900
pv% = 0
GOTO 1900

REM Four-square case
1210 IF aa >= 0 THEN 1900
IF pa >= 0 THEN 1900
cu% = INT(ao / 90)
IF pv% = 1 THEN IF pc% <> cu% THEN pv% = 0
1219 pc% = cu%
GOTO 1900

1900 po = ao: pa = aa

REM PRINT "Plotting:"; aa; ao; vo; dx

ON px% GOTO 1910, 1920, 1930, 1940, 1950, 1960

1910 ON py% GOTO 2010, 2020, 2030, 2040, 2050, 2060, 2070, 2080, 2090, 2100
1920 ON py% GOTO 2110, 2120, 2130, 2140, 2150, 2160, 2170, 2180, 2190, 2200
1930 ON py% GOTO 2210, 2220, 2230, 2240, 2250, 2260, 2270, 2280, 2290, 2300
1940 ON py% GOTO 2310, 2320, 2330, 2340, 2350, 2360, 2370, 2380, 2390, 2400
1950 ON py% GOTO 2410, 2420, 2430, 2440, 2450, 2460, 2470, 2480, 2490, 2500
1960 ON py% GOTO 2510, 2520, 2530, 2540, 2550, 2560

REM Mercator
2010 IF aa > ql THEN pv% = 0: RETURN
IF aa < -qs THEN pv% = 0: RETURN
y2 = qv + rd# * LOG((1 / COS(dr# * aa)) + TAN(dr# * aa))
x2 = ao
GOTO 2900

REM Gall's
2020 IF ABS(aa) > 89.9 THEN pv% = 0: RETURN
y2 = TAN(dr# * .5 * aa) * k1
x2 = ao
GOTO 2900

REM Miller's
2030 IF ABS(aa) > 89.9 THEN pv% = 0: RETURN
y2 = rd# * 1.25 * LOG((1 / COS(.8 * dr# * aa)) + TAN(.8 * dr# * aa))
x2 = ao
GOTO 2900

REM Plate Carree
2040 IF ABS(aa) > 89.9 THEN pv% = 0: RETURN
x2 = ao: y2 = aa
GOTO 2900

REM Cylindrical Equal-Area
2050 IF ABS(aa) > 89.9 THEN pv% = 0: RETURN
y2 = rd# * SIN(dr# * aa)
x2 = ao
GOTO 2900

REM Gnomonic
2060 IF ABS(aa) < gb THEN pv% = 0: RETURN
IF aa < 0 THEN 2065
rr = gs * TAN(dr# * (90 - aa))
x2 = rr * SIN(dr# * ao)
y2 = -rr * COS(dr# * ao)
IF (ABS(x2) > 90) OR (ABS(y2) > 90) THEN pv% = 0: RETURN
x2 = x2 - 90
GOTO 2900
2065 rr = gs * TAN(dr# * (90 + aa))
x2 = rr * SIN(dr# * ao)
y2 = -rr * COS(dr# * ao)
IF (ABS(x2) > 90) OR (ABS(y2) > 90) THEN pv% = 0: RETURN
x2 = 90 - x2
GOTO 2900

REM Stereographic
2070 IF aa < 0 THEN 2075
rr = 90 * TAN(dr# * (90 - aa) / 2)
x2 = rr * SIN(dr# * ao) - 90
y2 = -rr * COS(dr# * ao)
GOTO 2900
2075 rr = 90 * TAN(dr# * (90 + aa) / 2)
x2 = 90 - rr * SIN(dr# * ao)
y2 = -rr * COS(dr# * ao)
GOTO 2900

REM Azimuthal Equidistant
2080 IF aa < 0 THEN 2085
rr = 90 - aa
x2 = rr * SIN(dr# * ao) - 90
y2 = -rr * COS(dr# * ao)
GOTO 2900
2085 rr = 90 + aa
x2 = 90 - rr * SIN(dr# * ao)
y2 = -rr * COS(dr# * ao)
GOTO 2900

REM Lambert's Azimuthal Equal-Area
2090 IF aa < 0 THEN 2095
rr = 90 * SQR(1 - SIN(dr# * aa))
x2 = rr * SIN(dr# * ao) - 90
y2 = -rr * COS(dr# * ao)
GOTO 2900
2095 rr = 90 * SQR(1 + SIN(dr# * aa))
x2 = 90 - rr * SIN(dr# * ao)
y2 = -rr * COS(dr# * ao)
GOTO 2900

REM Orthographic
2100 IF aa < 0 THEN 2105
rr = 90 * SIN(dr# * (90 - aa))
x2 = rr * SIN(dr# * ao) - 90
y2 = -rr * COS(dr# * ao)
GOTO 2900
2105 rr = 90 * SIN(dr# * (90 + aa))
x2 = 90 - rr * SIN(dr# * ao)
y2 = -rr * COS(dr# * ao)
GOTO 2900

REM Azimuthal Equidistant (single circle)
2110 rr = 90 - aa
x2 = rr * SIN(dr# * ao)
y2 = -rr * COS(dr# * ao)
GOTO 2900

REM Lambert's Conic Conformal
2120 IF aa < -qs THEN pv% = 0: RETURN
r = k2 * TAN(.5 * dr# * (90 - aa)) ^ k1
a = dr# * k1 * ao
x2 = r * SIN(a)
y2 = 90 - r * COS(a)
GOTO 2900

REM Simple conic
2130 r = k2 + (sp - aa)
a = dr# * k1 * ao
x2 = r * SIN(a)
y2 = k2 - r * COS(a)
GOTO 2900

REM Albers' Conical Equal-Area
2140 r = k4 * SQR(k3 - SIN(dr# * aa))
a = dr# * k1 * ao
x2 = r * SIN(a)
y2 = k5 - r * COS(a)
GOTO 2900

REM Sinusoidal
2150 x2 = vo * COS(dr# * aa) + dx
y2 = aa
GOTO 2900

REM Bonne's projection
2160 r = k1 + (sp - aa)
a = k2 * dx + vo * COS(dr# * aa) / r
x2 = r * SIN(a)
y2 = (k1 + sp) - r * COS(a)
GOTO 2900

REM Mollweide
2170 pp = dr# * aa
k = pi# * .5 * SIN(pp)
IF aa = 0 THEN ps = 0: GOTO 2175
IF ABS(aa) = 90 THEN ps = pp: GOTO 2175
2172 sp = SIN(pp): cp = COS(pp)
ps = pp + (k - (pp + sp * cp)) / (1 + cp * cp - sp * sp)
IF ABS(pp - ps) > .0001 THEN pp = ps: GOTO 2172
2175 x2 = vo * COS(ps) + dx
y2 = 90 * SIN(ps)
GOTO 2900

REM Eckert IV Projection
2180 pp = dr# * aa
k = (.5 * pi# + 2) * SIN(dr# * aa)
IF aa = 0 THEN ps = 0: GOTO 2185
IF ABS(aa) = 90 THEN ps = pp: GOTO 2185
2182 sp = SIN(pp): cp = COS(pp)
ps = pp + (k - (pp + sp + sp + cp * sp)) / (1 + cp + cp + cp * cp - sp * sp)
IF ABS(pp - ps) > .0001 THEN pp = ps: GOTO 2182
2185 x2 = vo * .5 * (1 + COS(ps)) + dx
y2 = 90 * SIN(ps)
GOTO 2900

REM MacBryde's Flat Polar Quartic
2190 k = 1.707107 * SIN(dr# * aa)
IF aa = 0 THEN ps = 0: GOTO 2195
pp = dr# * aa: IF ABS(aa) = 90 THEN ps = pp: GOTO 2195
2192 ps = pp + (k - (SIN(pp) + SIN(.5 * pp))) / (COS(pp) + COS(pp / 2) / 2)
IF ABS(pp - ps) > .0001 THEN pp = ps: GOTO 2192
2195 x2 = vo * (1 + 2 * COS(ps) / COS(.5 * ps)) / 3 + dx
y2 = 2 * rd# * SIN(.5 * ps)
GOTO 2900

REM Adams' Equal-Area
2200 x2 = vo * (COS(dr# * aa) / COS(dr# * .5 * aa)) + dx
y2 = 2 * rd# * SIN(dr# * .5 * aa)
GOTO 2900

REM Parabolic Stereographic
2210 yp = TAN(dr# * .5 * aa)
y2 = 2 * rd# * yp
x2 = vo * (1 - yp * yp) + dx
GOTO 2900

REM Loximuthal
2220 y2 = aa
IF aa = sp THEN s = 1: GOTO 2225
IF ABS(aa) = 90 THEN s = 0: GOTO 2225
y = rd# * LOG((1 / COS(dr# * aa)) + TAN(dr# * aa))
s = (y2 - sp) / (y - k1)
2225 x2 = ao * s
GOTO 2900

REM Polyconic
2230 sa = SGN(aa): ma = ABS(aa)
y = dr# * ma
IF ma < 3 THEN 2235

r = rd# / TAN(y)
co = vo * SIN(y)
xs = dr# * co
xx = r * SIN(xs)
yy = ma + r * (1 - COS(xs))
GOTO 2237

REM Parabolic case
2235 x = dr# * vo
xx = vo * (1 - .6666667 * y * y) * (1 - .1666667 * x * x * y * y)
yy = ma * (1 + .5 * x * x) * (1 - .5 * y * y)

2237 x2 = xx + dx
y2 = sa * yy
GOTO 2900

REM Globular
2240 io = ao + 90: dx = -90: IF io > 90 THEN io = io - 180: dx = 90
qa = ABS(aa): sy = SGN(aa)
qo = ABS(io): sx = SGN(io)

REM code for edges of area
IF qa = 90 THEN x = 0: y = 90: GOTO 2247
IF qa = 0 THEN x = qo: y = 0: GOTO 2247
IF qo = 0 THEN x = 0: y = qa: GOTO 2247
IF qo = 90 THEN x = 90 * COS(dr# * qa): y = 90 * SIN(dr# * qa): GOTO 2247

REM parabolic cases
IF qo < 1 THEN y = qa: x = qo - qo * (qa / 90) ^ 2: GOTO 2247
IF qa < 1 THEN x = qo: y = qa + (90 * SIN(dr# * qa) - qa) * (qo / 90) ^ 2: GOTO 2247

REM calculate radii and centers of intersecting circles
a = qo: b = 90
c2 = (a * a + b * b) / (2 * a) - a
r2 = c2 + a
a = 90 * SIN(dr# * qa) - qa: b = 90 * COS(dr# * qa)
c1 = (a * a + b * b) / (2 * a) - a
r1 = c1 + a
c1 = r1 + qa

REM calculate intersection of circles
d2 = c1 * c1 + c2 * c2
D = SQR(d2)
rr = r1 * r1
zx = (d2 + rr - r2 * r2) / (2 * D)
zy = SQR(rr - zx * zx)

REM coordinate conversion
y = c1 - zx * c1 / D - zy * c2 / D
x = zy * c1 / D - zx * c2 / D

REM replace sign
2247 x2 = x * sx + dx
y2 = y * sy
GOTO 2900

REM Van der Grinten
2250  qa = ABS(aa): sy = SGN(aa)
qo = ABS(ao): sx = SGN(ao)

eh = 360 * qa / (360 - (qa + qa))
ch = 360 * qa / (180 + SQR(32400 - 4 * qa * qa))

REM code for edges of area
IF qa = 90 THEN x = 0: y = 90: GOTO 2257
IF qa = 0 THEN x = qo: y = 0: GOTO 2257
IF qo = 0 THEN x = 0: y = ch: GOTO 2257
IF qo = 180 THEN x = SQR(32400 - eh * eh): y = eh: GOTO 2257

REM parabolic cases
IF qo < 1 THEN y = ch: x = qo - qo * (ch / 180) ^ 2: GOTO 2257
IF qa < 9 THEN x = qo: y = qa + (eh - ch) * (qo / 180) ^ 2: GOTO 2257

REM calculate radii and centers of intersecting circles
a = qo: b = 180
c2 = (a * a + b * b) / (2 * a) - a
r2 = c2 + a
a = eh - ch: b = SQR(32400 - eh * eh)
c1 = (a * a + b * b) / (2 * a) - a
r1 = c1 + a
c1 = r1 + ch

REM calculate intersection of circles
d2 = c1 * c1 + c2 * c2
D = SQR(d2)
rr = r1 * r1
zx = (d2 + rr - r2 * r2) / (2 * D)
zy = SQR(rr - zx * zx)

REM coordinate conversion
y = c1 - zx * c1 / D - zy * c2 / D
x = zy * c1 / D - zx * c2 / D

REM replace sign
2257 x2 = x * sx
y2 = y * sy
GOTO 2900

REM Van der Grinten IV
2260 qa = ABS(aa): sy = SGN(aa)
qo = ABS(ao): sx = SGN(ao)

REM code for edges of area
IF qa = 90 THEN x = 0: y = 90: GOTO 2267
IF qa = 0 THEN x = qo: y = 0: GOTO 2267
IF qo = 0 THEN x = 0: y = qa: GOTO 2267

REM parabolic case 1
IF qo < 3 THEN y = qa: x = qo - qo * (qa / 90) ^ 2: GOTO 2267

REM calculate radii and centers of intersecting circles
a = qo: b = 90
c2 = (a * a + b * b) / (2 * a) - a
r2 = c2 + a
se = qa / (270 - qa)
s2 = se * se
xi = (SQR(129600 * s2 * s2 - 4 * (s2 + 1) * (32400 * s2 - 8100)) - 360 * s2) / (2 * (s2 + 1))
yi = SQR(8100 - xi * xi)

REM Parabolic case 2
IF qa < 8 THEN x = qo: y = qa + (yi - qa) * (qo / 90) ^ 2: GOTO 2267

a = yi - qa
b = xi
c1 = (a * a + b * b) / (2 * a) - a
r1 = c1 + a
c1 = r1 + qa

REM calculate intersection of circles
d2 = c1 * c1 + c2 * c2
D = SQR(d2)
rr = r1 * r1
zx = (d2 + rr - r2 * r2) / (2 * D)
zy = SQR(rr - zx * zx)

REM coordinate conversion
y = c1 - zx * c1 / D - zy * c2 / D
x = zy * c1 / D - zx * c2 / D

REM replace sign
2267 x2 = x * sx
y2 = y * sy
GOTO 2900

REM Savard Egg
2270 qa = ABS(aa): qo = ABS(ao)
tk = 8100 - qa * qa
x = SQR(tk) * qo / 90
yz = qa + qa * qa * qa / 24300
dl = qa / .75 - yz
IF qa < 1 THEN 2275
r = dl / 2 + 2 * tk / dl
y = yz + r - SQR(r * r - x * x)
GOTO 2277

REM Parabolic case
2275 y = yz + dl * qo * qo / 32400

2277 x2 = x * SGN(ao)
y2 = y * SGN(aa)
GOTO 2900

REM Modified Van der Grinten IV
2280 qa = 90 * TAN(dr# * .5 * ABS(aa)): sy = SGN(aa)
qo = .25 * pi# * ABS(ao): sx = SGN(ao)

REM code for edges of area
IF qa = 90 THEN x = 0: y = 90: GOTO 2287
IF qa = 0 THEN x = qo: y = 0: GOTO 2287
IF qo = 0 THEN x = 0: y = qa: GOTO 2287

REM parabolic case 1
IF qo < 3 THEN y = qa: x = qo - qo * (qa / 90) ^ 2: GOTO 2287

REM calculate radii and centers of intersecting circles
a = qo: b = 90
c2 = (a * a + b * b) / (2 * a) - a
r2 = c2 + a
se = qa / (270 - qa)
s2 = se * se
xi = (SQR(129600 * s2 * s2 - 4 * (s2 + 1) * (32400 * s2 - 8100)) - 360 * s2) / (2 * (s2 + 1))
yi = SQR(8100 - xi * xi)

REM Parabolic case 2
IF qa < 8 THEN x = qo: y = qa + (yi - qa) * (qo / 90) ^ 2: GOTO 2287

a = yi - qa
b = xi
c1 = (a * a + b * b) / (2 * a) - a
r1 = c1 + a
c1 = r1 + qa

REM calculate intersection of circles
d2 = c1 * c1 + c2 * c2
D = SQR(d2)
rr = r1 * r1
zx = (d2 + rr - r2 * r2) / (2 * D)
zy = SQR(rr - zx * zx)

REM coordinate conversion
y = c1 - zx * c1 / D - zy * c2 / D
x = zy * c1 / D - zx * c2 / D

REM replace sign
2287 x2 = 1.27324 * x * sx
y2 = 1.27324 * y * sy
GOTO 2900

REM Globular (whole world)
2290 qa = ABS(aa): sy = SGN(aa)
qo = ABS(ao): sx = SGN(ao)

REM code for edges of area
IF qa = 90 THEN x = 0: y = 90: GOTO 2297
IF qa = 0 THEN x = qo: y = 0: GOTO 2297
IF qo = 0 THEN x = 0: y = qa: GOTO 2297
IF qo = 90 THEN x = 90 * COS(dr# * qa): y = 90 * SIN(dr# * qa): GOTO 2297

REM parabolic cases
IF qo < 1 THEN y = qa: x = qo - qo * (qa / 90) ^ 2: GOTO 2297
IF qa < 1 THEN x = qo: y = qa + (90 * SIN(dr# * qa) - qa) * (qo / 90) ^ 2: GOTO 2297

REM calculate radii and centers of intersecting circles
a = qo: b = 90
c2 = (a * a + b * b) / (2 * a) - a
r2 = c2 + a
a = 90 * SIN(dr# * qa) - qa: b = 90 * COS(dr# * qa)
c1 = (a * a + b * b) / (2 * a) - a
r1 = c1 + a
c1 = r1 + qa

REM calculate intersection of circles
d2 = c1 * c1 + c2 * c2
D = SQR(d2)
rr = r1 * r1
zx = (d2 + rr - r2 * r2) / (2 * D)
zy = SQR(rr - zx * zx)

REM coordinate conversion
y = c1 - zx * c1 / D - zy * c2 / D
x = zy * c1 / D - zx * c2 / D

REM replace sign
2297 x2 = x * sx
y2 = y * sy
GOTO 2900

REM Lagrange Conformal
2300 oo = .5 * ao
w = SQR(COS(dr# * aa)) / COS(.5 * dr# * aa)
GOSUB 3000
ap = (90 - w) * SGN(aa)
ca = COS(dr# * ap)
dn = 1 + COS(dr# * oo) * ca
x2 = 180 * ca * SIN(dr# * oo) / dn
y2 = 180 * SIN(dr# * ap) / dn
GOTO 2900

REM August's Conformal Projection of the World on a Two-Cusped Epicycloid
2310 oo = .5 * ao
w = SQR(COS(dr# * aa)) / COS(.5 * dr# * aa)
GOSUB 3000
ap = (90 - w) * SGN(aa)
ca = COS(dr# * ap)
dn = 1 + COS(dr# * oo) * ca
xx = SIN(dr# * ap) / dn
yy = ca * SIN(dr# * oo) / dn
x2 = (yy * yy * yy + 3 * yy - 3 * xx * xx * yy) * 45
y2 = (3 * xx * yy * yy + 3 * xx - xx * xx * xx) * 45
GOTO 2900

REM Cosine-based
2320 oo = .5 * ao
w = SQR(COS(dr# * aa)) / COS(.5 * dr# * aa)
GOSUB 3000
ap = (90 - w) * SGN(aa)
ca = COS(dr# * ap)
dn = 1 + COS(dr# * oo) * ca
xx = .5 * pi# * ca * SIN(dr# * oo) / dn
yy = .5 * pi# * SIN(dr# * ap) / dn
r1 = EXP(xx): r2 = 1 / r1
x2 = COS(yy) * (r1 - r2) * 39.1
y2 = SIN(yy) * (r1 + r2) * 39.1
GOTO 2900

REM Eisenlohr
2330 oo = .5 * ao: oa = .5 * aa
s1 = SIN(dr# * oo): c1 = COS(dr# * oo)
c2 = COS(dr# * aa): c3 = COS(dr# * oa)
c4 = SQR(c2 / 2)
t = SIN(dr# * oa) / (c3 + 2 * c4 * c1)
c = SQR(2 / (1 + (t * t)))
v = SQR((c3 + c4 * (c1 + s1)) / (c3 + c4 * (c1 - s1)))
x2 = 166.9721 * (c * (v - 1 / v) - 2 * LOG(v))
y2 = 166.9721 * (c * t * (v + 1 / v) - 2 * ATN(t))
GOTO 2900

REM Guyou's doubly-periodic projection
2340 dp = -90: IF ao > 0 THEN dp = 90
oo = ao - dp
IF (ABS(aa) > 80) THEN 2345
im = SGN(aa) * LOG(TAN(dr# * .5 * (ABS(aa) + 90)))
re = dr# * oo

mm = .5: GOSUB 3200

x2 = 48.54173 * x + dp
y2 = 48.54173 * y
GOTO 2900

REM im goes to infinity as you approach the poles, so this is needed
2345 v1 = 0: v2 = 90: v3 = 0
wa = aa: wo = ao
GOSUB 3300

im = SGN(ca) * LOG(TAN(dr# * .5 * (ABS(ca) + 90)))
re = dr# * (co - dp)

mm = .5: GOSUB 3200

IF dp < 0 THEN 2347

x2 = 48.54173 * y + 90
y2 = -48.54173 * x
GOTO 2900

2347 x2 = -48.54173 * y - 90
y2 = 48.54173 * x
GOTO 2900

REM World on a rectangle
2350 IF (ABS(aa) > 85) THEN 2355
im = SGN(aa) * .5 * LOG(TAN(dr# * .5 * (ABS(aa) + 90)))
re = dr# * .5 * ao

mm = k1: GOSUB 3200

x2 = k2 * x
y2 = k2 * y
GOTO 2900

REM im goes to infinity as you approach the poles, so this is needed
2355 w = SQR(COS(dr# * aa)) / COS(.5 * dr# * aa)
GOSUB 3000
qa = (90 - w) * SGN(aa)
qo = ao / 2

v1 = 0: v2 = 90: v3 = 0
wa = qa: wo = qo
GOSUB 3300

im = SGN(ca) * LOG(TAN(dr# * .5 * (ABS(ca) + 90)))
re = dr# * co

mm = 1 - k1: GOSUB 3200

x2 = k2 * y
y2 = -k2 * x
GOTO 2900

REM World on a diamond
2360 w = SQR(COS(dr# * aa)) / COS(.5 * dr# * aa)
GOSUB 3000
qa = (90 - w) * SGN(aa)
REM qa = COS(dr# * ap)
qo = ao / 2

v1 = 90: v2 = 45: v3 = 90
wa = qa: wo = qo
GOSUB 3300

IF (ABS(ca) > 85) THEN 2365
im = SGN(ca) * LOG(TAN(dr# * .5 * (ABS(ca) + 90)))
re = dr# * co

mm = .5: GOSUB 3200

x2 = 48.54173 * (x - y)
y2 = 48.54173 * (x + y)
GOTO 2900

REM im goes to infinity as you approach the poles, so this is needed
2365 v1 = 90: v2 = 90: v3 = 90
wa = ca: wo = co
GOSUB 3300

im = SGN(ca) * LOG(TAN(dr# * .5 * (ABS(ca) + 90)))
re = dr# * co

mm = .5: GOSUB 3200

x2 = -48.54173 * (x + y)
y2 = -48.54173 * (y - x)
GOTO 2900

REM Hammer-Aitoff Projection
2370 di = SQR(1 + COS(dr# * aa) * COS(dr# * .5 * ao))
y2 = 90 * SIN(dr# * aa) / di
x2 = 180 * COS(dr# * aa) * SIN(dr# * .5 * ao) / di
GOTO 2900

REM Generalized Aitoff Equal-Area
2380 w = sh * SIN(dr# * aa): GOSUB 3000: ap = w
di = SQR(1 + COS(dr# * ap) * COS(dr# * .5 * ao))
y2 = 90 * ex * SIN(dr# * ap) / di
x2 = 180 * COS(dr# * ap) * SIN(dr# * .5 * ao) / di
GOTO 2900

REM Aitoff-Wagner
2390 ap = .7777778# * aa: op = .2777778# * ao
wx = SIN(dr# * op) * COS(dr# * ap)
wy = SIN(dr# * ap)
GOSUB 3020
az = w
wa = ap
wo = op
v1 = 0: v2 = -90: v3 = 0
GOSUB 3300
ra = 90 - ca
x2 = 3.6 * ra * COS(dr# * az)
y2 = ra * SIN(dr# * az) * 1.285714
GOTO 2900

REM Aitoff
2400 REM Calculate x and y for an orthographic projection
wx = SIN(.5 * dr# * ao) * COS(dr# * aa)
wy = SIN(dr# * aa)
GOSUB 3020
az = w
wa = aa
wo = .5 * ao
v1 = 0: v2 = -90: v3 = 0
GOSUB 3300
ra = 90 - ca
x2 = 2 * ra * COS(dr# * az)
y2 = ra * SIN(dr# * az)
GOTO 2900

REM Winkel's Tripel
2410 REM Calculate x and y for an orthographic projection
wx = SIN(.5 * dr# * ao) * COS(dr# * aa)
wy = SIN(dr# * aa)
GOSUB 3020
az = w
wa = aa
wo = .5 * ao
v1 = 0: v2 = -90: v3 = 0
GOSUB 3300
ra = 90 - ca
x2 = 1.132474 * (ra * COS(dr# * az) + .3830222 * ao)
y2 = 1.132474 * (.5 * ra * SIN(dr# * az) + aa / 2)
GOTO 2900

REM Not Quite Generalized Aitoff-Wagner
2420 ap = sh * aa: op = sh * ao
wx = SIN(.5 * dr# * op) * COS(dr# * ap)
wy = SIN(dr# * ap)
GOSUB 3020
az = w
wa = ap
wo = .5 * op
v1 = 0: v2 = -90: v3 = 0
GOSUB 3300
ra = 90 - ca
x2 = 2 * ra * COS(dr# * az) * ex
y2 = ra * SIN(dr# * az) * ex
GOTO 2900

REM Rectangular Polyconic
2430 sa = SGN(aa): ma = ABS(aa): so = SGN(vo): mo = ABS(vo)
oo = .5 * mo
IF sp = 0 THEN t = 2 * ATN(dr# * oo * SIN(dr# * ma)): co = mo: GOTO 2432
co = 2 * rd# * TAN(dr# * .5 * mo * k1) / k1
2432 IF ABS(ma) < 3 THEN 2435
t = 2 * ATN(.5 * dr# * co * SIN(dr# * ma))
r = 1 / TAN(dr# * ma)
xx = rd# * r * SIN(t) + dx
yy = ma + rd# * r * (1 - COS(t))
GOTO 2437

REM Parabolic case
2435 x = dr# * co: y = dr# * ma
xx = co * (1 - .6666667 * y * y) * (1 - .1666667 * x * x * y * y)
yy = ma * (1 + .5 * x * x) * (1 - .5 * y * y)

2437 x2 = xx * so + dx
y2 = yy * sa
GOTO 2900

REM Four Corners
2440 nn = TAN(dr# * (45 + .5 * aa))
qa = 2 * rd# * ATN(nn * nn) - 90
qo = ao + 180: dq = INT(qo / 90)
IF dq = 4 THEN dq = 3
qo = 2 * ((qo - 90 * dq) - 45)

v1 = 90: v2 = 45: v3 = 90
wa = qa: wo = qo
GOSUB 3300

IF (ABS(ca) > 85) THEN 2441
im = SGN(ca) * LOG(TAN(dr# * .5 * (ABS(ca) + 90)))
re = dr# * co

mm = .5: GOSUB 3200

x0 = x
y0 = y
GOTO 2442

REM im goes to infinity as you approach the poles, so this is needed
2441 v1 = 90: v2 = 90: v3 = 90
wa = ca: wo = co
GOSUB 3300

im = SGN(ca) * LOG(TAN(dr# * .5 * (ABS(ca) + 90)))
re = dr# * co

mm = .5: GOSUB 3200

x0 = -y
y0 = x

2442 ON dq + 1 GOTO 2448, 2445, 2446, 2447

2445 x2 = 48.54173 * x0 - 90
y2 = 48.54173 * y0 - 90
GOTO 2900

2446 x2 = 90 - 48.54173 * y0
y2 = 48.54173 * x0 - 90
GOTO 2900

2447 x2 = 90 - 48.54173 * x0
y2 = 90 - 48.54173 * y0
GOTO 2900

2448 x2 = 48.54173 * y0 - 90
y2 = 90 - 48.54173 * x0
GOTO 2900

REM Interrupted conformal
2450 nn = TAN(dr# * (45 + .5 * aa))
qa = 2 * rd# * ATN(nn * nn) - 90
qo = ao + 180: dq = INT(qo / 90)
IF dq = 4 THEN dq = 3
qo = 2 * ((qo - 90 * dq) - 45)

v1 = 90: v2 = 45: v3 = 90
wa = qa: wo = qo
GOSUB 3300

IF (ABS(ca) > 85) THEN 2451
im = SGN(ca) * LOG(TAN(dr# * .5 * (ABS(ca) + 90)))
re = dr# * co

mm = .5: GOSUB 3200

x0 = x
y0 = y
GOTO 2452

REM im goes to infinity as you approach the poles, so this is needed
2451 v1 = 90: v2 = 90: v3 = 90
wa = ca: wo = co
GOSUB 3300

im = SGN(ca) * LOG(TAN(dr# * .5 * (ABS(ca) + 90)))
re = dr# * co

mm = .5: GOSUB 3200

x0 = -y
y0 = x

2452 ON dq + 1 GOTO 2458, 2455, 2456, 2457

2455 x = .2696763 * x0 - .5
y = .2696763 * y0 - .5
GOTO 2459

2456 x = .5 - .2696763 * y0
y = .2696763 * x0 - .5
GOTO 2459

2457 x = .5 - .2696763 * x0
y = .5 - .2696763 * y0
GOTO 2459

2458 x = .2696763 * y0 - .5
y = .5 - .2696763 * x0
GOTO 2459

2459 x2 = 20 * (5 * x - x * x * x * x * x + 10 * x * x * x * y * y - 5 * x * y * y * y * y)
y2 = 20 * (5 * y - 5 * x * x * x * x * y + 10 * x * x * y * y * y - y * y * y * y * y)
GOTO 2900

REM Savard Conventional
2460 qa = ABS(aa): sy = SGN(aa)
qo = ABS(ao): sx = SGN(ao)

REM code for edges of area
IF qa = 90 THEN x = 45 * SIN(.5 * dr# * qo): y = 135 - 45 * COS(.5 * dr# * qo): GOTO 2467
IF qa = 0 THEN x = qo: y = 0: GOTO 2467
IF qo = 0 THEN x = 0: y = qa: GOTO 2467

REM parabolic cases
IF qo < 3 THEN y = qa: x = qo - qo * (qa / 115.489) ^ 2: GOTO 2467
IF qa < 5 THEN x = qo: y = qa: GOTO 2467

REM calculate radii and centers of intersecting circles
d1 = 45 * SIN(.5 * dr# * qo)
a = qo - d1: b = 135 - 45 * COS(.5 * dr# * qo)
c2 = (a * a + b * b) / (2 * a) - a
r2 = c2 + a
c2 = r2 - qo
z = SIN(dr# * qa)
w = (z / 1.5) + z * z * z * z * z / 6 + z * z * z / 6
GOSUB 3000
a = 135 * SIN(dr# * w) - qa: b = 45 + 135 * COS(dr# * w)
c1 = (a * a + b * b) / (2 * a) - a
r1 = c1 + a
c1 = r1 + qa

REM calculate intersection of circles
d2 = c1 * c1 + c2 * c2
D = SQR(d2)
rr = r1 * r1
zx = (d2 + rr - r2 * r2) / (2 * D)
zy = SQR(rr - zx * zx)

REM coordinate conversion
y = c1 - zx * c1 / D - zy * c2 / D
x = zy * c1 / D - zx * c2 / D

REM replace sign
2467 x2 = x * sx
y2 = y * sy
GOTO 2900

REM Savard Conventional II
2470 qa = ABS(aa): sy = SGN(aa)
qo = ABS(ao): sx = SGN(ao)

REM code for edges of area
IF qa = 90 THEN x = 45 * SIN(.5 * dr# * qo): y = 135 - 45 * COS(.5 * dr# * qo): GOTO 2477
IF qa = 0 THEN x = qo: y = 0: GOTO 2477
IF qo = 0 THEN x = 0: y = qa: GOTO 2477

REM parabolic cases
IF qo < 3 THEN y = qa: x = qo - qo * (qa / 115.489) ^ 2: GOTO 2477
IF qa < 5 THEN x = qo: y = qa: GOTO 2477

REM calculate radii and centers of intersecting circles
t = qo / 180
a = qo: b = 90 + 37.27922 * t
c2 = (a * a + b * b) / (2 * a) - a
r2 = c2 + a
c2 = r2 - qo
z = SIN(dr# * qa)
w = (z / 1.5) + z * z * z * z * z / 6 + z * z * z / 6
GOSUB 3000
a = 135 * SIN(dr# * w) - qa: b = 45 + 135 * COS(dr# * w)
c1 = (a * a + b * b) / (2 * a) - a
r1 = c1 + a
c1 = r1 + qa

REM calculate intersection of circles
d2 = c1 * c1 + c2 * c2
D = SQR(d2)
rr = r1 * r1
zx = (d2 + rr - r2 * r2) / (2 * D)
zy = SQR(rr - zx * zx)

REM coordinate conversion
y = c1 - zx * c1 / D - zy * c2 / D
x = zy * c1 / D - zx * c2 / D

REM replace sign
2477 x2 = x * sx
y2 = y * sy
GOTO 2900

REM Modified Van der Grinten IV
2480 qa = ABS(aa): sy = SGN(aa)
qo = ABS(ao): sx = SGN(ao)

REM code for edges of area
IF qa = 90 THEN x = 0: y = 90: GOTO 2487
IF qa = 0 THEN x = qo: y = 0: GOTO 2487
IF qo = 0 THEN x = 0: y = qa: GOTO 2487

REM parabolic case 1
IF qo < 3 THEN y = qa: x = qo - qo * (qa / 90) ^ 2: GOTO 2487

REM calculate radii and centers of intersecting circles
a = qo: b = 90
c2 = (a * a + b * b) / (2 * a) - a
r2 = c2 + a
se = qa / (k2 - qa)
s2 = se * se
xi = (SQR(k4 * s2 * s2 - 4 * (s2 + 1) * (k3 * s2 - 8100)) - k1 * s2) / (2 * (s2 + 1))
yi = SQR(8100 - xi * xi)

REM Parabolic case 2
IF qa < 9 THEN x = qo: y = qa + (yi - qa) * (qo / 90) ^ 2: GOTO 2487

a = yi - qa
b = xi
c1 = (a * a + b * b) / (2 * a) - a
r1 = c1 + a
c1 = r1 + qa

REM calculate intersection of circles
d2 = c1 * c1 + c2 * c2
D = SQR(d2)
rr = r1 * r1
zx = (d2 + rr - r2 * r2) / (2 * D)
zy = SQR(rr - zx * zx)

REM coordinate conversion
y = c1 - zx * c1 / D - zy * c2 / D
x = zy * c1 / D - zx * c2 / D

REM replace sign
2487 x2 = x * sx
y2 = y * sy
GOTO 2900

REM Squashed Conformal
2490 oo = .5 * ao
w = SQR(COS(dr# * aa)) / COS(.5 * dr# * aa)
GOSUB 3000
ap = (90 - w) * SGN(aa)
ca = COS(dr# * ap)
dn = 1 + COS(dr# * oo) * ca
xx = SIN(dr# * ap) / dn
yy = ca * SIN(dr# * oo) / dn
x = (yy * yy * yy + 3 * yy - 3 * xx * xx * yy) * .18
y = (3 * xx * yy * yy + 3 * xx - xx * xx * xx) * .18
yy = (x * x * x * x * x - 10 * x * x * x * y * y + 5 * x * y * y * y * y - 5 * x) * .25
xx = (5 * x * x * x * x * y - 10 * x * x * y * y * y + y * y * y * y * y - 5 * y) * .25
IF k1 = 2 THEN x2 = -213 * yy: y2 = -213 * xx: GOTO 2900
x = (yy * yy * yy + 3 * yy - 3 * xx * xx * yy) * .18
y = (3 * xx * yy * yy + 3 * xx - xx * xx * xx) * .18
IF k1 = 3 THEN x2 = -317 * x: y2 = -317 * y: GOTO 2900
yy = (x * x * x * x * x - 10 * x * x * x * y * y + 5 * x * y * y * y * y - 5 * x) * .18
xx = (5 * x * x * x * x * y - 10 * x * x * y * y * y + y * y * y * y * y - 5 * y) * .18
IF k1 = 4 THEN x2 = 360 * yy: y2 = 360 * xx: GOTO 2900
y2 = -(xx * xx * xx * xx * xx - 10 * xx * xx * xx * yy * yy + 5 * xx * yy * yy * yy * yy - 5 * xx) * 73
x2 = -(5 * xx * xx * xx * xx * yy - 10 * xx * xx * yy * yy * yy + yy * yy * yy * yy * yy - 5 * yy) * 73
GOTO 2900

REM Ginzoid
2500 oa = aa / k3: oo = vo / k1
sa = SGN(oa): ma = ABS(oa)
ve = k5 * ABS(aa / k1)
rr = dr# * ABS(aa / k2)
y = dr# * ma
IF ma < 3 THEN 2505

t = TAN(rr)
r = rd# / t
co = oo * (k4 + (1 - k4) * COS(y)) * t
xs = dr# * co
xx = r * SIN(xs)
yy = ve + r * (1 - COS(xs))
GOTO 2507

REM Parabolic case
2505 xe = dr# * ABS(aa / (k1 * k2))
xx = vo * (1 - .6666667 * y * y) * (1 - .1666667 * xe * xe * y * y)
yy = ma * (1 + .5 * xe * xe) * (1 - .5 * y * y)

2507 x2 = xx * k1 + dx
y2 = sa * yy * k1
GOTO 2900

REM Larrivee
2510 x2 = ao * (.5 + .5 * SQR(COS(dr# * aa)))
y2 = aa / (COS(dr# * .5 * aa) * COS(dr# * ao / 6))
GOTO 2900

REM Laskowski Tri-Optimal
2520 a = dr# * aa: o = dr# * ao
x2 = rd# * o * (.975534 - .119161 * a * a - .0547009# * a * a * a * a - .0143059 * o * o * a * a)
y2 = rd# * a * (1.00384 + .0802894 * o * o + .000199025# * o * o * o * o + a * a * (.0998909# - .02855 * o * o - a * a * .0491032))
GOTO 2900

REM More Generalized Aitoff-Wagner
2530 ap = k1 * aa: op = k3 * ao
wx = SIN(.5 * dr# * op) * COS(dr# * ap)
wy = SIN(dr# * ap)
GOSUB 3020
az = w
wa = ap
wo = .5 * op
v1 = 0: v2 = -90: v3 = 0
GOSUB 3300
ra = 90 - ca
x2 = 2 * ra * COS(dr# * az) * k4
y2 = ra * SIN(dr# * az) * k2 * k5
GOTO 2900

REM Imitation Canters' Minimum-Error Polyconic
2540 a = aa / 30: o = ao / 30
x2 = 30 * o * (1 - .0557 * a * a - .00031 * a * a * o * o - .000016 * a * a * a * a * o * o)
y2 = 30 * a * (1.159 + o * o * (.0117 - .00106 * a * a) + .00001 * a * a * o * o * o * o)
GOTO 2900

REM Canters' Minimum-Error Polyconic
2550 a = dr# * aa: o = dr# * ao
x2 = 67.581717# * o * (.8478 - .1702 * a * a - .0062 * a * a * o * o + .0036 * a * a * a * a)
y2 = 67.581717# * a * (.9964 + o * o * (.0313 - .01 * a * a) + .0009 * o * o * o * o)
GOTO 2900

REM Canters' Minimum-Error Polyconic II
2560 a = dr# * aa: o = dr# * ao
x2 = 61.57526 * o * (.9305 - .1968 * a * a - .0067 * a * a * o * o + .0076 * a * a * a * a)
y2 = 61.57526 * a * (.9305 + o * o * (.0394 - .0115 * a * a) + .0005 * o * o * o * o)
GOTO 2900

2900 GOSUB 4000
RETURN

REM mathematical subroutines

REM arcsine
3000 IF ABS(w) > .7071 THEN 3010
w = rd# * ATN(w / SQR(1 - w * w))
RETURN
3010 w = SGN(w) * (90 - rd# * ATN(SQR(1 - w * w) / ABS(w)))
RETURN

REM atan2
3020 IF ABS(wy) < ABS(wx) THEN 3030
IF wy = 0 THEN w = 0: RETURN
w = 90 - rd# * ATN(wx / wy)
IF wy < 0 THEN w = w - 180: IF w < -180 THEN w = w + 360
RETURN
3030 w = rd# * ATN(wy / wx)
IF wx < 0 THEN w = w + 180: IF w > 180 THEN w = w - 360
RETURN

REM arithmetic-geometric mean routine
REM calculates F(w|m) where m = k^2 = (sin(alpha))^2

3100 a = 1: b = SQR(1 - m): c = SQR(m)
t1 = 1
3110 t2 = (rd# * w) / 90: IF t2 = INT(t2) THEN wa = w: GOTO 3120
wa = ATN(b * TAN(w) / a)
IF wa < 0 THEN wa = wa + pi#
wa = wa + pi# * INT(w / pi#)
3120 w = w + wa
c = .5 * (a - b): aa = .5 * (a + b): b = SQR(a * b): a = aa
t1 = t1 + t1
IF ABS(c) > .0001 THEN 3110
w = w / (t1 * a)
RETURN

REM Elliptic integral subroutine
REM calculates F( re + i*im | mm )

3200 sx = SGN(re): sy = SGN(im)
zr = ABS(re): zi = ABS(im)
t2 = EXP(zi)
t2 = .5 * (t2 - 1 / t2)
IF zr = 0 THEN x = 0: mu = ATN(t2): GOTO 3277
t3 = 1 / SIN(zr): t1 = COS(zr) * t3: t1 = t1 * t1
b = -(t1 + mm * (t2 * t2 * t3 * t3 + 1) - 1)
c = (mm - 1) * t1
t2 = .5 * (-b + SQR(b * b - 4 * c))
t3 = t2 / t1: IF t3 < 1 THEN t3 = 1
mu = ATN(SQR((t3 - 1) / mm))
lm = ATN(1 / SQR(t2))
w = lm: m = mm: GOSUB 3100: x = w
3277 w = mu: m = 1 - mm: GOSUB 3100: y = w
x = x * sx
y = y * sy
RETURN

REM Coordinate conversion as a subroutine
3300 IF v2 = 0 THEN 3350
ja = wa
jo = wo - v1
IF jo < -180 THEN jo = jo + 360
IF jo > 180 THEN jo = jo - 360

REM Convert to Cartesian coordinates
cy = SIN(dr# * ja)
vv = COS(dr# * ja)
cx = COS(dr# * jo) * vv
cz = SIN(dr# * jo) * vv

REM This matrix multiplication performs the actual rotation
vc = COS(dr# * v2)
vs = SIN(dr# * v2)
qy = cy * vc - cx * vs
qx = cx * vc + cy * vs
qz = cz

rr = SQR(qx * qx + qz * qz)
IF rr > .005 THEN ca = rd# * ATN(qy / rr): GOTO 3310
ca = 90: IF qy < 0 THEN ca = -90
co = 0
REM Near the poles, computing longitude is meaningless
GOTO 3340

3310 IF ABS(qz) > ABS(qx) THEN 3320
co = rd# * ATN(qz / qx)
IF qx < 0 THEN co = co + 180
GOTO 3340

3320 co = 90 - rd# * ATN(qx / qz)
IF qz < 0 THEN co = co - 180

3340 co = co + v3
GOTO 3370

3350 co = (wo - v1) + v3
ca = wa

3370 IF co < -180 THEN co = co + 360
IF co > 180 THEN co = co - 360

RETURN

REM Symmetric interruption input subroutine
3800 ip% = 1
INPUT "Number of gores? "; gn
IF gn < 1 THEN 3850
ni(0) = -180
si(0) = -180
ht = 180 / gn
th = ht + ht
sm = ht - 180
ir = th - 180
FOR ic% = 1 TO gn
ns(ic%) = sm
ss(ic%) = sm
ni(ic%) = ir
si(ic%) = ir
sm = sm + th
ir = ir + th
NEXT ic%
RETURN

3850 ic% = 1
PRINT "First interruption, in map coordinates, at -180."
ni(0) = -180
si(0) = -180
PRINT "Final interruption will be 180, enter when finished."
PRINT "Enter 0, then 180, for normal uninterrupted projection."
3860 INPUT "Standard meridian? "; ns(ic%)
ss(ic%) = ns(ic%)
INPUT "Interruption? "; ni(ic%)
si(ic%) = ni(ic%)
ic% = ic% + 1
IF ni(ic% - 1) < 180 THEN 3860
RETURN

REM Interruption input subroutine
3900 ip% = 1
PRINT "Northern Hemisphere"
ic% = 1
PRINT "First interruption, in map coordinates, at -180."
ni(0) = -180
PRINT "Final interruption will be 180, enter when finished."
PRINT "Enter 0, then 180, for normal uninterrupted projection."
3910 INPUT "Standard meridian? "; ns(ic%)
INPUT "Interruption? "; ni(ic%)
ic% = ic% + 1
IF ni(ic% - 1) < 180 THEN 3910
PRINT "Southern Hemisphere"
ic% = 1
PRINT "First interruption, in map coordinates, at -180."
si(0) = -180
PRINT "Final interruption will be 180, enter when finished."
3980 INPUT "Standard meridian? "; ss(ic%)
INPUT "Interruption? "; si(ic%)
ic% = ic% + 1
IF si(ic% - 1) < 180 THEN 3980
RETURN

REM line drawing and point plotting subroutines

4000 IF pp% = 1 THEN pv% = 1

REM PRINT "Attempting to plot"; pp%; pv%

IF pv% = 0 THEN pv% = 1: GOTO 4020
GOSUB 4100
4020 x1 = x2
y1 = y2
RETURN

4100 IF rt% = 0 THEN 4110

x4 = x2 * rc + y2 * rs
y4 = y2 * rc - x2 * rs

x2% = hw% + INT(sf * x4 + hd + .5)
y2% = hh% + INT(sf * y4 + vd + .5)
IF pp% = 1 THEN 4600

x3 = x1 * rc + y1 * rs
y3 = y1 * rc - x1 * rs

x1% = hw% + INT(sf * x3 + hd + .5)
y1% = hh% + INT(sf * y3 + vd + .5)
GOTO 4120

4110 x2% = hw% + INT(sf * x2 + hd + .5)
y2% = hh% + INT(sf * y2 + vd + .5)
IF pp% = 1 THEN 4600
x1% = hw% + INT(sf * x1 + hd + .5)
y1% = hh% + INT(sf * y1 + vd + .5)

REM print x1,y1, x2,y2
REM print x1%,y1%, x2%,y2%

4120 xd% = ABS(x1% - x2%)
xs% = SGN(x2% - x1%)
yd% = ABS(y1% - y2%)
ys% = SGN(y2% - y1%)
IF yd% > xd% THEN 4500
IF xd% = 0 THEN 4900

y = y1%
yi = (y2% - y1%) / xd%
FOR x% = x1% TO x2% STEP xs%
y% = INT(y + .5)
GOSUB 5000
y = y + yi
NEXT x%
RETURN

4500 x = x1%
xi = (x2% - x1%) / yd%
FOR y% = y1% TO y2% STEP ys%
x% = INT(x + .5)
GOSUB 5000
x = x + xi
NEXT y%
RETURN

4600 ii% = 1

REM PRINT "Plotting "; x2%; y2%; mg; r1(mg, ii%); ii%

4610 IF r1(mg, ii%) > 99 THEN 4690
x% = x2% + r1(mg, ii%)
y% = y2% + r2(mg, ii%)

REM PRINT "Point "; x%; y%

GOSUB 5000
ii% = ii% + 1
GOTO 4610
4690 RETURN

REM Plot a single point
4900 x% = x1%: y% = y1%: GOSUB 5000
RETURN

5000 IF x% < 1 THEN RETURN
IF x% >= w% THEN RETURN
IF y% < 1 THEN RETURN
IF y% >= h% THEN RETURN
xb% = x% \ 16
mi% = ma%(x% - 16 * xb%)
im%(xb% + 1, h% - y%) = im%(xb% + 1, h% - y%) OR mi%
RETURN

REM .xbm file writing subroutines

10000 IF wi% < 0 THEN 10010
os$ = os$ + " " + la$ + ","
10010 la$ = a$
wi% = wi% + 1
IF wi% < 12 THEN 10020
PRINT #2, os$
os$ = "  "
wi% = 0
10020 RETURN

10100 os$ = os$ + " " + la$ + " } ;"
PRINT #2, os$
RETURN

REM binary values for creating graphics file
DATA 256,512,1024,2048,4096,8192,16384,-32768,1,2,4,8,16,32,64,128

REM displacements for plotting points
DATA -3,0, -3,1, -2,2, -1,3, 0,3, 1,3, 2,2, 3,1, 3,0, 3,-1
DATA 2,-2, 1,-3, 0,-3, -1,-3, -2,-2, -3,-1
DATA -2,0, 2,0, 0,-2, 0,2
DATA -1,0, 1,0, 0,-1, 0,1, 0,0, 100,100

DATA -3,0, -2,-1, -1,-2, 0,-3, 1,-2, 2,-1, 3,0, 2,1, 1,2, 0,3
DATA -1,2, -2,1, -2,0, -1,0, 0,0, 100, 100

DATA -2,0, -2,1, -1,2, 0,2, 1,2, 2,1, 2,0, 2,-1, 1,-2, 0,-2, -1,-2
DATA -2,-1, 0,1, 0,0, 100,100

DATA -1,0, 1,0, 0,-1, 0,1, 0,0, 100,100

DATA 0,0, 100,100