include Julian

-- graph phase of moon for a month

t = JD(1996,2,X)

plot (1-cos(l2-lambda_sun))/2
Xmin = 0; Xmax = 31
Ymin = 0; Ymax = 1; Ydiv = .25
---------------------------------------------------
-- from "Practical astronomy with your calculator"
-- by Peter Duffet-Smith
-- 3rd edition 1988

D = t-2447891.5+55/86400  -- TDT = UT+55 in 1986
l = 13.1763966*D+l0
M_moon = l-0.1114041*D-P0
Ev = 1.2739*sin(2*C-M_moon)
Ae = 0.1858*sin(M_sun)
A3 = 0.37*sin(M_sun)
C = l-lambda_sun
M_moon1 = M_moon+Ev-Ae-A3
Ec = 6.2886*sin(M_moon1)
A4 = 0.214*sin(2*M_moon1)
l1 = l+Ev+Ec-Ae+A4
V = 0.6583*sin(2*(l1-lambda_sun))
l2 = l1+V

l0 = 318.351648
P0 = 36.340410
N0 = 318.510107

---------------------------------------------------
M_sun = 360/365.242191*D+epsilong-omegabarg
lambda_sun = v+omegabarg
v = M_sun+360/pi*e_sun*sin(M_sun)

epsilong = 279.403303
omegabarg = 282.768422
e_sun = 0.016713