pro planet_coords, date, ra, dec, planet=planet, jd = jd, jpl = jpl
;+
; NAME:
;    PLANET_COORDS
; PURPOSE:  
;    Find low or high precision RA and DEC for the planets given a date
;
; EXPLANATION:
;    For low precision this routine uses HELIO to get the heliocentric ecliptic
;    coordinates of the planets at the given date, then converts these to 
;    geocentric ecliptic coordinates ala "Astronomical Algorithms" by Jean 
;    Meeus (1991, p 209). These are then converted to RA and Dec using EULER.
;    The accuracy between the years 1800 and 2050 is better than 1 arcminute 
;    for  the terrestial planets, but reaches 10 arcminutes for Saturn.    
;    Before 1850 or after 2050 the accuracy can get much worse.   
;
;    For high precision use the /JPL option ito use the full JPL ephemeris.
; CALLING SEQUENCE:
;    PLANET_COORDS, DATE, RA, DEC, [ PLANET = , /JD, /JPL]
;
; INPUTS:
;       DATE - If /JD is not set, then date is a 3-6 element vector containing
;              year,month (1-12), day, and optionally hour, minute, & second.
;              If /JD is set then DATE is a Julian date.   An advantage of the
;              /JD option is that it allows the use of vector dates.
; OUTPUTS:
;       RA - right ascension of planet(s), J2000 degrees, double precision
;       DEC - declination of   planet(s), J2000 degrees, double precision
;
; OPTIONAL INPUT KEYWORD:
;       PLANET - scalar string giving name of a planet, e.g. 'venus'. Default 
;               is to compute coords for all of them (except Earth).
;       /JD - If set, then the date parameter should be supplied as Julian date
;       JPL - if /JPL set, then PLANET_COORDS will call the procedure 
;             JPLEPHINTERP to compute positions using the full JPL ephemeris.   
;             The JPL ephemeris FITS file JPLEPH.405 must exist in either the 
;             current directory, or in the directory specified by the 
;             environment variable ASTRO_DATA.   Alternatively, the JPL keyword
;             can be set to the full path and name of the ephemeris file.
;             A copy of the JPL ephemeris FITS file JPLEPH.405 is available in
;                 http://idlastro.gsfc.nasa.gov/ftp/data/         
; EXAMPLES:
;    (1)  Find the RA, Dec of Venus on 2018 Dec 20
;          IDL> planet_coords, [2018,12,20], ra,dec    ;Compute for all planets
;          IDL> print,adstring(ra[1],dec[1],1)         ;Venus is second planet
;     ====> RA = 14 42 41.17  Dec = -12 30 53.5
;    This position is 61" from the full JPL ephemeris position of
;          RA = 14 42 45.32       -12 30 50.6
;
;    (2) Return the current RA and Dec of all 8 planets using JPL ephemeris
;          IDL> get_juldate, jd                 ;Get current Julian Date
;          IDL> planet_coords,jd,ra,dec,/jd,/jpl  ;Find positions of all planets
;          IDL> forprint,adstring(ra,dec,0)     ;Display positions   
;
;    (3) Plot the declination of Mars for every day in the year 2001
;          IDL> jdcnv,2001,1,1,0,jd      ;Get Julian date of midnight on Jan 1 
;               Now get Mars RA,Dec for 365 consecutive days
;          IDL> planet_coords,jd+indgen(365),ra,dec,/jd, planet = 'mars'     
;          IDL> plot,indgen(365)+1,dec
; NOTES:
;          HELIO is based on the two-body problem and neglects interactions 
;           between the planets.   This is why the worst results are for
;           Saturn.   Use the /JPL option or the online ephemeris generator 
;           http://ssd.jpl.nasa.gov/horizons.cgi for more accuracy. 
;
;           The procedure returns astrometric coordinates, i.e. no correction
;           for aberration.   A correction for light travel time is applied
;           when /JPL is set, but not for the default low-precision calculation.
; PROCEDURES USED:
;        JULDATE 
;        EULER, HELIO  - if /JPL is not set
;        JPLEPHREAD, JPLEPHINTERP - if /JPL is set
; REVISION HISTORY:
;        Written P.Plait & W. Landsman     August 2000
;        Fixed Julian date conversion   W. Landsman August 2000
;        Added /JPL keyword  W. Landsman   July 2001
;        Allow vector Julian dates with JPL ephemeris W. Landsman December 2002
;-   
; On_error,2
 if N_params() LT 1 then begin
     print,'Syntax - PLANET_COORDS, date, ra,dec, [PLANET =, /JD , JPL= ]'
     print,'      date - either 3-6 element date or Julian date (if /JD is set)'
     print,'      ra,dec - output ra and dec in degrees'
     print,'      PLANET - name of planet (optional)'
     return
 endif

 radeg = 180.0d/!DPI
 c = 2.99792458d5

;convert input date to real JD
  
  if keyword_set(jd) then begin 
       jj = date
       if N_elements(jj) GT 0 then if N_elements(planet) GT 1 then $
       message,'ERROR - A planet name must be supplied for vector dates'
  endif else begin 
           juldate,date,jj
            jj = jj + 2400000.0d
  endelse 

;make output arrays to include each planet
; note that we need Earth to convert from heliocentric
; ecliptic coordinates to geocentric and then to RA and DEC

   if keyword_set(planet) then begin
       planetlist = ['MERCURY','VENUS','MARS', $
        'JUPITER','SATURN','URANUS','NEPTUNE','PLUTO']
        index = 1+ where(planetlist eq strupcase(strtrim(planet,2)), Nfound) 
        if index[0] GE 3 then index = index + 1
       if Nfound EQ 0 then message,'Unrecognized planet of ' + planet
   endif else index = [1,2,4,5,6,7,8,9]

   if keyword_set(JPL) then begin
       if size(jpl,/TNAME) EQ 'STRING' then jplfile = jpl else $
            jplfile = find_with_def('JPLEPH.405','ASTRO_DATA')

        if jplfile EQ '' then message,'ERROR - Cannot find JPL ephemeris file' 
;Read ephemeris FITS file
        JPLEPHREAD,jplfile, pinfo, pdata, [long(min(jj)-1), long(max(jj)+1)]
        np = N_elements(index)
        njd = n_elements(jj)
        ra = dblarr(njd,np)   & dec = dblarr(njd,np)

        for i=0, Np-1 do begin
        JPLEPHINTERP, pinfo, pdata, jj, x,y,z, $
                       objectname=index[i],center='EARTH'
;  Compute distance to planet(s) and adjust Julian date for light travel time
; and recompute planet positions
       dis = sqrt(x^2 + y^2 + z^2)
       jj1 = jj - dis/c/86400.0d

; Compute position of Earth at current time, but position of planet at time
; light started traveling 
       JPLEPHINTERP, pinfo, pdata, jj, xe,ye,ze, /EARTH
       JPLEPHINTERP, pinfo, pdata, jj1, x,y,z, objectname=index[i]
       x = x-xe & y = y-ye & z = z-ze
         ra[0,i] = atan(y,x) * radeg
       g = where(ra LT 0, Ng)
       if Ng GT 0 then ra[g] = ra[g] + 360.0d
       dec[0,i]   = atan(z,sqrt(x*x + y*y)) * radeg
       endfor
       ra = reform(ra) & dec = reform(dec)
       return
   endif

   helio,jj,index,rad,lon,lat,/radian

; extract Earth-Moon barycenter info

   helio,jj,3,rade,lone,late,/radian

;get rectangular coords of planets

   x = rad * cos(lat) * cos(lon) - rade * cos(late) * cos(lone)
   y = rad * cos(lat) * sin(lon) - rade * cos(late) * sin(lone)
   z = rad * sin(lat)            - rade * sin(late)

;get geocentric longitude lambda and geo latitude, beta

   lambda = atan(y,x) * radeg
   beta   = atan(z,sqrt(x*x + y*y)) * radeg

;convert to Ra and Dec

   euler, lambda, beta, ra, dec, 4

   return
   end