function spacecraft_sza, spacecraft_name, jd
;
; Returns the solar zenith angle for the given spacecraft 
; (observer) and Julian Day Number(s).
;
; B. Knapp, 2004-01-29
;
; $Revision$ $Date$
;
; $Header$
;
; $Log$
;
; Constants:
  SPEED_OF_LIGHT_IN_VACUUM_KM_PER_S = 299792.458d0
  KM_PER_ASTRONOMICAL_UNIT = 149597870.66d0
;
; Print usage?
  if n_params() lt 2 then begin
    print,""
    print," Returns the Solar zenith angle for the given spacecraft "
    print," (e.g., 'uars', 'sorce') and Julian Date(s) (e.g, 2452640.5"
    print," for 2003 Jan 1.0). The argument spacecraft_name must be a"
    print," scalar string, but the argument jd (double) may be a "
    print," scalar or array."
    print,""
    print," Usage:"
    print,""
    print,"     corr = spacecraft_sza(spacecraft_name, jd)"
    return,""
  endif
;
; Get spacecraft ID
  norad_id = [21701, 27651]
  p = strpos('UARS  SORCE ',strtrim(strupcase(spacecraft_name),2))
  if p lt 0 then begin
    message, 'Unknown spacecraft: '+spacecraft_name, /info
    return,''
  endif
  spacecraft_id = norad_id[p/6]

; Get the ECI Solar position (km, km/sec)
  ephem, jd, 0, earthSunXyz, status, posnid=2, framid=1, coorid=2
  earthSunXyz = earthSunXyz*KM_PER_ASTRONOMICAL_UNIT
  for j=3,5 do earthSunXyz[j,*] = earthSunXyz[j,*]/86400.0d0
;
; Get the ECI observer position (km, km/sec)
  spacecraft_pv, spacecraft_id, jd, earthObsXyz, coorid=2, $
    tle_path='/mizar/proj/sorce/science_data_system/data/', $
    pv_uncertainty=dEarthObsXyz, interpolation_fraction=f, status=status
;
; Get obsSunXyz (km, km/sec)
  obsSunXyz = earthSunXyz - earthObsXyz

; cos(SZA) = dot(earthObsDir, obsSunDir)
  earthObsDir = earthObsXyz[0:2,*]
  earthObsMag = $
    (sqrt(earthObsDir[0,*]^2+earthObsDir[1,*]^2+earthObsDir[2,*]^2))[*]
  obsSunDir = obsSunXyz[0:2,*]
  obsSunMag = $
    (sqrt(obsSunDir[0,*]^2+obsSunDir[1,*]^2+obsSunDir[2,*]^2))[*]
  for j=0L,n_elements(obsSunMag)-1 do begin
    earthObsDir[0,j] = earthObsDir[0:2,j]/earthObsMag[j]
    obsSunDir[0,j] = obsSunDir[0:2,j]/obsSunMag[j]
  endfor
  cosSza = (earthObsDir[0,*]*obsSunDir[0,*] + $
            earthObsDir[1,*]*obsSunDir[1,*] + $
            earthObsDir[2,*]*obsSunDir[2,*])[*]
  return, acos(cosSza)*180.d0/!dpi
end