;+
; NAME:
;   oc_eph_et2state.pro
;
; PURPOSE:
;   returns [x,y,z,x',y',z'] for a eph at a given et
;
; DESCRIPTION:
;   returns [x,y,z,x',y',z'] for an ephemeris structure at a given et
;   where xyz are J2000 cartesian coordinates
;
; CALLING SEQUENCE:
;   state = oc_eph_et2state(et,eph)
;
; INPUTS:
;   et - ephemeris time (seconds TDB after J2000)
;   eph - array of eph structures
; OPTIONAL INPUT PARAMETERS:
;
; KEYWORDS
;  method
;    'piececub' - piecewise cubic interpolation of the state
;    'naif'     - calls to naif routines, using eph.naif
; OUTPUTS:
;  If et is a scalar, return [ra, dec].
;
;  If et is an array, return [ radec_0, radec_1, ... ] where radec_i is for et[i]
;
;  If et is a scalar and the eph is an array of structures
;  then return [ radec_0, radec_1, ... ] where radec_i is for eph[i]
;  
;  If both et and the eph are arrays
;  then return [ radec_0, radec_1, ... ] where radec_i is for et[i], eph[i]
;  (if they are different lengths, then do up to the minimum)
; 
;  REVISON HISTORY:
;   2006 Jan 02 Leslie Young.  
;-

function oc_eph_et2state, et, eph, method = method

    net = n_elements(et)
    neph = n_elements(eph)
    
    if net eq 1 then begin
      ; Determine index of data-points from which interpolation is made
        if neph eq 1 then begin
            i0 = 0
            i1 = 0
            ii = 0
        endif else begin
            i0 = value_locate(eph.et, et)
            i0 = 0 > i0 < (neph-2)
            i1 = i0 + 1
            ii = i0
            if (eph[i1].et-et) lt (et - eph[i0].et) then ii = 0
        endelse
        
        if not keyword_set(method) then method = 'naif'

        if method eq 'naif' then begin
            if total(strcmp(tag_names(eph),'naif', /fold)) ne 0 then begin
                naif = eph[ii].naif
                state = oc_naifstate(et, naif.targetCode,$
                                     obs=naif.observerCode,$
                                     abcorr=naif.abcorr,$
                                     frame=naif.frame,$
                                     lonlatalt=naif.lonlatalt)
                return, state
            endif
            ; we're here if there's no naif structure.
            ; next best is piecewise cubic interpolation
            method = 'piececub'
        endif
        
        if method eq 'piececub' then begin
            if neph gt 1 then begin
                state = dblarr(6)
                for j=0,2 do begin
                    piececub_interp, eph[i0:i1].et, $
                      reform((eph.state)[j,i0:i1]), $
                      reform((eph.state)[j+3,i0:i1]), $
                      et, xj, dxj
                    state[j] = xj
                    state[j+3] = dxj
                endfor
                return, state
            endif
            ; we're here if there's only one ephemeris point
            ; next best is linear extrapolation
            method = 'linext'
        endif

        if method eq 'linext' then begin
            s = eph[ii].state
            state = [s[0:2] + (et-eph[ii].et) * s[3:5], s[3:5] ]
            return, state
        endif
            
        
    endif
    
    state = dblarr(6,net)
    for i = 0, net - 1 do begin
        state[*,i] = oc_eph_et2state(et[i],eph, method=method)
    endfor
    return, state
end