Pro guvi_pvat,pvat,guvi_time,geom,guvi_doy=guvi_doy
; Provides PVAT data on GUVI time grid. simple linear interpolation
; INPUTS
; pvat - PVAT file data structure read with read_ncdf,pvatfilename,pvat.
; Note that if the GUVI time crosses a day boundary, the pvat
; structures should be combined into an array. For example,
; read_ncdf,pvatfilename_1,p1
; read_ncdf,pvatfilename_2,p2
; pvat=[p1,p2]
; this routine will operate correctly if pvat is given as an array
; and a day boundary is not crossed
; guvi_time - time of day in GUVI representation (milliseconds of day).
; May be an array. When the timed crosses a day boundary,
; there are 2 options. The user can correct for the day crossing,
; for example
; read_ncdf,l1bfilename,d
; guvi_time=d.time+(d.doy-d.doy[0])*86400d3
; or the guvi_doy keyword may be used (see below).
; OUTPUTS
; geom - structure containing PVAT data. Interpolated values for
; position_eci
; position_ecef
; position_lat_lon_hgt
; roll_pitch_yaw
; refer to PVAT definitions
; KEYWORDS
; guvi_doy - array same number of elements as guvi_time specifying
; the day of year of the observation. For example
; read_ncdf,l1bfilename,d
; guvi_time=d.time
; guvi_pvat,pvat,guvi_time,geom,guvi_doy=d.doy
; ***Note - the use of this keyword WILL FAIL for orbits that cross
; from one year to another. The user must correct the year boundary
; time in the guvi_time input***
;-----------------------------------------
; -Demajistre 02
; - modified 2/03 to handle day boundaries
; - modified 2/04 to handle zero padded pvat structures
;-----------------------------------------
; find julian time (in days) from pvat and guvi time
gpsday=julday(1,6,1980,0.,0.,0.)
np=n_elements(pvat.spacecraft_time)
tp=reform(pvat.spacecraft_time-pvat.leap_seconds,np)
if keyword_set(guvi_doy) then $
guvi_time=guvi_time+(guvi_doy-guvi_doy[0])*86400d3
;
p_eci=reform(pvat.position_eci_cis,3,np)
p_ecef=reform(pvat.position_ecef_cts,3,np)
p_llh=reform(pvat.position_lat_lon_hgt,3,np)
p_rpy=reform(pvat.roll_pitch_yaw,3,np)
; cut out elements where spacecraft time <= 0
ii=where(pvat.spacecraft_time gt 0)
np=n_elements(ii)
tp=tp[ii]
p_eci=p_eci[*,ii]
p_ecef=p_ecef[*,ii]
p_llh=p_llh[*,ii]
p_rpy=p_rpy[*,ii]
jdg=gpsday+ $
double(tp[0]+guvi_time/1.d3)/86400d
jdp=gpsday+double(tp)/86400.
; find bracketing pvat time range
i0=max(where(jdp lt min(jdg)))
i1=min(where(jdp gt max(jdg)))
jd1=jdp[i0:i1]
; position ECI
position_eci=[[interpol(p_eci[0,i0:i1],jd1,jdg)], $
[interpol(p_eci[1,i0:i1],jd1,jdg)], $
[interpol(p_eci[2,i0:i1],jd1,jdg)]]
position_ecef=[[interpol(p_ecef[0,i0:i1],jd1,jdg)], $
[interpol(p_ecef[1,i0:i1],jd1,jdg)], $
[interpol(p_ecef[2,i0:i1],jd1,jdg)]]
position_lat_lon_hgt=[[interpol(p_llh[0,i0:i1],jd1,jdg)], $
[interpol(p_llh[1,i0:i1],jd1,jdg)], $
[interpol(p_llh[2,i0:i1],jd1,jdg)]]
roll_pitch_yaw=[[interpol(p_rpy[0,i0:i1],jd1,jdg)], $
[interpol(p_rpy[1,i0:i1],jd1,jdg)], $
[interpol(p_rpy[2,i0:i1],jd1,jdg)]]
geom={position_eci:position_eci,position_ecef:position_ecef, $
position_lat_lon_hgt:position_lat_lon_hgt, $
roll_pitch_yaw:roll_pitch_yaw}
return
end