REM Copyright (C) 2001, 2019 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#

90 PRINT "Select projection or group of projections (0 to finish):"

PRINT " 1) Cylindrical projections"
PRINT " 2) Azimuthal projections (in hemispheres unless otherwise noted)"
PRINT " 3) Conic projections"
PRINT " 4) Pseudocylindrical and pseudoconic projections"
PRINT " 5) Polyconic projections"
PRINT " 6) Conventional projections"
PRINT " 7) Other conformal projections"
PRINT " 8) Other equal-area projections"
PRINT " 9) Other projections"
PRINT "99) List projection groups again"

100 INPUT "Projection: ", pj%
IF pj% = 0 THEN 3990
IF pj% = 99 THEN 90
IF pj% > 99 THEN 700

ON pj% GOTO 110, 120, 130, 140, 150, 160, 170, 180, 190

110 PRINT "101) Mercator             102) Gall                   103) Miller"
PRINT "104) Plate Carree         105) Cylindrical Equal-Area"
GOTO 100

120 PRINT "106) Gnomonic (incomplete)"
PRINT "108) Azimuthal Equidistant"
PRINT "107) Stereographic        109) Lambert Azimuthal Equal-area"
PRINT "110) Orthographic         111) Azimuthal Equidistant (single circle)"
PRINT "164) Stereographic (single circle)"
PRINT "165) Airy's Minimum-Error Azimuthal Projection (single circle)"
GOTO 100

130 PRINT "112) Lambert conformal    113) Simple conic           114) Albers equal-area"
PRINT "162) Simple Conic with Two Standard Parallels"
GOTO 100

140 PRINT "115) Sinusoidal           116) Bonne's                117) Mollweide's"
PRINT "118) Eckert IV            119) MacBryde-Thomas Flat Polar Quartic"
PRINT "120) Adams' Equal-Area    121) Parabolic stereographic"
PRINT "122) Loximuthal"
PRINT "159) Bromley-Mollweide    160) Modified Gall        161) Equal Earth"
GOTO 100

150 PRINT "123) Polyconic            143) Rectangular Polyconic"
GOTO 100

160 PRINT "124) Globular             125) Van der Grinten's      126) Van der Grinten IV"
PRINT "127) Savard Egg           128) Ster Modified VdG IV   129) Globular (whole world)"
PRINT "146) New Conventional"
PRINT "147) Conventional II      148) Modified VdG IV        150) Ginzoid"
PRINT "169) Expanding Parabolic  172) Henri Bouthillier de Beaumont's Projection"
GOTO 100

170 PRINT "130) Lagrange's conformal 131) August's conformal     132) Cosine-based"
PRINT "133) Eisenlohr            134) Guyou's doubly-periodic"
PRINT "135) World in rectangle   136) World in diamond"
PRINT "144) Four-square world    145) Quadfoil 1"
PRINT "149) Corrected Adams Ellipse                              157) Dixon elliptic"
PRINT "168) Generalized Lagrange Conformal"
GOTO 100

180 PRINT "137) Hammer-Aitoff        138) Generalized Hammer-Aitoff"
PRINT "158) Generalized Aitoff Equal-Area"
PRINT "166) Dietrich-Kitada          171) Generalized Hammer-Wagner"
GOTO 100

190 PRINT "139) Aitoff-Wagner        140) Aitoff projection    141) Winkel's Tripel"
PRINT "142) Pole-line Aitoff"
PRINT "151) Larrivee             152) Laskowski Tri-Optimal"
PRINT "153) Generalized Aitoff   154) Approximate Canters"
PRINT "155) Canters' Minimum-Error Polyconic              156) Canters' Equal-Scale"
PRINT "163) Regional, Kite, Lotus                         167) Modified Regional"
GOTO 100

700 px% = (pj% - 101) \ 10
py% = (pj% - 100) - 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 900
tc = COS(dr# * a2)
ts = SIN(dr# * a2)

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

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

ON pa% GOTO 950, 960, 970, 980

950 ON pb% GOTO 1010, 1020, 1030, 1040, 1050, 1060, 1070, 1080, 1090, 1100, 1110, 1120, 1130, 1140, 1150, 1160, 1170, 1180, 1190, 1200
960 ON pb% GOTO 1210, 1220, 1230, 1240, 1250, 1260, 1270, 1280, 1290, 1300, 1310, 1320, 1330, 1340, 1350, 1360, 1370, 1380, 1390, 1400
970 ON pb% GOTO 1410, 1420, 1430, 1440, 1450, 1460, 1470, 1480, 1490, 1500, 1510, 1520, 1530, 1540, 1550, 1560, 1570, 1580, 1590, 1600
980 ON pb% GOTO 1610, 1620, 1630, 1640, 1650, 1660, 1670, 1680, 1690, 1700, 1710, 1720

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 12) Check as above and for crossing 180 in northern hemisphere
REM 13) Interrupted projection with central band

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
REM 14) Octahedron unfolded
REM 15) Projection-native with mandatory dual interruptions (i.e., regional)

1010 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 3665
1020 tp% = 2: bo% = 3
k1 = rd# * (1 + SQR(2))
ar = 180 / k1
GOTO 3665
1030 tp% = 2: bo% = 3: ar = 1.363886
GOTO 3665
1040 tp% = 2: bo% = 3: ar = 2
GOTO 3665
1050 INPUT "Standard parallel: ", sp
k1 = COS(dr# * sp): k2 = 1 / (k1 * k1)
tp% = 2: bo% = 3: ar = k1 * k1 * 3.141592
GOTO 3665
1060 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 3665
1070 tp% = 3: bo% = 2: ar = 2
GOTO 3665
1080 tp% = 3: bo% = 2: ar = 2
GOTO 3665
1090 tp% = 3: bo% = 2: ar = 2
GOTO 3665
1100 tp% = 3: bo% = 2: ar = 2
GOTO 3665
1110 tp% = 1: bo% = 1: ar = 1
GOTO 3665
1120 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 3665
1130 INPUT "Standard parallel: ", sp
tp% = 2: bo% = 9: ar = 2
k1 = SIN(dr# * sp)
k2 = rd# / TAN(dr# * sp)
GOTO 3665
1140 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 3665
1150 sp% = 0: GOSUB 7900
tp% = 8: bo% = 6: ar = 2
GOTO 3665
1160 PRINT "(If integer, interruptions allowed)"
INPUT "Standard parallel: ", sp
sp% = sp
IF sp% = sp THEN 1165
sp% = 0: ip% = 1
ni(0) = -180: ni(1) = 180: ns(1) = 0
si(0) = -180: si(1) = 180: ss(1) = 0
GOTO 1167
1165 GOSUB 7900
1167 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 3665
1170 sp% = 0: GOSUB 7900
tp% = 8: bo% = 6: ar = 2
REM Border was 1 when not interrupted
GOTO 3665
1180 sp% = 0: GOSUB 7900
tp% = 8: bo% = 8: ar = 2
REM Border was 5 when not interrupted
GOTO 3665
1190 sp% = 0: GOSUB 7900
tp% = 8: bo% = 8: ar = 2.221442
GOTO 3665
1200 sp% = 0: GOSUB 7900
tp% = 8: bo% = 6: ar = 2.221442
GOTO 3665
1210 sp% = 0: GOSUB 7900
tp% = 8: bo% = 6: ar = 1.570796
GOTO 3665
1220 INPUT "Standard parallel? ", sp
k1 = rd# * LOG((1 / COS(dr# * sp)) + TAN(dr# * sp))
tp% = 2: bo% = 6: ar = 2
GOTO 3665
1230 GOSUB 7800
tp% = 8: bo% = 6: ar = 2
GOTO 3665
1240 tp% = 4: bo% = 2: ar = 2
GOTO 3665
1250 tp% = 2: bo% = 1: ar = 1
GOTO 3665
1260 tp% = 2: bo% = 6: ar = 1.6
GOTO 3665
1270 tp% = 2: bo% = 1: ar = 1.5
GOTO 3665
1280 tp% = 2: bo% = 6: ar = 1.63
GOTO 3665
1290 tp% = 2: bo% = 6: ar = 1.6
GOTO 3665
1300 tp% = 2: bo% = 1: ar = 1
GOTO 3665
1310 tp% = 2: bo% = 6: ar = 1.42
GOTO 3665
1320 tp% = 2: bo% = 6: ar = 1.5
GOTO 3665
1330 tp% = 2: bo% = 6: ar = 1.5
GOTO 3665
1340 tp% = 9: bo% = 10: ar = 2
GOTO 3665
1350 INPUT "Parameter (m, k^2, (sin(alpha))^2) (.5: square): ", k1
w = .5 * pi#: m = k1: GOSUB 7100: wi = w
w = .5 * pi#: m = 1 - k1: GOSUB 7100: ht = w
tp% = 2: bo% = 3: ar = wi / ht
k2 = 180 / wi
PRINT "Aspect ratio will be: "; ar
GOTO 3665
1360 tp% = 2: bo% = 11: ar = 1
GOTO 3665
1370 tp% = 2: bo% = 1: ar = 2
GOTO 3665
1380 PRINT "Stretch factor (or its reciprocal)"
INPUT ".875: Briesemeister, 1.5: 3:2 ellipse: ", p
IF p < 1 THEN sh = p: ex = 1 / p: GOTO 1385
sh = 1 / p: ex = p
1385 tp% = 2: bo% = 6: ar = 2 * sh
GOTO 3665
1390 tp% = 2: bo% = 9: ar = 1.6
GOTO 3665
1400 tp% = 2: bo% = 1: ar = 2
GOTO 3665
1410 tp% = 2: bo% = 8: ar = 1.766044
GOTO 3665
1420 INPUT "Compression factor (or its reciprocal): ", p
IF p < 1 THEN sh = p: ex = 1 / p: GOTO 1425
sh = 1 / p: ex = p
1425 tp% = 10: th% = 16: bo% = 9: ar = 2 * sh
GOTO 3665
1430 INPUT "Standard parallel: ", sp
GOSUB 7800
k1 = SIN(dr# * sp)
tp% = 8: bo% = 6: ar = 2
GOTO 3665
1440 tp% = 11: bo% = 3: ar = 1
GOTO 3665
1450 tp% = 11: bo% = 7: ar = 1
GOTO 3665
1460 tp% = 10: th% = 15: bo% = 13: ar = 1.333333
GOTO 3665
1470 tp% = 10: th% = 15: bo% = 13: ar = 1.333333
GOTO 3665
1480 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 3665
1490 tp% = 2: bo% = 1
INPUT "Parameter (m, k^2, (sin(alpha))^2) (.5: square): ", k1
w = .5 * pi#: m = k1: GOSUB 7100: wi = w
w = .5 * pi#: m = 1 - k1: GOSUB 7100: ht = w
ag = wi / ht
k2 = 180 / wi
PRINT "Aspect ratio of underlying rectangle will be: "; ag
k4 = ag
REM sin(z) = sin(x) * cosh(y) + i * cos(x) * sinh(y)
REM interested in x = 0 or x = pi/2
REM interested in y = pi/(2*ag)
ht = pi# / (ag + ag)
xh = EXP(ht)
he = xh - 1 / xh
we = xh + 1 / xh
k3 = 360 / we
ar = we / he
PRINT "Aspect ratio of ellipse will be: "; ar
GOTO 3665
1500 GOSUB 7800
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 3665
1510 tp% = 10: th% = 8: bo% = 9: ar = 1.3
GOTO 3665
1520 tp% = 10: th% = 24: bo% = 9: ar = 1.8
GOTO 3665
1530 INPUT "Latitude compression factor (or its reciprocal): ", p
IF p < 1 THEN k1 = p: k2 = 1 / p: GOTO 1533
k1 = 1 / p: k2 = p
1533 PRINT "(If negative, relative to Azimuthal Equidistant rather than Aitoff)"
INPUT "Longitude compression factor (or its reciprocal): ", p
IF p < 0 THEN p = -p: GOTO 1534
IF p < 1 THEN k3 = p: k4 = 1 / p: GOTO 1535
k3 = 1 / p: k4 = p: GOTO 1535
1534 IF p < 1 THEN k3 = 2 * p: k4 = 1 / k3: GOTO 1535
k4 = 0.5 * p: k3 = 1 / k4
1535 INPUT "Vertical stretch: ", p
k5 = p: IF p < 1 THEN k5 = 1 / p
tp% = 10: th% = 16: bo% = 9: ar = 2 * k1 / k5
GOTO 3665
1540 tp% = 10: th% = 16: bo% = 9: ar = 1.5
GOTO 3665
1550 tp% = 10: th% = 16: bo% = 9: ar = 1.5
GOTO 3665
1560 tp% = 10: th% = 16: bo% = 9: ar = 1.5
GOTO 3665
1570 tp% = 12: bo% = 14: ar = 1.7320508
k1 = 103.92305
k2 = 90 / 1.1129127
k3 = COS(dr# * 30)
GOTO 3665
1580 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 1585
sh = 1 / p: ex = p
1585 tp% = 10: th% = 16: bo% = 9: ar = 2 * sh
GOTO 3665
1590 sp% = 0: GOSUB 7900
tp% = 8: bo% = 6: ar = 2.4674011
k1 = 0.81056947
k2 = 90 * k1
GOTO 3665
1600 tp% = 2: bo% = 8
k1 = rd# * (1 + SQR(2))
ar = 180 / k1
GOTO 3665
1610 tp% = 2: bo% = 8: ar = 2.205458
k1 = SQR(3) / 2.0
k2 = 4.020792 / (4.0 * k1)
GOTO 3665
1620 INPUT "Standard parallels: ", s1, s2
tp% = 2: bo% = 9: ar = 2
IF s2 < s1 THEN t = s2: s2 = s1: s1 = t
IF s2 - s1 < .001 THEN 1625
c1 = COS(dr# * s1)
c2 = COS(dr# * s2)
k2 = (c2 / (c1 - c2)) * (s2 - s1)
k1 = rd# * c2 / k2
GOTO 3665
1625 s1 = .5 * (s1 + s2)
k1 = SIN(dr# * s1)
k2 = rd# / TAN(dr# * s1)
GOTO 3665
1630 GOSUB 7900
PRINT "23.5, 66.5 in original"
INPUT "Standard parallels: ", s1, s2
tp% = 13: bo% = 15: ar = 2
IF s2 < s1 THEN t = s2: s2 = s1: s1 = t
s2% = s2 + 1
s1% = s1 - 1
IF s2 - s1 < .001 THEN 1635
c1 = COS(dr# * s1)
c2 = COS(dr# * s2)
k2 = (c2 / (c1 - c2)) * (s2 - s1)
k1 = rd# * c2 / k2
PRINT "Angular scale factor: "; k1
GOTO 3665
1635 s1 = .5 * (s1 + s2)
k1 = SIN(dr# * s1)
k2 = rd# / TAN(dr# * s1)
GOTO 3665
1640 tp% = 1: bo% = 3: ar = 1
INPUT "Angular limit of projection: ", rg
u = 90 * TAN(dr# * rg / 2)
k1 = 180 / u
GOTO 3665
1650 tp% = 1: bo% = 3: ar = 1
INPUT "Limit to which error is to be minimized: ", ch
INPUT "Angular limit of projection: ", rg
ag = dr# * 0.5 * ch
sz = SIN(ag)
cz = COS(ag)
cc = cz / sz
k1 = 0 - (cc * cc * LOG(cz))
aq = dr# * 0.5 * rg
sr = SIN(aq)
cr = COS(aq)
ck = cr / sr
u = (0 - (ck * LOG(cr))) + k1 / ck
k2 = 180 / u
GOTO 3665
1660 tp% = 2: bo% = 6: ar = 1.933
k2 = 1.21262
k1 = k2 * pi# * pi# / 8
GOTO 3665
1670 tp% = 13: bo% = 8: ar = 2
GOSUB 7850
REM No interruptions on the northern side
REM doing this allows use of existing border type 8
ns%(1) = 0
ni%(1) = 180
INPUT "Northern standard parallel: ", s1
INPUT "Latitude where interruptions begin: ", s2
INPUT "Width of transition band: ", tr
s3 = s2 + 0.5 * tr
s4 = s2 + tr
s2% = s2 + 1
c1 = COS(dr# * s1)
c3 = COS(dr# * s3)
k2 = (c3 / (c1 - c3)) * (s3 - s1)
k1 = rd# * c3 / k2
PRINT "Angular scale factor: "; k1
k3 = SIN(dr# * s1)
k4 = (k2 + (s3 - s1)) / (TAN(.5 * dr# * (90 + s1)) ^ k3)
GOTO 3665
1680 INPUT "Angular range on Stereographic to use: ", x
k1 = x/360
tp% = 2: bo% = 6: ar = 1
GOTO 3665
1690 INPUT "Center excess in degrees (e.g. 7.5):"; k1
INPUT "Edge excess in degrees (e.g. 20):"; k2 
tp% = 2: bo% = 8: ar = 2
GOTO 3665
1700 INPUT "Ellipse height excess in degrees (e.g. 7.5):"; k2
INPUT "Center intersection excess in degrees (e.g. 4):"; k1
INPUT "Center height excess in degrees (e.g. 30):"; k3
tp% = 2: bo% = 8: ar = 2
GOTO 3665
1710 PRINT "Latitude compression factor (or its reciprocal)"
INPUT ".9: original Aitoff, .90630779: Wagner, 1.5: 3:2 ellipse: ", p
IF p < 1 THEN sh = p: ex = 1 / p: GOTO 1712
sh = 1 / p: ex = p
1712 PRINT "Longitude compression factor (or its reciprocal)"
INPUT "2: original Hammer-Aitoff, 3: Wagner: ", q
IF q < 1 THEN hs = q: he = .5/q: GOTO 1714
hs = 1/q: he = q/2
1714 PRINT "Final vertical stretch (or its reciprocal)"
INPUT "1.2650925: Wagner: "; r
IF r < 1 THEN vh = r: ve = 1/r: GOTO 1718
vh = 1/r: ve = r
1718 tp% = 10: th% = 16: bo% = 9: ar = 2 * sh
GOTO 3665
1720 tp% = 2: bo% = 6: ar = 1.6
GOTO 3665

3665 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 3670
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 3670
rt% = 1
rc = COS(dr# * a)
rs = SIN(dr# * a)

3670 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 3675
LINE INPUT "Input file name? "; fj$
OPEN fj$ FOR INPUT AS #1

3675 IF (ac% AND 2) = 0 THEN 3680

ac% = ac% OR 24

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

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

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

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

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

REM draw map border

3700 pv% = 0

REM PRINT "Drawing border of type "; bo%

IF bo% = 0 THEN 3793
ON bo% GOTO 3705, 3710, 3715, 3720, 3725, 3730, 3735, 3740, 3745, 3750, 3755, 3760, 3765, 3770, 3775

3705 FOR i = 0 TO 360
    ag = dr# * i
    x2 = 180 * COS(ag)
    y2 = 180 * SIN(ag) / ar
    GOSUB 8000
NEXT i
GOTO 3793

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

3715 xb = 180: yb = 180 / ar
x2 = xb: y2 = yb
GOSUB 8000
x2 = xb: y2 = -yb
GOSUB 8000
x2 = -xb: y2 = -yb
GOSUB 8000
x2 = -xb: y2 = yb
GOSUB 8000
x2 = xb: y2 = yb
GOSUB 8000
GOTO 3793

3720 REM

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

3730 IF ip% = 1 THEN 3731
ix = -180
FOR iy = -90 TO 90
    aa = iy: ao = ix
    GOSUB 4080
NEXT iy
ix = 180
FOR iy = 90 TO -90 STEP -1
    aa = iy: ao = ix
    GOSUB 4080
NEXT iy
GOTO 3793

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

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

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

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

ic% = ic% + 1: GOTO 3732

REM Southern hemisphere

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

IF si(ic%) > 179.9 THEN 3793

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

ic% = ic% + 1: GOTO 3733

3735 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 8000
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 8000
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 8000
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 8000
NEXT ix
GOTO 3793

3740 IF ip% = 1 THEN 3741
ix = -180
FOR iy = -90 TO 90
    aa = iy: ao = ix
    GOSUB 4080
NEXT iy
ix = 180
FOR iy = 90 TO -90 STEP -1
    aa = iy: ao = ix
    GOSUB 4080
NEXT iy
ao = -180: aa = -90
GOSUB 4080
GOTO 3793

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

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

IF ni(ic%) > 179.9 THEN 3743

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

ic% = ic% + 1: GOTO 3742

REM Southern hemisphere

3743 ic% = 1
pv% = 0

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

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

IF si(ic%) > 179.9 THEN 3793

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

ic% = ic% + 1: GOTO 3744

3745 ix = -180
FOR iy = -90 TO 90
    aa = iy: ao = ix
    GOSUB 4080
NEXT iy
iy = 90
FOR ix = -179 TO 180
    aa = iy: ao = ix
    GOSUB 4080
NEXT ix
ix = 180
FOR iy = 89 TO -90 STEP -1
    aa = iy: ao = ix
    GOSUB 4080
NEXT iy
iy = -90
FOR ix = 179 TO -180 STEP -1
    aa = iy: ao = ix
    GOSUB 4080
NEXT ix
GOTO 3793

3750 xb = 180: yb = 180 / ar
x2 = xb: y2 = yb
GOSUB 8000
x2 = xb: y2 = -yb
GOSUB 8000
x2 = -xb: y2 = -yb
GOSUB 8000
x2 = -xb: y2 = yb
GOSUB 8000
x2 = xb: y2 = yb
GOSUB 8000
pv% = 0
x2 = 0: y2 = yb
GOSUB 8000
x2 = 0: y2 = -yb
GOSUB 8000
GOTO 3793

3755 xb = 180: yb = 180 / ar
x2 = 0: y2 = yb
GOSUB 8000
x2 = xb: y2 = 0
GOSUB 8000
x2 = 0: y2 = -yb
GOSUB 8000
x2 = -xb: y2 = 0
GOSUB 8000
x2 = 0: y2 = yb
GOSUB 8000
GOTO 3793

3760 x2 = 0: y2 = 90
GOSUB 8000
FOR ix = -180 TO 180
    x2 = 180 * SIN(dr# * ix * k1)
    y2 = 90 - 180 * COS(dr# * ix * k1)
    GOSUB 8000
NEXT ix
x2 = 0: y2 = 90
GOSUB 8000
GOTO 3793

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

3770 x2 = -180: y2 = .5 * k1
GOSUB 8000
x2 = -90: y2 = k1
GOSUB 8000
x2 = 0: y2 = .5 * k1
GOSUB 8000
x2 = 90: y2 = k1
GOSUB 8000
x2 = 180: y2 = .5 * k1
GOSUB 8000
x2 = 90: y2 = 0
GOSUB 8000
x2 = 90: y2 = -k1
GOSUB 8000
x2 = 0: y2 = -.5 * k1
GOSUB 8000
x2 = -90: y2 = -k1
GOSUB 8000
x2 = -90: y2 = 0
GOSUB 8000
x2 = -180: y2 = .5 * k1
GOSUB 8000
GOTO 3793

3775 ix = -180
ao = ix
vo = ix - ss(1)
dx = ss(1)
FOR iy = -90 TO s1%
    aa = iy
    GOSUB 4800
NEXT iy
ao = ix
FOR iy = s1% + 1 TO s2% - 1
    aa = iy
    GOSUB 4800
NEXT iy
ao = ix
vo = ix - ns(1)
dx = ns(1)
FOR iy = s2% TO 90
    aa = iy
    GOSUB 4800
NEXT iy

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

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

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

ic% = ic% + 1: GOTO 3777

REM Southern hemisphere

3778 ix = si(ic%)
ao = ix
vo = ix - ss(ic%)
dx = ss(ic%)
FOR iy = -90 TO s1%
    aa = iy
    GOSUB 4800
NEXT iy
aa = s1
GOSUB 4800

IF si(ic%) > 179.9 THEN 3779

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

ic% = ic% + 1: GOTO 3778

3779 ao = 180
FOR iy = s1% + 1 TO s2% - 1
    aa = iy
    GOSUB 4800
NEXT iy
aa = s2
GOSUB 4800
GOTO 3793

REM draw latitude and longitude lines

3793 IF (ac% AND 16) = 0 THEN 3795

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

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

3795 IF (ac% AND 8) = 0 THEN 3797

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

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

3797 IF (ac% AND 64) = 0 THEN 3800

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 4000
    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 4000
    NEXT ix
NEXT iy

REM PRINT "Beginning actual draw."

3800 IF (ac% AND 1) = 0 THEN 3820
pv% = 0
pp% = 0

REM read map file and draw

3810 IF EOF(1) THEN 3820
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 4000

NEXT ix%

GOTO 3810

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

pv% = 0
pp% = 0

REM read map file and draw

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

REM PRINT tx$

n1 = INSTR(tx$, " ")

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

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 4000

GOTO 3825

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

REM PRINT "Plotting point layer"

pv% = 0
pp% = 1
3860 IF EOF(3) THEN 3900
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 3860
IF mg < 1 THEN mg = 1

REM PRINT la, lo, mg

GOSUB 4000

GOTO 3860

3900 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 100

3990 END

REM map drawing routine

REM coordinate transform

4000 IF ti% = 0 THEN 4050
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 4070

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 4005

REM Secondary rotation (for ease of matching)

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

qy = cy: qz = cz

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

4010 IF ABS(qz) > ABS(qx) THEN 4020
ao = rd# * ATN(qz / qx)
IF qx < 0 THEN ao = ao + 180
GOTO 4040

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

4040 ao = ao + a3
GOTO 4070

4050 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

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

4080 IF pp% = 1 THEN pv% = 1: GOTO 4800
IF (tp% <> 8) AND (tp% <> 13) AND (pv% = 0) THEN 4800
ON tp% GOTO 4110, 4120, 4130, 4140, 4150, 4160, 4170, 4180, 4190, 4200, 4210, 4220, 4230

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

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

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

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

4150 REM
4160 REM
4170 REM
GOTO 4800

REM Interrupted projection
4180 cu% = 1

REM PRINT "Was:"; aa; ao

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

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

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

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

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

4195 IF ao * po > 0 THEN 4800
pv% = 0
GOTO 4800

REM Aitoff case
4200 IF ABS(ao) < 90 THEN 4205
IF ABS(po) < 90 THEN 4205
IF ao * po > 0 THEN 4205
pv% = 0
GOTO 4800
4205 IF aa < -70 THEN 4207
IF aa > 70 THEN 4207
GOTO 4800
4207 uu = ABS(ao - po)
IF uu < th% THEN 4800
IF uu > 360 - th% THEN 4800
pv% = 0
GOTO 4800

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

REM Gnomonic or Conformal Butterfly case
4220 IF pa >= 0 THEN 4800
cu% = INT(ao / 90)
IF aa >= 0 THEN 4225
IF pv% = 1 THEN IF pc% <> cu% THEN pv% = 0
GOTO 4229
4225 IF ABS(ao) < 90 THEN 4800
IF ABS(po) < 90 THEN 4800
IF ao * po > 0 THEN 4800
pv% = 0
4229 pc% = cu%
GOTO 4800

REM Regional projection
4230 cu% = 1

REM PRINT "Was:"; aa; ao

IF aa < s1 THEN 4235
IF aa > s2 THEN 4232
GOTO 4239

4232 IF ao > ni(cu%) THEN cu% = cu% + 1: GOTO 4232
vo = ao - ns(cu%)
dx = ns(cu%)
IF pv% = 1 THEN IF pa > s2 THEN IF pc% <> cu% THEN pv% = 0
GOTO 4237

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

4237 pc% = cu%

REM PRINT "Became:"; aa; ao; vlo1; dx1

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

4800 po = ao: pa = aa

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

ON px% GOTO 4810, 4820, 4830, 4840, 4850, 4860, 4870, 4880

4810 ON py% GOTO 5010, 5020, 5030, 5040, 5050, 5060, 5070, 5080, 5090, 5100
4820 ON py% GOTO 5110, 5120, 5130, 5140, 5150, 5160, 5170, 5180, 5190, 5200
4830 ON py% GOTO 5210, 5220, 5230, 5240, 5250, 5260, 5270, 5280, 5290, 5300
4840 ON py% GOTO 5310, 5320, 5330, 5340, 5350, 5360, 5370, 5380, 5390, 5400
4850 ON py% GOTO 5410, 5420, 5430, 5440, 5450, 5460, 5470, 5480, 5490, 5500
4860 ON py% GOTO 5510, 5520, 5530, 5540, 5550, 5560, 5570, 5580, 5590, 5600
4870 ON py% GOTO 5610, 5620, 5630, 5640, 5650, 5660, 5670, 5680, 5690, 5700
4880 ON py% GOTO 5710, 5720

REM Mercator
5010 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 6900

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

REM Miller's
5030 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 6900

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

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

REM Gnomonic
5060 IF ABS(aa) < gb THEN pv% = 0: RETURN
IF aa < 0 THEN 5065
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 6900
5065 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 6900

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

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

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

REM Orthographic
5100 IF aa < 0 THEN 5105
rr = 90 * SIN(dr# * (90 - aa))
x2 = rr * SIN(dr# * ao) - 90
y2 = -rr * COS(dr# * ao)
GOTO 6900
5105 rr = 90 * SIN(dr# * (90 + aa))
x2 = 90 - rr * SIN(dr# * ao)
y2 = -rr * COS(dr# * ao)
GOTO 6900

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

REM Lambert's Conic Conformal
5120 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 6900

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

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

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

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

REM Mollweide
5170 pp = dr# * aa
k = pi# * .5 * SIN(pp)
IF aa = 0 THEN ps = 0: GOTO 5175
IF ABS(aa) = 90 THEN ps = pp: GOTO 5175
5172 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 5172
5175 x2 = vo * COS(ps) + dx
y2 = 90 * SIN(ps)
GOTO 6900

REM Eckert IV Projection
5180 pp = dr# * aa
k = (.5 * pi# + 2) * SIN(dr# * aa)
IF aa = 0 THEN ps = 0: GOTO 5185
IF ABS(aa) = 90 THEN ps = pp: GOTO 5185
5182 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 5182
5185 x2 = vo * .5 * (1 + COS(ps)) + dx
y2 = 90 * SIN(ps)
GOTO 6900

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

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

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

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

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

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

REM Parabolic case
5235 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)

5237 x2 = xx + dx
y2 = sa * yy
GOTO 6900

REM Globular
5240 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 5247
IF qa = 0 THEN x = qo: y = 0: GOTO 5247
IF qo = 0 THEN x = 0: y = qa: GOTO 5247
IF qo = 90 THEN x = 90 * COS(dr# * qa): y = 90 * SIN(dr# * qa): GOTO 5247

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

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
5247 x2 = x * sx + dx
y2 = y * sy
GOTO 6900

REM Van der Grinten
5250 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 5257
IF qa = 0 THEN x = qo: y = 0: GOTO 5257
IF qo = 0 THEN x = 0: y = ch: GOTO 5257
IF qo = 180 THEN x = SQR(32400 - eh * eh): y = eh: GOTO 5257

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

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
5257 x2 = x * sx
y2 = y * sy
GOTO 6900

REM Van der Grinten IV
5260 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 5267
IF qa = 0 THEN x = qo: y = 0: GOTO 5267
IF qo = 0 THEN x = 0: y = qa: GOTO 5267

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

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 5267

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
5267 x2 = x * sx
y2 = y * sy
GOTO 6900

REM Savard Egg
5270 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 5275
r = dl / 2 + 2 * tk / dl
y = yz + r - SQR(r * r - x * x)
GOTO 5277

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

5277 x2 = x * SGN(ao)
y2 = y * SGN(aa)
GOTO 6900

REM Modified Van der Grinten IV
5280 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 5287
IF qa = 0 THEN x = qo: y = 0: GOTO 5287
IF qo = 0 THEN x = 0: y = qa: GOTO 5287

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

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 5287

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
5287 x2 = 1.27324 * x * sx
y2 = 1.27324 * y * sy
GOTO 6900

REM Globular (whole world)
5290 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 5297
IF qa = 0 THEN x = qo: y = 0: GOTO 5297
IF qo = 0 THEN x = 0: y = qa: GOTO 5297
IF qo = 90 THEN x = 90 * COS(dr# * qa): y = 90 * SIN(dr# * qa): GOTO 5297

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

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
5297 x2 = x * sx
y2 = y * sy
GOTO 6900

REM Lagrange Conformal
5300 oo = .5 * ao
w = SQR(COS(dr# * aa)) / COS(.5 * dr# * aa)
GOSUB 7000
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 6900

REM August's Conformal Projection of the World on a Two-Cusped Epicycloid
5310 oo = .5 * ao
w = SQR(COS(dr# * aa)) / COS(.5 * dr# * aa)
GOSUB 7000
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 6900

REM Cosine-based
5320 oo = .5 * ao
w = SQR(COS(dr# * aa)) / COS(.5 * dr# * aa)
GOSUB 7000
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 6900

REM Eisenlohr
5330 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 6900

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

mm = .5: GOSUB 7200

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

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

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

mm = .5: GOSUB 7200

IF dp < 0 THEN 5347

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

5347 x2 = -48.54173 * y - 90
y2 = 48.54173 * x
GOTO 6900

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

mm = k1: GOSUB 7200

x2 = k2 * x
y2 = k2 * y

GOTO 6900

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

v1 = 90: v2 = 90: v3 = 90
wa = qa: wo = qo
GOSUB 7300

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

mm = 1 - k1: GOSUB 7200

x2 = -k2 * y
y2 = k2 * x

GOTO 6900

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

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

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

mm = .5: GOSUB 7200

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

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

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

mm = .5: GOSUB 7200

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

REM Hammer-Aitoff Projection
5370 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 6900

REM Generalized Hammer-Aitoff
5380 di = SQR(1 + COS(dr# * aa) * COS(dr# * .5 * ao))
y2 = 90 * ex * SIN(dr# * aa) / di
x2 = 180 * COS(dr# * aa) * SIN(dr# * .5 * ao) / di
GOTO 6900

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

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

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

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

REM Rectangular Polyconic
5430 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 5432
co = 2 * rd# * TAN(dr# * .5 * mo * k1) / k1
5432 IF ABS(ma) < 3 THEN 5435
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 5437

REM Parabolic case
5435 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)

5437 x2 = xx * so + dx
y2 = yy * sa
GOTO 6900

REM Four Corners
5440 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 7300

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

mm = .5: GOSUB 7200

x0 = x
y0 = y
GOTO 5442

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

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

mm = .5: GOSUB 7200

x0 = -y
y0 = x

5442 ON dq + 1 GOTO 5448, 5445, 5446, 5447

5445 x2 = 48.54173 * x0 - 90
y2 = 48.54173 * y0 - 90
GOTO 6900

5446 x2 = 90 - 48.54173 * y0
y2 = 48.54173 * x0 - 90
GOTO 6900

5447 x2 = 90 - 48.54173 * x0
y2 = 90 - 48.54173 * y0
GOTO 6900

5448 x2 = 48.54173 * y0 - 90
y2 = 90 - 48.54173 * x0
GOTO 6900

REM Interrupted conformal
5450 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 7300

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

mm = .5: GOSUB 7200

x0 = x
y0 = y
GOTO 5452

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

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

mm = .5: GOSUB 7200

x0 = -y
y0 = x

5452 ON dq + 1 GOTO 5458, 5455, 5456, 5457

5455 x = .2696763 * x0 - .5
y = .2696763 * y0 - .5
GOTO 5459

5456 x = .5 - .2696763 * y0
y = .2696763 * x0 - .5
GOTO 5459

5457 x = .5 - .2696763 * x0
y = .5 - .2696763 * y0
GOTO 5459

5458 x = .2696763 * y0 - .5
y = .5 - .2696763 * x0
GOTO 5459

5459 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 6900

REM Savard Conventional
5460 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 5467
IF qa = 0 THEN x = qo: y = 0: GOTO 5467
IF qo = 0 THEN x = 0: y = qa: GOTO 5467

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

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 7000
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
5467 x2 = x * sx
y2 = y * sy
GOTO 6900

REM Savard Conventional II
5470 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 5477
IF qa = 0 THEN x = qo: y = 0: GOTO 5477
IF qo = 0 THEN x = 0: y = qa: GOTO 5477

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

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 7000
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
5477 x2 = x * sx
y2 = y * sy
GOTO 6900

REM Modified Van der Grinten IV
5480 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 5487
IF qa = 0 THEN x = qo: y = 0: GOTO 5487
IF qo = 0 THEN x = 0: y = qa: GOTO 5487

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

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 5487

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
5487 x2 = x * sx
y2 = y * sy
GOTO 6900

REM Lee/Adams Ellipse
5490 dp = -180 / k4: IF aa < 0 THEN dp = -dp

oo = .5 * ao
w = SQR(COS(dr# * aa)) / COS(.5 * dr# * aa)
GOSUB 7000

wa = (90 - w): IF aa < 0 THEN wa = -wa

v1 = 0: v2 = 90: v3 = 0
wo = oo
GOSUB 7300

IF aa >= 0 THEN 5491
ca = -ca
IF co < 0 THEN co = -180 - co ELSE co = 180 - co

5491 ca = -ca: co = -co
IF (ABS(aa) > 88) THEN 5493

im = LOG(TAN(dr# * .5 * (ABS(ca) + 90)))
IF aa < 0 THEN im = -im
re = dr# * co

mm = k1: GOSUB 7200

xx = .5 * dr# * k2 * x
yy = .5 * dr# * (k2 * y + dp)
GOTO 5495

REM im goes to infinity as you approach the poles, so this is needed
5493 qo = .5 * co
w = SQR(COS(dr# * ca)) / COS(.5 * dr# * ca)
GOSUB 7000
qa = (90 - w): IF aa < 0 THEN qa = -qa

v1 = 90: v2 = 90: v3 = 90
wa = qa: wo = qo
GOSUB 7300

im = LOG(TAN(dr# * .5 * (ABS(ca) + 90)))
IF ca < 0 THEN im = -im
re = dr# * co

mm = 1 - k1: GOSUB 7200

xx = .5 * dr# * -k2 * y
yy = .5 * dr# * (k2 * x + dp)

5495 REM now for sine part

r1 = EXP(yy): r2 = 1 / r1
y2 = -COS(xx) * (r1 - r2) * k3 * 0.5
x2 = -SIN(xx) * (r1 + r2) * k3 * 0.5
GOTO 6900

REM Ginzoid
5500 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 < 1.5 THEN 5505

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 5507

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

5507 x2 = xx * k1 + dx
y2 = sa * yy * k1
GOTO 6900

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

REM Laskowski Tri-Optimal
5520 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 6900

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

REM Imitation Canters' Minimum-Error Polyconic
5540 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 6900

REM Canters' Minimum-Error Polyconic
5550 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 6900

REM Canters' Minimum-Error Polyconic II
5560 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 6900

REM Adams-Cahill Dixon elliptic
5570 oc = INT((ao + 180) / 90)
zo = (ao - 90 * oc) - 45
oh = 1: IF aa < 0 THEN aa = -aa: oh = 2

REM if aa > 50 then tr = 1: ca = aa: co = zo: GOTO 5575
tr = 1: ca = aa: co = zo: GOTO 5575

wa = aa: wo = zo
tr = 2: IF zo < 0 THEN tr = 3: GOTO 5573

v1 = 45: v2 = -90: v3 = 45
GOSUB 7300
GOTO 5575

5573 v1 = -45: v2 = -90: v3 = -45
GOSUB 7300

5575 w = (COS(dr# * ca)) ^ 0.6666667 / COS(2 * dr# * ca / 3)
GOSUB 7000
zo = 2 * co / 3
za = w

IF zo >= 0 THEN wo = 30 - zo: si = 1 ELSE wo = 30 + zo: si = 2
p = 90 - za
li = 6 * wo
li = li - 360 * INT(li / 360)

m = 3: n = 1: pp = 1: nu = 1

tp = TAN(dr# * p / 2)
tm = tp ^ 6

u = tp * COS(dr# * wo)
v = tp * SIN(dr# * wo)

ll = wo
ta = tp

REM for i = 2 to 19
FOR i = 2 TO 43
    pp = pp * (nu / m)
    n = n + 6
    mz = pp / n
    ta = ta * tm

    ll = ll + li
    IF ll > 360 THEN ll = ll - 360

    u = u + mz * ta * COS(dr# * ll)
    v = v + mz * ta * SIN(dr# * ll)

NEXT i

u = u * k2: v = v * k2

x2 = u
y2 = v

GOTO 6900

REM Generalized Aitoff Equal-Area
5580 w = sh * SIN(dr# * aa): GOSUB 7000: 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 6900

REM Bromley-Mollweide
5590 pp = dr# * aa
k = pi# * .5 * SIN(pp)
IF aa = 0 THEN ps = 0: GOTO 5595
IF ABS(aa) = 90 THEN ps = pp: GOTO 5595
5592 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 5592
5595 x2 = vo * COS(ps) + dx
y2 = k2 * SIN(ps)
GOTO 6900

REM Modified Gall's
5600 IF ABS(aa) > 89.9 THEN pv% = 0: RETURN
y2 = TAN(dr# * .5 * aa) * k1
sz = aa * aa
x2 = ao * ((1.0 - (9.75E-6 * sz)) - (1.419753E-9 * sz * sz))
GOTO 6900

REM Equal Earth
5610 w = k1 * SIN(dr# * aa)
GOSUB 7000
sz = dr# * dr# * w * w
ca = COS(dr# * w)
x2 = ao * 4.020792 * ca / (4.020792 - 0.729954 * sz + sz * sz * sz * (0.018753 + 0.102492 * sz))
y2 = w * k2 * (1.340264 - 0.081106 * sz + sz * sz * sz * (0.000893 + 0.003796 * sz))
GOTO 6900

REM Simple Conic with Two Standard Parallels
5620 r = k2 + (s2 - aa)
a = dr# * k1 * ao
x2 = r * SIN(a)
y2 = k2 - r * COS(a)
GOTO 6900

REM Regional
5630 r = k2 + (s2 - aa)
IF aa > s2 THEN 5633
IF aa < s1 THEN 5635
a = dr# * k1 * ao
GOTO 5639
5633 a = dr# * k1 * dx + vo * COS(dr# * aa) / r
GOTO 5639
5635 a = dr# * k1 * dx + vo * COS(dr# * aa) / r
5639 x2 = r * SIN(a)
y2 = k2 - r * COS(a)
GOTO 6900

REM Stereographic
5640 rr = k1 * 90 * TAN(dr# * (90 - aa) / 2)
x2 = rr * SIN(dr# * ao)
y2 = -rr * COS(dr# * ao)
GOTO 6900

REM Airy Minimum-Error
5650 IF 90 - aa > rg THEN pv% = 0: RETURN
aq = dr# * 0.5 * (90 - aa)
sr = SIN(aq)
cr = COS(aq)
IF ABS(sr) < .01 THEN 5652
ck = cr / sr
rr = k2 * ((0 - (ck * LOG(cr))) + k1 / ck)
GOTO 5655
5652 tk = sr / cr
rr = k2 * tk * k1
5655 x2 = rr * SIN(dr# * ao)
y2 = -rr * COS(dr# * ao)
GOTO 6900

REM Dietrich-Kitada
5660 oa = ABS(ao)
IF oa < 2 THEN 5669
a = oa / 90
IF a < 1 THEN 5661
IF a > 1 THEN 5664

c = 0
r = 1
rr = k2 * 90
ag = 180

GOTO 5668

5661 xi = a
xo = 0: ap = 0
5662 r = .5 * (xi * xi + 1) / xi
w = 1 / r: GOSUB 7000

ca = pi# * r * r
pc = (w / 180) * ca
tr = r - xi
ak = 2 * (pc - tr)/pi#

IF ABS(ak - a) < .002 THEN 5663
xt = xi
xi = xi + (xi - xo) * (a - ak) / (ak - ap)
xo = xt: ap = ak
IF xi > 1 THEN xi = 1
GOTO 5662
5663 ag = 2 * w
c = k2 * 90 * (xi - r)
rr = k2 * 90 * r
GOTO 5668

5664 xi = a
IF xi > 1.653079 THEN xi = 1.653079
xo = 0: ap = 0
5665 r = .5 * (xi * xi + 1) / xi
w = 1 / r: GOSUB 7000
y = r + r - xi
ak = 2 * (((r - y)/pi#) + r * r * (1 - (w / 180)))
IF ABS(ak - a) < .002 THEN 5666
xt = xi
xi = xi + (xi - xo) * (a - ak) / (ak - ap)
xo = xt: ap = ak
IF xi > 1.653079 THEN xi = 1.653079
IF xi < 1 THEN xi = 1
GOTO 5665
5666 ag = 360 - 2 * w
c = k2 * 90 * (xi - r)
rr = k2 * 90 * r
5668 ng = dr# * aa * ag / 180
x2 = SGN(ao) * ((rr * COS(ng)) + c)
y2 = rr * SIN(ng)
GOTO 6900

REM parabolic case
5669 x2 = k1 * ao * COS(dr# * aa)
y2 = k2 * aa
GOTO 6900

REM Modified Regional
5670 IF aa < s1 THEN 5674
r = k2 + (s3 - aa)
REM IF aa > s4 THEN 5672
REM IF aa > s2 THEN 5676
REM temporary until parabolic is implemented
IF aa > s3 THEN 5672
REM Simple Conic
a = dr# * k1 * ao
GOTO 5639
REM Bonne's
5672 a = dr# * k1 * dx + vo * COS(dr# * aa) / r
GOTO 5639
REM Conic conformal
5674 r = k4 * TAN(.5 * dr# * (90 - aa)) ^ k3
a = dr# * k3 * ao
x2 = r * SIN(a)
y2 = 90 - r * COS(a)
GOTO 6900
REM Parabolic

5679 x2 = r * SIN(a)
y2 = k2 - r * COS(a)
GOTO 6900

REM Lagrange Conformal
5680 oo = k1 * ao
w = ((COS(dr# * aa))^k1) / COS(.5 * dr# * aa)
GOSUB 7000
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 6900

REM Expanding Parabolic
5690 x = ABS(ao/180): xs = 1 - x*x
k = k1 * (1-xs) + k2 * xs
y = ABS(aa/(90+k))
x2 = ao * (1 - y * y)
y2 = aa
GOTO 6900

REM Expanding Elliptical
5700 x = ABS(ao/180): xe = 1 - x*x: xm = 1 - xe * xe * xe
k = k3 * (1-xm) + k2 * xm
y = ABS(aa/(90+k))
IF ABS(ao) < 1.5 THEN 5703
kk = k1 * (1-xm) + k2 * xm
REM PRINT "Intersection: ";kk;"Curve height: ";k
yh = (90+kk)/(90+k)
p1 = SQR(1 - yh * yh)
p2 = SQR(1 - y * y)
x2 = ao * ( 1 - ((1 - p2)/(1 - p1)) )
GOTO 5707

REM Parabolic case
5703 x2 = ao * (1 - .5 * y * y)
5707 y2 = aa
GOTO 6900

REM Generalized Wagner Equal-Area
5710 w = sh * SIN(dr# * aa): GOSUB 7000: ap = w
di = SQR(1 + COS(dr# * ap) * COS(dr# * hs * ao))
y2 = 90 * ve * ex * SIN(dr# * ap) / di
x2 = 180 * he * COS(dr# * ap) * SIN(dr# * hs * ao) / di
GOTO 6900

REM Henri Bouthillier de Beaumont's Projection
5720 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 5727
IF qa = 0 THEN x = qo: y = 0: GOTO 5727
IF qo = 0 THEN x = 0: y = qa: GOTO 5727

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

REM calculate radii and centers of intersecting circles
a = qo: b = 90
c2 = (a * a + b * b) / (2 * a) - a
r2 = c2 + a
a = 112.5 * SIN(dr# * 1.4096655 * qa) - qa: b = 112.5 * COS(dr# * 1.4096655 * qa) + 67.5
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
5727 x2 = x * sx
y2 = y * sy
GOTO 6900

6900 GOSUB 8000
RETURN

REM mathematical subroutines

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

REM atan2
7020 IF ABS(wy) < ABS(wx) THEN 7030
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
7030 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

7100 a = 1: b = SQR(1 - m): c = SQR(m)
t1 = 1
7110 t2 = (rd# * w) / 90: IF t2 = INT(t2) THEN wa = w: GOTO 7120
wa = ATN(b * TAN(w) / a)
IF wa < 0 THEN wa = wa + pi#
wa = wa + pi# * INT(w / pi#)
7120 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 7110
w = w / (t1 * a)
RETURN

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

7200 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 7277
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 7100: x = w
7277 w = mu: m = 1 - mm: GOSUB 7100: y = w
x = x * sx
y = y * sy
RETURN

REM Coordinate conversion as a subroutine
7300 IF v2 = 0 THEN 7350
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 7310
ca = 90: IF qy < 0 THEN ca = -90
co = 0
REM Near the poles, computing longitude is meaningless
GOTO 7340

7310 IF ABS(qz) > ABS(qx) THEN 7320
co = rd# * ATN(qz / qx)
IF qx < 0 THEN co = co + 180
GOTO 7340

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

7340 co = co + v3
GOTO 7370

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

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

RETURN

REM Symmetric interruption input subroutine
7800 ip% = 1
INPUT "Number of gores? "; gn
IF gn < 1 THEN 7850
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

7850 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."
7860 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 7860
RETURN

REM Interruption input subroutine
7900 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."
7910 INPUT "Standard meridian? "; ns(ic%)
INPUT "Interruption? "; ni(ic%)
ic% = ic% + 1
IF ni(ic% - 1) < 180 THEN 7910
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."
7980 INPUT "Standard meridian? "; ss(ic%)
INPUT "Interruption? "; si(ic%)
ic% = ic% + 1
IF si(ic% - 1) < 180 THEN 7980
RETURN

REM line drawing and point plotting subroutines

8000 IF pp% = 1 THEN pv% = 1

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

IF pv% = 0 THEN pv% = 1: GOTO 8020
GOSUB 8100
8020 x1 = x2
y1 = y2
RETURN

8100 IF rt% = 0 THEN 8110

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 8600

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 8120

8110 x2% = hw% + INT(sf * x2 + hd + .5)
y2% = hh% + INT(sf * y2 + vd + .5)
IF pp% = 1 THEN 8600
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%

8120 xd% = ABS(x1% - x2%)
xs% = SGN(x2% - x1%)
yd% = ABS(y1% - y2%)
ys% = SGN(y2% - y1%)
IF yd% > xd% THEN 8500
IF xd% = 0 THEN 8900

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

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

8600 ii% = 1

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

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

REM PRINT "Point "; x%; y%

GOSUB 9000
ii% = ii% + 1
GOTO 8610
8690 RETURN

REM Plot a single point
8900 x% = x1%: y% = y1%: GOSUB 9000
RETURN

9000 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