PRO ephem,jd,ut,gha,dec,eqtime

    ;INPUT: 
    ;jd  Julian day number
    ;ut  Universal Time in hours

    ;OUTPUT:
    ;gha     Greenwich hour angle, the angle between the Greenwich
    ;        meridian and the meridian containing the subsolar pt
    ;dec     solar declination, latitude of subsolar point
    ;eqtime  equation of time, longitudinal correction to mean sun
    ;        position

    ; REF:  Explanatory Supplement to the Astronomical Almanac
    ;       ISBN 0-935702-68-7
    ;       University Science Books
    ;       1992
    ;       pp. 484-485

    ;calculate number of centuries from J2000
    t = (jd + (ut/24.) - 2451545.0) / 36525.

    ;mean longitude corrected for aberration
    l = (280.460 + 36000.770 * t) MOD 360

    ;mean anomaly
    g = 357.528 + 35999.050 * t

    ;ecliptic longitude
    lm = l + 1.915 * SIN(g*!DTOR) + 0.020 * SIN(2*g*!DTOR)

     ;obliquity of the ecliptic
    ep = 23.4393 - 0.01300 * t


    ;equation of time
    eqtime = -1.915*SIN(g*!DTOR) - 0.020*SIN(2*g*!DTOR) $
            + 2.466*SIN(2*lm*!DTOR) - 0.053*SIN(4*lm*!DTOR)

    ;Greenwich hour angle
    gha = 15*ut - 180 + eqtime

    ;declination of sun
    dec = ASIN(SIN(ep*!DTOR) * SIN(lm*!DTOR)) * !RADEG

END

;----------------------------------------------------------------

PRO greg_to_jdaynum,y,m,d,jd

    ; Converts a Gregorian calendar date (month, day, year)
    ; to Julian day number (number of days since 1 Jan 4713 B.C.

    ; INPUT:  
    ; y   year
    ; m   month (1-12)
    ; d   day of month (1-31)

    ; OUTPUT: 
    ; jd  Julian day number

    ; REF:  Explanatory Supplement to the Astronomical Almanac
    ;       ISBN 0-935702-68-7
    ;       University Science Books
    ;       1992
    ;       p. 604

    t1 = (1461.*(y+4800+(m-14)/12.))/4.
    t2 = (367.*(m-2-12.*((m-14)/12.)))/12.
    t3 = -(3.*((y+4900+(m-14)/12.)/100.))/4.

    jd = t1+t2+t3 + d - 32075

END

;----------------------------------------------------------------

PRO get_sza,ut,year,month,day,lat,lon,sza

    ;Calculate solar zenith angles

    ;INPUT:
    ;ut      Universal Time in hours
    ;year    year (1000-9999)
    ;month   month of year (1-12)
    ;day     day of month (1-31)
    ;lat     geodetic latitude
    ;lon     geodetic longitude

    ;OUTPUT:
    ;sza     solar zenith angle, the angle between the local
    ;        vertical and the sun, in degrees

    ;get Julian daynumber, needed for the call to ephem
    greg_to_jdaynum,year,month,day,jd

    ;get equation of time and Greenwich hour angle
    ephem,jd,ut,gha,dec,eqtime

    ;calculate hour angle in degrees
    ha = lon + gha

    ;check range of ha 
    ndx=WHERE(ha LT 0,count)
    IF(count GT 0) THEN ha(ndx)=ha(ndx)+360.
    ndx=WHERE(ha GT 360.,count)
    IF(count GT 0) THEN ha(ndx)=ha(ndx)-360.

    ;calculate solar zenith angle in radians
    t1  = sin(dec*!DTOR)*sin(lat*!DTOR)
    t2  = cos(dec*!DTOR)*cos(lat*!DTOR)*cos(ha*!DTOR)
    sza = ACOS(t1 + t2)

    ;convert to degrees
    sza = sza * !RADEG

END