function theta = solarZenithAngle (hour, doy, lat)
% lat is in degree, (- 90, 90).
% doy: date of year, taking any value between 1 and 365
% hour: local hour, taking any value btw 0 and 24

    Lamda=lat*pi/180;
    Zeta=2*pi*(doy-1)/365;

% Delta is the sun's declination angle
    Delta=(0.006918-0.399912*cos(Zeta)+0.070257*sin(Zeta) ...
        - 0.006758*cos(2*Zeta) +0.000907*sin(2*Zeta) ...
        -0.002697*cos(3*Zeta) + 0.00148*sin(3*Zeta));
    cosTheta=sin(Lamda).*sin(Delta)+cos(Lamda).*cos(Delta).*cos(0.2618*(hour-12));
    theta=acos(cosTheta)*180/pi;

