function zenith,jday,hour,lat,lon

;Purpose: 	Function to calculate the solar senith angle from known position and time of day.
;				Adapted from NASA code.
;
;Inputs:	jday: Day of the year
;			hour: Decimal hour of the day.
;			lat:  Decimal latitude (0 to +/-90, positive in N hemisphere)
;			lon:  Decimal longitude (0 to +/-180, positive west of the prime meridian
;
;Outputs:	result:  the solar zenith angle in degrees.
;
;I/O:		result=zenith(jday,dhour,lat,lon)
;
;Author:	Matthew Shupe
;Date:		06/05/01
;--------------------------------------------------------------------------------------------
sinob=0.3978
dpy=365.242
dph=15.0

rpd=!pi/180.0
dpr=1.0/rpd
dang=2.0*!pi*(jday-1.0)/dpy
homp=12.0 + 0.12357*sin(2.0*dang) - 0.004289*cos(dang) + 0.153809*sin(2.0*dang) + 0.060783*cos(2.0*dang)
hang=dph*(hour-homp)-lon
ang=279.9348*rpd+dang
sigma=(ang*dpr+0.4087*sin(ang)+1.8724*cos(ang)-0.0182*sin(2.0*ang)+0.0083*cos(2.0*ang))*rpd
sindlt=sinob*sin(sigma)
cosdlt=sqrt(1.0-sindlt^2)
cozena=sindlt*sin(rpd*lat)+cosdlt*cos(rpd*lat)*cos(rpd*hang)


zen=acos(cozena)*180.0/!pi
return,zen

end  ;function