function [S,Ta,sun_elevation, declination]=solar_radiation(solar_constant, doy, lat, hour, minute)
%Usage
%[S,Ta,sun_elevation, declination]=solar_radiation(solar_constant, doy, lat, hour, minute)
%Input
%solar_constant [1]   solar constant in W m-2, usually about 1367
%doy        [m]   day of year for m days
%lat        [m]   degree latitude, e.g. 30 for northern, -30 for southern hemisphere
%hour       [m]   solar time in 24h format
%minute     [m]   solar minutes
%
%Output
%S                estimated downwelling solar radiation in W m-2 
%Ta               transmissivity of atmosphere, using cos(90-elevation)
%sun_elevation    sun elevation from the horizon in degrees
%declination      latitude where sun is at 90 degree at solar noon
%Tom Schut, 2008

nd=length(doy);
interval=1;% 0.3; %0.3 means once every 20 mins
TC=4;    

rf=1/interval;
rad=pi()./180;
Tday=[];

hourmin=hour+minute./60;

for i=1:nd    
    declrad = - asin(sin(23.45 * rad) .* cos (2.* pi() .* (doy(i) + 10)./365));
    sinld= sin(rad .* lat(i)) .* sin(declrad);
    cosld= cos(rad .* lat(i)) .* cos(declrad);
    aob=sinld./cosld;
    sinb=sinld + cosld .* cos(2.*pi().*(hourmin(i)-12)/24);

    declination(i)=declrad ./ rad;
    sun_elevation(i)=asin(sinb) ./ rad;
    Ta(i)=cos( (90-sun_elevation(i)) .* rad);
    S(i)=solar_constant * Ta(i) * sinb .* (1+0.033.*cos(2*pi()*(doy(i)-10) ./ 365 ));             
end

    