pro aa_sun,jd,sun
; Use low precision formulas for sun position from the
; Astronomical Almanac (1995) [C24] to calculate various
; solar quantities. Takes fractional julian day (jd) as input
; (jd may be an array). Good to .01 degree between 1950-2050
n=jd-2451544.5
L=(280.466+.9856474*n) mod 360. ;mean solar longitude (deg)
g=(357.528+.9856003*n) mod 360. ;mean anomaly (deg)
;place in range 0 to 360.
il=where(L lt 0)
if il[0] ge 0 then L[il]=L[il]+360.
ig=where(g lt 0)
if ig[0] ge 0 then g[ig]=g[ig]+360.
lam=L+1.915*sin(g*!dtor)+.02*sin(2.*g*!dtor) ;ecliptic longitude
eps=23.44-.0000004*n ;obliquity of the ecliptic (degrees)
; right asc. (degrees)
ra=atan(cos(eps*!dtor)*tan(lam*!dtor))*!radeg
ir=where(L ge 90 and L le 270)
if ir[0] ge 0 then ra[ir]=ra[ir]+180.
ra=ra mod 360.
ir2=where(ra lt 0)
if ir2[0] ge 0 then ra[ir2]=ra[ir2]+360.
dec=asin(sin(eps*!dtor)*sin(lam*!dtor))*!radeg ; dec. (degrees)
dist=1.00014-.01671*cos(g*!dtor)-.00014*cos(2.*g*!dtor) ; distance from earth (au)
x=dist*cos(lam*!dtor) ; rectangular position (au)
y=dist*cos(eps*!dtor)*sin(lam*!dtor)
z=dist*sin(eps*!dtor)*sin(lam*!dtor)
sun={jd:jd,n:n,L:L,g:g,lam:lam,eps:eps,ra:ra,dec:dec,dist:dist, $
x:x,y:y,z:z}
return
end