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