zenith.pro

;+
; NAME:
; zenith
;
; PURPOSE:
; Calculate the solar zenith angle (and optionally the omega angle (see below))
; 
; CATEGORY:
; FUNCTION
;
; CALLING SEQUENCE:
; zenith(utc,lat,lon)
;
; EXAMPLE:
; zenith(201.35,20,10.)
;
; INPUTS: 
; utc flt or fltarr: the time (UTC) as days since Jan 1 00:00 of the year
;   e.g. the time 0.5 is noon of Jan 1st.
;   (presently it also works for the following year, if the first year
;    had 365 days. I.e it works for the NOXAR times for 1995 and 1996 but
;    not yet for 1997
; lat flt or fltarr: the latitude in deg (pos. for northern hemisphere)
; lon flt or fltarr: the longitude in deg (pos. for eastern longitudes)
;
; OPTIONAL INPUT PARAMETERS:
;
; KEYWORD INPUT PARAMETERS:
;
; OUTPUTS
; the solar zenith angle in radian
; omega : the hour angle in radian
;     the local time (in fractional hours) can be approximated by
;     calculating 12+omega*180/!pi*24/360.
; COMMON BLOCKS:
;
; SIDE EFFECTS: 
;
; RESTRICTIONS:
; 
; PROCEDURE:
; 
; MODIFICATION HISTORY:
; first implementation Jan, 17 1997 by Dominik Brunner
; adapted from the JNO2 program by Wiel Wauben, KNMI
;-

FUNCTION zenith,utc,lat,lon,omega=omega

; if utc,lat and lon do not have the same number of elements
; then resize the arrays to the minimal number of elements
dim=Min([n_elements(utc),n_elements(lat),n_elements(lon)])
if n_elements(utc) gt dim then utc=utc(0:dim-1)
if n_elements(lat) gt dim then lat=lat(0:dim-1)
if n_elements(lon) gt dim then lon=lon(0:dim-1)

; calculate the number of days since the 1.1. of the specific year
daynum=floor(utc)+1

; calculate the relative SUN-earth distance for the given day
; resulting from the elliptic orbit of the earth
eta=2.*!pi*daynum/365.
fluxfac=1.000110 $
  +0.034221*cos(eta)+0.000719*cos(2.*eta) $
  +0.001280*sin(eta)+0.000077*sin(2.*eta)
dist=1./sqrt(fluxfac)

; calculate the solar declination for the given day
; the declination varies due to the fact, that the earth rotation axis
; is not perpendicular to the ecliptic plane
delta=0.006918$
  -0.399912*cos(eta)-0.006758*cos(2.*eta)-0.002697*cos(3.*eta)$
  +0.070257*sin(eta)+0.000907*sin(2.*eta)+0.001480*sin(3.*eta)

; equation of time, used to compensate for the earth's elliptical orbit
; around the sun and its axial tilt when calculating solar time
; eqt is the correction in hours
et=2.*!pi*daynum/366.
eqt= 0.0072*cos(et)-0.0528*cos(2.*et)-0.0012*cos(3.*et)$
    -0.1229*sin(et)-0.1565*sin(2.*et)-0.0041*sin(3.*et)

; calculate the solar zenith angle
dtr=!pi/180. ; degrees to radian conversion factor
time=(utc+1.-daynum)*24 ; time in hours
omega=(360./24.)*(time+lon/15.+eqt-12.)*dtr
sinh=sin(delta)*sin(lat*dtr)+cos(delta)*cos(lat*dtr)*cos(omega)
solel=asin(sinh)
sza=!pi/2.-solel

return,sza

end