' Moon algorithm from Meeus Astronomical Formulae for Calculators ch 30
DEFDBL A-Z
pr$ = "####.####"
pr2$ = "##.######"
CLS
pi = 4 * ATN(1)
tpi = 2 * pi
degs = 180 / pi
rads = pi / 180
DEF FNipart (x) = SGN(x) * INT(ABS(x))
'   the function below returns an angle in the range
'   0 to 360 degs
DEF FNr (x)
    B1 = x / 360
    a1 = 360 * (B1 - FNipart(B1))
    IF a1 < 0 THEN a1 = 360 + a1
    FNr = a1
END DEF
DEF FNatn2 (y, x)
    a = ATN(y / x)
    IF x < 0 THEN a = a + pi
    IF (y < 0) AND (x > 0) THEN a = a + tpi
    FNatn2 = a
END DEF
DEF FNasin (x)
    c = SQR(1 - x * x)
    T = x / c
    FNasin = ATN(T)
END DEF
DEF FNday (y, M, D, h) = 367 * y - 7 * (y + (M + 9) \ 12) \ 4 + 275 * M \ 9 + D - 694006.5 + h / 24
CLS
OPEN "O", #1, "moontab.txt"
' one saros
FOR days = 29590 TO 36525
PRINT ".";
h = h + mins / 60
T = days / 36525
'  Mean arguments
 Om = FNr(259.183275# - 1934.142# * T + .002078 * T * T + .0000022# * T * T * T) * rads
 Lm = FNr(270.434164# + 481267.8831# * T - .001133 * T * T + .0000019# * T * T * T)
 Lm = Lm + .000233 * SIN((51.2 + 20.2 * T) * rads)
 Lm = Lm + .003964 * SIN((346.56 + 132.87 * T - .0091731# * T * T) * rads)
 Lm = Lm + .001964 * SIN(Om)
 Lm = Lm * rads
 M = FNr(358.475833# + 35999.0498# * T - .00015 * T * T - .0000033# * T * T * T)
 M = M - .001778 * SIN((51.2 + 20.2 * T) * rads)
 M = M * rads
 Mm = FNr(296.104608# + 477198.8491# * T + .009192 * T * T + .0000144# * T * T * T)
 Mm = Mm + .000817 * SIN((51.2 + 20.2 * T) * rads)
 Mm = Mm + .003964 * SIN((346.56 + 132.87 * T - .0091731# * T * T) * rads)
 Mm = Mm + .002541 * SIN(Om)
 Mm = Mm * rads
 D = FNr(350.737486# + 445267.1142# * T - .001436 * T * T + .0000019# * T * T * T)
 D = D + .002011 * SIN((51.2 + 20.2 * T) * rads)
 D = D + .003964 * SIN((346.56 + 132.87 * T - .0091731# * T * T) * rads)
 D = D + .001964 * SIN(Om)
 D = D * rads
 F = FNr(11.250889# + 483202.0251# * T - .003211 * T * T - .0000003# * T * T * T)
 F = F + .003964 * SIN((346.56 + 132.87 * T - .0091731# * T * T) * rads)
 F = F - .024691 * SIN(Om)
 F = F - .004328 * SIN(Om + (275.05 - 2.3 * T) * rads)
 F = F * rads
 e = 1 - .002495 * T - .00000752# * T * T
 e2 = e * e
 ' Correction for geocentric ecliptic longitude for Moon
 dl = 6.28875 * SIN(Mm)
 dl = dl + 1.274018 * SIN(2 * D - Mm)
 dl = dl + .658309 * SIN(2 * D)
 dl = dl + .213616 * SIN(2 * Mm)
 dl = dl - e * .185596 * SIN(M)
 dl = dl - .114336 * SIN(2 * F)
 dl = dl + .058793 * SIN(2 * D - 2 * Mm)
 dl = dl + e * .057212 * SIN(2 * D - M - Mm)
 dl = dl + .05332 * SIN(2 * D + Mm)
 dl = dl + e * .045874 * SIN(2 * D - M)
 dl = dl + e * .041024 * SIN(Mm - M)
 dl = dl - .034718 * SIN(D)
 dl = dl - e * .030465 * SIN(M + Mm)
 dl = dl + .015326 * SIN(2 * D - 2 * F)
 dl = dl - .012528 * SIN(2 * F + Mm)
 dl = dl - .01098 * SIN(2 * F - Mm)
 dl = dl + .010674 * SIN(4 * D - Mm)
 dl = dl + .010034 * SIN(3 * Mm)
 dl = dl + .008548 * SIN(4 * D - 2 * Mm)
 dl = dl - e * .00791 * SIN(M - Mm + 2 * D)
 dl = dl - e * .006783 * SIN(2 * D + M)
 dl = dl + .005162 * SIN(Mm - D)
 dl = dl + e * .005 * SIN(M + D)
 dl = dl + e * .004049 * SIN(Mm - M + 2 * D)
 dl = dl + .003996 * SIN(2 * Mm + 2 * D)
 dl = dl + .003862 * SIN(4 * D)
 dl = dl + .003665 * SIN(2 * D - 3 * Mm)
 dl = dl + e * .002695 * SIN(2 * Mm - M)
 dl = dl + .002602 * SIN(Mm - 2 * F - 2 * D)
 dl = dl + e * .002396 * SIN(2 * D - M - 2 * Mm)
 dl = dl - .002349 * SIN(Mm + D)
 dl = dl + e2 * .002249 * SIN(2 * D - 2 * M)
 dl = dl - e * .002125 * SIN(2 * Mm + M)
 dl = dl - e2 * .002079 * SIN(2 * M)
 dl = dl + e2 * .002059 * SIN(2 * D - Mm - 2 * M)
 dl = dl - .001773 * SIN(Mm + 2 * D - 2 * F)
 dl = dl - .001595 * SIN(2 * F + 2 * D)
 dl = dl + e * .00122 * SIN(4 * D - M - Mm)
 dl = dl - .00111 * SIN(2 * Mm + 2 * F)
 dl = dl + .000892 * SIN(Mm - 3 * D)
 dl = dl - e * .000811 * SIN(M + Mm + 2 * D)
 dl = dl + e * .000761 * SIN(4 * D - M - 2 * Mm)
 dl = dl + e2 * .000717 * SIN(Mm - 2 * M)
 dl = dl + e2 * .000704 * SIN(Mm - 2 * M - 2 * D)
 dl = dl + e * .000693 * SIN(M - 2 * Mm + 2 * D)
 dl = dl + e * .000598 * SIN(2 * D - M - 2 * F)
 dl = dl + .00055 * SIN(Mm + 4 * D)
 dl = dl + .000538 * SIN(4 * Mm)
 dl = dl + e * .000521 * SIN(4 * D - M)
 dl = dl + .000486 * SIN(2 * Mm - D)
 lambda = Lm + dl * rads
 '  B series for Moon - very nealy the geocentric ecliptic latitude
 B = 5.128189 * SIN(F)
 B = B + .280606 * SIN(Mm + F)
 B = B + .277693 * SIN(Mm - F)
 B = B + .173238 * SIN(2 * D - F)
 B = B + .055413 * SIN(2 * D + F - Mm)
 B = B + .046272 * SIN(2 * D - F - Mm)
 B = B + .032573 * SIN(2 * D + F)
 B = B + .017198 * SIN(2 * Mm + F)
 B = B + .009267 * SIN(2 * D + Mm - F)
 B = B + .008823 * SIN(2 * Mm - F)
 B = B + e * .008247 * SIN(2 * D - M - F)
 B = B + .004323 * SIN(2 * D - F - 2 * Mm)
 B = B + .0042 * SIN(2 * D + F + Mm)
 B = B + e * .003372 * SIN(F - M - 2 * D)
 B = B + e * .002472 * SIN(2 * D + F - M - Mm)
 B = B + e * .002222 * SIN(2 * D + F - M)
 B = B + e * .002072 * SIN(2 * D - F - M - Mm)
 B = B + e * .001877 * SIN(F - M + Mm)
 B = B + .001828 * SIN(4 * D - F - Mm)
 B = B - e * .001803 * SIN(F + M)
 B = B - .00175 * SIN(3 * F)
 B = B + e * .00157 * SIN(Mm - M - F)
 B = B - .001487 * SIN(F + D)
 B = B - e * .001481 * SIN(F + M + Mm)
 B = B + e * .001417 * SIN(F - M - Mm)
 B = B + e * .00135 * SIN(F - M)
 B = B + .00133 * SIN(F - D)
 B = B + .001106 * SIN(F + 3 * Mm)
 B = B + .00102 * SIN(4 * D - F)
 B = B + .000833 * SIN(F + 4 * D - Mm)
 B = B + .000781 * SIN(Mm - 3 * F)
 B = B + .00067 * SIN(F + 4 * D - 2 * Mm)
 B = B + .000606 * SIN(2 * D - 3 * F)
 B = B + .000597 * SIN(2 * D + 2 * Mm - F)
 B = B + e * .000492 * SIN(2 * D + Mm - M - F)
 B = B + .00045 * SIN(2 * Mm - F - 2 * D)
 B = B + .000439 * SIN(3 * Mm - F)
 B = B + .000423 * SIN(F + 2 * D + 2 * Mm)
 B = B + .000422 * SIN(2 * D - F - 3 * Mm)
 B = B - e * .000367 * SIN(M + F + 2 * D - Mm)
 B = B - e * .000353 * SIN(M + F + 2 * D)
 B = B + .000331 * SIN(F + 4 * D)
 B = B + e * .000317 * SIN(2 * D + F - M + Mm)
 B = B + e2 * .000306 * SIN(2 * D - 2 * M - F)
 B = B - .000283 * SIN(Mm + 3 * F)
 B = B * rads
 W1 = .0004664# * COS(Om)
 W2 = .0000754# * COS(Om + (275.05 - 2.3 * T) * rads)
 beta = B * (1 - W1 - W2)
 '  Parallax of moon
 Pm = .950724#
 Pm = Pm + .051818# * COS(Mm)
 Pm = Pm + .009531# * COS(2 * D - Mm)
 Pm = Pm + 7.842999999999999D-03 * COS(2 * D)
 Pm = Pm + .002824# * COS(2 * Mm)
 Pm = Pm + .000857# * COS(2 * D + Mm)
 Pm = Pm + e * 5.330000000000001D-04 * COS(2 * D - M)
 Pm = Pm + e * .000401# * COS(2 * D - M - Mm)
 Pm = Pm + e * .00032# * COS(Mm - M)
 Pm = Pm - .000271# * COS(D)
 Pm = Pm - e * .000264# * COS(M + Mm)
 Pm = Pm - .000198# * COS(2 * F - Mm)
 Pm = Pm + .000173# * COS(3 * Mm)
 Pm = Pm + .000167# * COS(4 * D - Mm)
 Pm = Pm - e * .000111# * COS(M)
 Pm = Pm + .000103# * COS(4 * D - 2 * Mm)
 Pm = Pm - .000084# * COS(2 * Mm - 2 * D)
 Pm = Pm - e * .000083# * COS(2 * D + M)
 Pm = Pm + .000079# * COS(2 * D + 2 * Mm)
 Pm = Pm + .000072# * COS(4 * D)
 Pm = Pm + e * .000064# * COS(2 * D - M + Mm)
 Pm = Pm - e * .000063# * COS(2 * D + M - Mm)
 Pm = Pm + e * .000041# * COS(M + D)
 Pm = Pm + e * .000035# * COS(2 * Mm - M)
 Pm = Pm - .000033# * COS(3 * Mm - 2 * D)
 Pm = Pm - .00003# * COS(Mm + D)
 Pm = Pm - .000029# * COS(2 * F - 2 * D)
 Pm = Pm - e * .000029# * COS(2 * Mm + M)
 Pm = Pm + e2 * .000026# * COS(2 * D - 2 * M)
 Pm = Pm - .000023# * COS(2 * F - 2 * D + Mm)
 Pm = Pm + e * .000019# * COS(4 * D - M - Mm)
 '  nutation in longitude from Duffett-Smith 'Astronomy with
 '  your calculator' 3rd ed section 35, p60
 L = FNr(279.6967 + 36000.7689# * T + .000303 * T * T) * rads
 dphi = (-17.2 * SIN(Om) - 1.3 * SIN(2 * L)) * rads / 3600
 depsilon = (9.2 * COS(Om) + .5 * COS(2 * L)) * rads / 3600
 lambda = lambda + dphi
 PRINT #1, days + 2415020; " , "; lambda * degs; " , "; degs * beta; " , "; Pm
 NEXT days
 CLOSE #1
 END