;+
; NAME:
;   oc_star_radec.pro
;
; PURPOSE:
;   returns [ra, dec] for a star at a given et
;
; DESCRIPTION:
;   returns [ra, dec] for a star at a given et
;
; CALLING SEQUENCE:
;   [ra,dec] = oc_star_radec(et,star, dradec, radecerr)
;
; INPUTS:
;   et - ephemeris time (seconds TDB after J2000)
;   star - star structure or array of structures
; OPTIONAL INPUT PARAMETERS:
;
; KEYWORDS
; OUTPUTS:
;  If et is a scalar and star is a single structure then
;  return [ra, dec].
;
;  If et is an array and the star is a single structure
;  then return [ radec_0, radec_1, ... ] where radec_i is for et[i]
;
;  If et is a scalar and the star is an array of structures
;  then return [ radec_0, radec_1, ... ] where radec_i is for star[i]
;  
;  If both et and the star are arrays
;  then return [ radec_0, radec_1, ... ] where radec_i is for et[i], star[i]
;  (if they are different lengths, then do up to the minimum)
; 
;  dradec is the derivatives of cos(dec) ra and dec w/ respect to time
;  (radian/sec).  Note: cos(dec) d ra/dt, not dra/dt.
;
;  REVISON HISTORY:
;   2006 Jan 02 Leslie Young.  
;   2006 Mar 08 LAY.  Corrected docs on dradec, added 2-element et,
;                     added array of radecerr.
;-

function oc_star_radec, et, star, dradec, radecerr

    net = n_elements(et)
    nstar = n_elements(star)

    if net eq 1 and nstar eq 1 then begin
        if star.ra lt -998.d or star.dec lt -998.d then begin
            dradec = [-999.d, -999.d]
            radecerr = [-999.d, -999.d]
            return, [-999.d,-999.d]
        endif
        if n_elements(star.et) eq 1 then begin 
            dtra = et - star.et
            dtdec = et - star.et
        endif else begin
            dtra = et - star.et[0]
            dtdec = et - star.et[1]
        endelse
        if star.radot lt -998.d then radot = 0.d else radot = star.radot
        ra = star.ra + dtra * radot / cos(star.dec)
        if star.decdot lt -998.d then decdot = 0.d else decdot = star.decdot
        dec = star.dec + dtdec * decdot
        dradec = [radot, decdot]
        mas = !dpi/(180.d * 3600.d * 1.d3)
        pmrms = 100*mas/31557600.  ; 100 mas/year
        if star.raerr lt -998. then begin
            raerr = -999.d
        endif else begin
            if star.radot lt -998. or star.radoterr lt -998.d then begin
                radoterr = pmrms
            endif else begin
                radoterr =  star.radoterr
            endelse
            raerr = sqrt(star.raerr^2 + (dtra * radoterr)^2)
        endelse
        if star.decerr lt -998. then begin
            decerr = -999.d
        endif else begin
            if star.decdot lt -998. or star.decdoterr lt -998.d then begin
                decdoterr = pmrms
            endif else begin
                decdoterr =  star.decdoterr
            endelse
            decerr = sqrt(star.decerr^2 + (dtdec * decdoterr)^2)
        endelse
        radecerr = [raerr, decerr]
        return, [ra, dec]
    endif

    if nstar gt 1 and net eq 1 then begin
        radec = dblarr(2,nstar)
        dradec = dblarr(2,nstar)
        radecerr = dblarr(2,nstar)
        for i = 0L, nstar - 1 do begin
            radec[*,i] = oc_star_radec(et,star[i], dradec1, radecerr1)
            dradec[*,i] = dradec1
            radecerr[*,i] = radecerr1
        endfor
        return, radec
    endif

    if nstar eq 1 and net gt 1 then begin
        radec = dblarr(2,net)
        dradec = dblarr(2,net)
        radecerr = dblarr(2,net)
        for i = 0L, net - 1 do begin
            radec[*,i] = oc_star_radec(et[i],star, dradec1,radecerr1)
            dradec[*,i] = dradec1
            radecerr[*,i] = radecerr1
        endfor
        return, radec
    endif

    if nstar gt 1 and net gt 1 then begin
        n = min( [nstar, net] )
        radec = dblarr(2,n)
        dradec = dblarr(2,n)
        radecerr = dblarr(2,n)
        for i = 0L, n - 1 do begin
            radec[*,i] = oc_star_radec(et[i],star[i], dradec1,radecerr1)
            dradec[*,i] = dradec1
            radecerr[*,i] = radecerr1
        endfor
        return, radec
    endif

end