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



Page Last Modified: February 24, 2014