;+
; NAME:
;     DELAY
; PURPOSE:
;     Given the current time and hardcoded values for baseline info and
;     source position info (RA and Dec), calculate the time delay (TAU)
;     and its time derivative (TAUDOT) for each baselines 1-2, 1-4, and 2-4.
; CATEGORY:
;     OVRO APC INTERFEROMETRY
; CALLING SEQUENCE:
;     delay,time,tcycle,tau,taudot,geometry,uv
; INPUTS:
;     time      the current time [msec]
;     tcycle    the duration of a cycle [msec]
;     geometry  the geometry structure as returned from NEWSCAN.
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
; ROUTINES CALLED:
;     bdots, lst
; OUTPUTS:
;     tau       an NANT x NANT array containing the delay in nsec for
;                 each of the baselines
;     taudot    a corresponding array containing the rate of change of
;                 tau.
;     uv        an NANT x NANT array containing the u,v coordinates in
;                 nsec.  U is contained in the upper non-diagonal elements
;                 and V is in the lower non-diagonal elements.
; COMMENTS:
; SIDE EFFECTS:
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 20-Jun-1998 by Dale Gary
;     24-Jun-1998  DG
;       Made it somewhat more general by using baseline coordinates
;       from config segment structure
;     10-Jan-1999  DG
;       Made more general--it is clear that this routine will become a
;       standard one, so must be made fully precise.
;     01-Jul-1999  DG
;       Removed some parameters in the input list that were no longer
;       used.
;     01-Nov-1999  DG
;       Added UV argument, to return u,v coordinates calculated by
;       BDOTS for middle of interval.
;     28-Mar-2000  DG
;       Change to new LST routine, OVSA_LST, based on slalib routines
;       by Patrick Wallace.
;     15-Mar-2005  DG
;       Added code to deal with new GPS satellite observing scheme.
;-
pro delay,time,tcycle,tau,taudot,geometry,uv

   ; Some parameters
;   C2SIDSEC = 1.002738           ; Convert seconds to sidereal seconds
;   CBAM2MAS = 0.30174851D        ; Convert BAM to milliarcsec (masec)
;   CBAM2MS  = 0.020116568D       ; Convert BAM to millisec (msec)
;   CMM      = 2.997925D2         ; Speed of light [mm/nsec]
;   CSLAT    = 0.79619623         ; Cosine of OVRO latitude (37d13'53.8")
;   SNLAT    = 0.60503848         ; Sine of OVRO latitude (37d13'53.8")

   ; Time of the start of the data cycle, relative to the time of
   ; the first entry in the Ephemeris Table

   tnow = time + 64000.D  ; Add approx ET - UT1

   ; Calculate the RA, Dec, and HA at the current time, in absence of
   ; horizontal parallax.  Multiply powers of tnow by fit coefficients and
   ; use TOTAL() function to obtain sum.

   tarray = [1D,tnow,tnow^2]

   ; Convert date and time to single MJD double-precision number, for use in
   ; OVSA_LST() routine.
   mjd = cvdoy(geometry.year,geometry.doy[0],/JUL) + (time+geometry.dut)/8.64D7
;   hnow = lst(geometry.year,geometry.doy(0),time,geometry.dut) - ranow

   if (geometry.srcname eq 'GPS         ') then begin
      if (size(*geometry.psattab,/type) eq 5) then begin
         ; SATTAB is a double, so probably contains correct data
         tlist = (tnow-reform((*geometry.psattab)[0,*,0]))
         nlist = n_elements(tlist)
         later = where((tlist<0d) ne 0d,nlater)
         if (nlater ne 0) then inow = ((later-1)>0)[0] else inow = nlist-1
         tarray = [1D,tlist[inow],(tlist[inow])^2]
         geometry.a_ra = (*geometry.psattab)[1,inow,*]
         geometry.a_dec = (*geometry.psattab)[2,inow,*]
      endif
      ; For GPS satellite, RA is really HA)
      hnow = total(geometry.a_ra*tarray)
      decnow = total(geometry.a_dec*tarray)
      taunow = bdots(hnow,decnow,geometry,uvnow)
   endif else begin
      ranow = total(geometry.a_ra*tarray)
      decnow = total(geometry.a_dec*tarray)
      hnow = ovsa_lst(mjd) - ranow
      ; Time delay [msec] relative to current time
      taunow = bdots(hnow,decnow,geometry,uvnow)
   endelse

   ; Repeat these steps for next time
   tnext = tnow + tcycle
   tarray = [1D,tnext,tnext^2]

   ; Add TCYCLE, converted to fraction of a day
   mjd = mjd + tcycle/8.64D7
;   hnext = lst(geometry.year,geometry.doy(0),time+tcycle,geometry.dut) - ranext

   if (geometry.srcname eq 'GPS         ') then begin
      if (size(*geometry.psattab,/type) eq 5) then begin
         ; SATTAB is a double, so probably contains correct data
         tlist = (tnext-reform((*geometry.psattab)[0,*,0]))
         nlist = n_elements(tlist)
         later = where((tlist<0d) ne 0d,nlater)
         if (nlater ne 0) then inext = ((later-1)>0)[0] else inext = nlist-1
         tarray = [1D,tlist[inext],(tlist[inext])^2]
         geometry.a_ra = (*geometry.psattab)[1,inext,*]
         geometry.a_dec = (*geometry.psattab)[2,inext,*]
      endif
      ; For GPS satellite, RA is really HA)
      hnext = total(geometry.a_ra*tarray)
      decnext = total(geometry.a_dec*tarray)
      taunext = bdots(hnext,decnext,geometry,uvnext)
   endif else begin
      ranext = total(geometry.a_ra*tarray)
      decnext = total(geometry.a_dec*tarray)
      hnext = ovsa_lst(mjd) - ranext
      ; Time delay [msec] relative to current time
      taunext = bdots(hnext,decnext,geometry,uvnext)
   endelse

   ; Calculate delay and rate of change of delay at center of data cycle

   taudot = (taunext - taunow)/tcycle
   tau    = (taunow + taunext)*0.5
   uv     = (uvnow + uvnext)*0.5

return
end