Pro timed_sun_angle,pvat,guvi_time,azimuth,zenith, $
                    coordinate_system=coordinate_system
; calculates sun position at a given time from the PVAT file in various
; coordinate systems
; INPUTS
;     pvat - PVAT file data structure read with read_ncdf,pvatfilename,pvat
;     guvi_time - time of day in GUVI representation (milliseconds of day). May be an array.
; OUTPUTS
;     azimuth - azimuth angle (in degrees) of sun in specified coordinates
;     zenith - polar (zenith) angle (in degrees) of sun in specified coordinates
; Keywords
;     coordinate_system - specifies coordinate system for output, can have the following
;			values
;	0 (or absent)  - spacecraft coordinates
;	1 - nominal spacecraft coordinates (determined by sat. position, velocity and yaw state)
;	2 - East north up coordinates
;-----------------------------------------
  if not keyword_set(coordinate_system) then coordinate_system=0
; find julian time (in days) from pvat and guvi time
  gpsday=julday(1,6,1980,0.,0.,0.)
  jd=gpsday+ $
     double(pvat.spacecraft_time[0]-pvat.leap_seconds[0]+ $
            guvi_time/1.d3)/86400d
  jdp=gpsday+double(pvat.spacecraft_time-pvat.leap_seconds[0])/86400.
; find bracketing pvat time range
  i0=max(where(jdp lt min(jd)))
  i1=min(where(jdp gt max(jd)))

  jd1=jdp[i0:i1]
  njd1=n_elements(jd1)
; sun postion in eci
  aa_sun,jd1,sun_eci
  svec_eci=transpose([[sun_eci.x],[sun_eci.y],[sun_eci.z]])
  case coordinate_system of
     0: rmat=q2mat(pvat.quaternion[*,i0:i1])
     1: begin
        rmat=dblarr(3,3,njd1)
        for ic=0,njd1-1 do  $
           rmat[*,*,ic]=transpose( $
           nor2eci(pvat.position_eci_cis[*,i0+ic], $
                   pvat.velocity_eci_cis[*,i0+ic], $
                   svec_eci[*,ic]))
     end
     2: begin
        rmat=dblarr(3,3,njd1)
        zhat=-pvat.position_eci_cis[1,*]/total(pvat.position_eci_cis^2,1)
        for ic=0,njd1-1 do  begin
           zhat=pvat.position_eci_cis[*,i0+ic]/ $
                sqrt(total(pvat.position_eci_cis[*,i0+ic]^2))
           xhat=crossp([0.,0.,1.],zhat)
           xhat=xhat/sqrt(total(xhat^2))
           yhat=crossp(zhat,xhat)
           rmat[*,*,ic]=transpose([[xhat],[yhat],[zhat]])
        endfor
     end
     else: begin
        print,'TIMED_SUN_ANGLE: coordinate option not implemented: '+ $
              string(coordinate_system)
        return
     end
  endcase
  svec=svec_eci*0.
  for ic=0,njd1-1 do svec[*,ic]=rmat[*,*,ic]#svec_eci[*,ic]
  zenith1=!radeg*acos(svec[2,*]/sqrt(total(svec^2,1)))
  azimuth1=!radeg*atan(svec[1,*],svec[0,*])
  zenith=interpol(zenith1,jd1,jd)
  azimuth=interpol(azimuth1,jd1,jd)
  return
end