Viewing contents of file '../idllib/contrib/buie/ephem.pro'
;+
; NAME:
;     ephem
; PURPOSE: (one line)
;     Ephemeris generator for solar system objects.
; DESCRIPTION:
;     This program pipes the input data to an external FORTRAN program
;     written by Larry Wasserman.  The external program looks up the orbital
;     elements of each object and computes the positions requested by the
;     output code.  The output code options are:
;
;        Code  Meaning
;        -----------------------------------------------------------------------
;         0    Topocentric RA, Dec, of date of object.
;         1    Topocentric RA, Dec,  B1950  of object.
;         2    Topocentric RA, Dec,  J2000  of object.
;         3    Topocentric x, y, z, of date of object.
;         4    Topocentric x, y, z,  B1950  of object.
;         5    Topocentric x, y, z,  J2000  of object.
;         6    Topocentric x, y, z, of date of Sun.
;         7    Topocentric x, y, z,  B1950  of Sun.
;         8    Topocentric x, y, z,  J2000  of Sun.
;         9    B1950 x, y, z, of object with respect to solar system barycenter.
;        10    J2000 x, y, z, of object with respect to solar system barycenter.
;        11    J2000 osculating orbital elements for the object
;                 M, arg peri, node, i, e, a, q, Q, Period, epoch,
;                 VRA, VDec, Epherr, JDlast, arc, nobs
;                 (16 values)
;                 VRA, VDec, line of variations for +1 deg change in M (radians)
;                 JDlast - date of last astrometric measurement (nearest day)
;                 arc - Arc length of astrometric observations (days)
;                 nobs - Number of observations used to generate orbit.
;
;        20    Code 0 and Code 3 and Code 6 information.
;        21    Code 1 and Code 4 and Code 7 information.
;        22    Code 2 and Code 5 and Code 8 information.
;        23    Code 2, 5, 8, and 11 information.
;
;     The object code depends on the type of object requested.
;     All asteroids are preceeded by A, comets by C, planets and satellites
;     by P, and topocentric locations by T.  For asteroids, either the number
;     of the asteroid or its provisional designation will work.  Comets
;     are recognized by either provisional or final designations, not names.
;     Planets and satellites use the NAIF coding scheme for objects.  For
;     example, the Pluto-Charon barycenter is 9, Charon is 901 and the
;     center of Pluto is 999.  Topocentric locations use the same code as
;     for the observatory codes, for example, Anderson Mesa would be T688.
;     None of the names allow spaces.  If the object isn't recognized or
;     supported by the external databases, RA and Dec values (if requested)
;     are returned as -99 and the x,y,z values are returned as 0.  Elements
;     will return as zero for invalid input.
;
; CATEGORY:
;  Astronomy
; CALLING SEQUENCE:
;     ephem,jd,obs,code,object,ephemeris
; INPUTS:
;     jd     - Julian Date (UT) of position to calculate (scalar or vector).
;     obs    - Observatory code, Marsden's IAUC codes (scalar).
;     code   - Output format code (see description) (scalar).
;     object - Object code (see description) (string or string array).
;
; OPTIONAL INPUT PARAMETERS:
; KEYWORD PARAMETERS:
; OUTPUTS:
;     epehemeris - 2-D double precision array of values requested by code.
;
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
; MODIFICATION HISTORY:
;     3/24/93 - IDL interface written by Marc W. Buie, Lowell Observatory.
;-
pro ephem,jd,obs,code,object,ephemeris, DEBUG=debug

on_ioerror,bail_out

;Validate number of parameters.
if n_params() ne 5 then begin
   print,'Usage:  ephem,jd,obs,code,object,ephemeris'
   return
end

; Validate input arguments.
if badpar(jd,5,[0,1],caller='EPHEM: (jd) ' ,npts=npts) then return
if badpar(obs,2,[0,1],caller='EPHEM: (obs) ',npts=nobs) then return
if badpar(code,2,0,caller='EPHEM: (code) ') then return
if badpar(object,7,[0,1], $
            caller='EPHEM: (object) ',RANK=objrank,NPTS=nobj) then return

; Ensure match between object list and time.
if nobj ne 1 and npts ne nobj then begin
   print,'EPHEM: Error.  The length of object must be 1 or match the length of JD.'
   return
endif

; Ensure match between observatory list and time.
if nobs ne 1 and npts ne nobs then begin
   print,'EPHEM: Error.  The length of obs must be 1 or match the length of JD.'
   return
endif

; Validate output code requested.
if code ge 50 then chk_code=code-50 else chk_code=code
case chk_code of 
   0 : nvals = 2
   1 : nvals = 2
   2 : nvals = 2
   3 : nvals = 3
   4 : nvals = 3
   5 : nvals = 3
   6 : nvals = 3
   7 : nvals = 3
   8 : nvals = 3
   9 : nvals = 3
  10 : nvals = 3
  11 : nvals = 16
  20 : nvals = 8
  21 : nvals = 8
  22 : nvals = 8
  23 : nvals = 24
  else : nvals = 0
endcase

if nvals eq 0 then begin
   print,'EPHEM: Error.  ',code,' is an illegal output code value.'
   return
endif

; Start up Larry's ephemeris generator (this assumes it's in the path).
spawn,'geteph',unit=pipe

; Create output and temporary arrays.
ephemeris=dblarr(nvals,npts)
tmp = dblarr(nvals)

;kludge for schwartz account
;printf,pipe,jd(0),obs[0],code,object(0),format='(f13.5,1x,i3,1x,i2,1x,a)'
aaaa=' '
;readf,pipe,aaaa

; Write all the values to the generator.
for i=0L,npts-1 do begin

   if nobj eq 1 then obj=object else obj=object[i]
   if nobs eq 1 then obse=obs else obse=obs[i]

   printf,pipe,jd[i],obse,code,obj,format='(f13.5,1x,i3,1x,i2,1x,a)'
   flush,pipe
   readf,pipe,aaaa

   if keyword_set(debug) then begin
      print,jd[i],obse,code,obj,format='(f13.5,1x,i3,1x,i2,1x,a)'
      print,aaaa
   endif

   reads,aaaa,tmp
   ephemeris[*,i] = tmp

endfor

; Scan the data to see if valid ephemerides were returned for all.
badflags = 0
if chk_code eq 0 or chk_code eq 1 or chk_code eq 2 or $
   chk_code eq 20 or chk_code eq 21 or chk_code eq 22 then $
   z = where(ephemeris[0,*] lt -90.0, badflags) $
else $
   z = where(ephemeris[0,*] eq 0.0 and $
             ephemeris[1,*] eq 0.0 and $
             ephemeris[2,*] eq 0.0, badflags)

if badflags ne 0 then begin
   bad = object[z]
   bad = bad[uniq(bad,sort(bad))]
   if n_elements(bad) eq 1 then begin
      print,'EPHEM:  Warning!  Object code ',bad[0],' did not return a valid ephemeris.'
   endif else begin
      print,'EPHEM:  Warning!  No ephemeris was returned for the following object codes.'
      for i=0L,n_elements(bad)-1 do print,'   ',bad[i]
   endelse
endif

; Close the pipe.
close,pipe
free_lun,pipe
return

; This for unexpected output from geteph
bail_out:
print,'EPHEM:  Illegal information returned by geteph:'
print,aaaa
free_lun,pipe
return

end