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