;+
;FUNCTION:
;   moon_ecl_pos
;PURPOSE:
;   Calculates the approximate Geocentric ecliptic coordinates of the Moon
;   (in deg) for a given time.  The method was copied from the Astronomical
;   Almanac for 1988 (pg. E46) and tested on tabulated values in that volume.
;
;   Returns and anonymous structure with the following elements:
;
;     lon  : ecliptic longitude  [degrees]
;     lat  : ecliptic latitude   [degrees]
;     dist : Earth-Moon distance [Earth radii]
;
;   Precision:  0.3 deg in longitude
;               0.2 deg in latitude
;               0.2 Earth radii in distance
;
;USAGE:
;   moonpos = moon_ecl_pos(time)
;INPUTS:
;       time:   Any valid time specification for 'time_double'.
;
;KEYWORDS:
;       RAD:    Return the angles in radians instead of degrees.
;
;CREATED BY:	David L. Mitchell  03-25-98
;FILE:  moon_ecl_pos.pro
;VERSION:  1.1
;LAST MODIFICATION:  03-25-98
;-
function moon_ecl_pos, time, rad=rad, equatorial=equatorial,pos=pos

  rtod = 180D/!dpi

  time = time_double(time)
  J2000 = time_double('2000-01-01/12:00:00')

  dt = ((time - J2000)/86400D)/36525D

  lon = 218.32D + 481267.883D*dt                        $
        + 6.29D*sin((134.9D + 477198.85D*dt)/rtod)      $
        - 1.27D*sin((259.2D - 413335.38D*dt)/rtod)      $
        + 0.66D*sin((235.7D + 890534.23D*dt)/rtod)      $
        + 0.21D*sin((269.9D + 954397.70D*dt)/rtod)      $
        - 0.19D*sin((357.5D +  35999.05D*dt)/rtod)      $
        - 0.11D*sin((186.6D + 966404.05D*dt)/rtod)

  nwraps = floor(lon/360D)
  lon = lon - (360D*double(nwraps))

  lat =   5.13D*sin(( 93.3D + 483202.03D*dt)/rtod)      $
        + 0.28D*sin((228.2D + 960400.87D*dt)/rtod)      $
        - 0.28D*sin((318.3D +   6003.18D*dt)/rtod)      $
        - 0.17D*sin((217.6D - 407332.20D*dt)/rtod)

  plx =   0.9508D                                       $
        + 0.0518D*cos((134.9D + 477198.85D*dt)/rtod)    $
        + 0.0095D*cos((259.2D - 413335.38D*dt)/rtod)    $
        + 0.0078D*cos((235.7D + 890534.23D*dt)/rtod)    $
        + 0.0028D*cos((269.9D + 954397.70D*dt)/rtod)

  dist = 1D/sin(plx/rtod)

    
    x = [dist * cos(lat/rtod) * cos(lon/rtod) ]
    y = [dist * cos(lat/rtod) * sin(lon/rtod) ]
    z = [dist * sin(lat/rtod)]
    pos = [[x],[y],[z]]
    
    
  if  keyword_set(equatorial) then begin
    tilt = 23.45/rtod
    ct = cos(tilt)
    st = sin(tilt)
    rotmat=[[1,0,0],[0,ct,-st],[0,st,ct]]
    pos = pos # rotmat
    x = pos[*,0] &  y = pos[*,1]  &  z = pos[*,2]
    cart_to_sphere,x,y,z,dummy,lat,lon
    if n_elements(dist) le 1 then begin & lat=lat[0] & lon=lon[0] & endif
  endif

  if keyword_set(rad) then begin
    lon = lon/rtod
    lat = lat/rtod
  endif
  
  str = replicate({lon : 0.d, lat : 0.d, dist : 0.d},n_elements(lon))
  str.lon = lon
  str.lat = lat
  str.dist = dist
  return, str

end