;+
; NAME:
;     GET_SOLEPH
; PURPOSE:
;     Return a three-day consecutive solar ephemeris, given the
;     Julian date of the first day.
; CATEGORY:
;     OVRO APC EPHEMERIS
; CALLING SEQUENCE:
;     solar = get_soleph(jd)
; INPUTS:
;     jd      Julian date of the first day for which the solar
;               ephemeris data are wanted.
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
; ROUTINES CALLED:
;     rd_soleph (internal routine), pandb, get_sol_struct
; OUTPUTS:
;     solar   A solar ephemeris structure of the form
;                 { RA:    LONARR(3),  ; Right ascension for each day [0.1 mdeg]
;                   DEC:   LONARR(3),  ; Declination for each day [0.1 mdeg]
;                   RAD:   INTARR(3),  ; Radius for each day      [0.1 mdeg]
;                   PA:    LONARR(3),  ; P-Angle of solar disk    [0.1 mdeg]
;                   B0:    LONARR(3),  ; B0-Angle of solar disk   [0.1 mdeg]
;                   HP:    INTARR(3),  ; Horizontal parallax      [0.1 mdeg]
;                   REFT:  INTARR(3),  ; Reference time (always 0)[sec]
;                   DOY:   INTARR(3)}  ; Reference day of year
; COMMENTS:
;   Requires the file !DEFAULTS.EPHEM+'SUN1994-2020.DAT, which is the same file from
;   which the CPC (FASAR) file SOLAREPH.DAT was generated to determine the
;   solar coordinates.
; SIDE EFFECTS:
; RESTRICTIONS:
;     Negative declination is determined by dash in column 41!  Beware if
;     structure changes in future solar ephemeris files.
; MODIFICATION HISTORY:
;     Written 12-Jan-1999 by Dale E. Gary
;     30-Jan-1999  DG
;       Slight change to use GET_SOL_STRUCT() routine
;     11-Jan-2000  DG
;       Eliminated hardwired directory locations.
;     25-Jul-2000  DG
;       Discovered a bug in the solar declination near the equinoxes!
;       When, e.g., declination -00 54 13 was read, -00 was read as zero,
;       not negative zero, so sign was not transferred properly.  The routine
;       now keys on the dash representing the minus sign, so should be
;       correct, but if the structure of future files change, this will
;       have to be changed also.  Note that all files within one hour
;       declination of 0 (near the equinoxes) have incorrect ephemeris
;       (about 5 days in March, and again in September).
;     18-Jan-2005  DG
;       Inserted another 15 years worth of solar data and changed the
;       file name to 'sun1994-2020.dat'.
;-

function rd_soleph,jd

   ; Open the ephemeris file
   openr,/get_lun,lun,!defaults.ephemdir+'sun1994-2020.dat'
   line = ''
   lines = strarr(3)

   ON_IOERROR,bail
   ; Read header to get it out of the way
   readf,lun,lines

   ; Read a data line
   readf,lun,line

   ; Loop until Julian date is found
   while (strmid(line,1,7) ne long(jd)) do begin
      readf,lun,line
   endwhile

   ; Put this line into the output array and read two more lines
   lines(0) = line
   readf,lun,line
   lines(1) = line
   readf,lun,line
   lines(2) = line

   ; Close the file and free the logical unit number
   free_lun,lun
   return,lines
bail:
   return,-1
end

function get_soleph,jdin

   ; Convert Julian Date from current system time to JD of same date, but at 0 UT
   jd = double(long(jdin - 0.5))+0.5

   ; Open and read the solar ephemeris lines for this date
   lines = rd_soleph(jd)

   ; Return -1 if an error
   if (n_elements(lines) eq 1) then return,-1

   ; Get a copy of a SOLAR structure
   solar = get_sol_struct()

   ; Fill in SOLAR structure using data from Ephemeris file
   for i = 0, 2 do begin
      reads,lines[i],rah,ram,ras,decd,decm,decs,r,$
          format='(25x,i2,1x,i2,1x,f7.4,2x,i3,1x,i2,1x,f6.3,4x,f11.9)'
      ; Solar radius (semi-diameter) in degrees
      semid = 0.266183/r
      ; Horizontal parallax, in degrees
      hp = 0.0024417/r

      ; Keep track of sign of declination
      sgn = +1
      if (strmid(lines[i],40,1) eq '-') then begin
         sgn = -1
         decd = abs(decd)
      endif

      ; Convert input variables to integers in units of 0.1 mdeg
      solar.ra(i) = long((rah+ram/60.+ras/3600.)*15.*10000.D)
      solar.dec(i) = long(sgn*(decd+decm/60.+decs/3600.)*10000.D)
      solar.rad(i) = fix(semid*10000.)
      solar.hp(i) = fix(hp*10000.)

      ; The PANDB routine requires the Modified Julian Date
      mjd = strmid(lines(i),0,11)-2400000.
      pandb,mjd,p,b,rad
      solar.pa(i) = long(p*10000.)
      solar.b0(i) = long(b*10000.)

      ; Determine day of year from date string
      datstr = strmid(lines(i),11,11)
      doy = dayofyr(strmid(datstr,0,4),strmid(datstr,5,3),strmid(datstr,9,2))
      solar.doy(i) = doy
   endfor
return,solar
end