function [Ra Rso theta] = R(time,Lat,Lon,Z)
% computes extraterrestrial and clear sky surface irradiance
% INPUT
% time: matlab time vector (HAS TO BE A COLUMN VECTOR!)
% Z: station altitude above sea level [m]
% Lat,Lon: Latitude and Longitude of station [degrees, West is negative]
% OUTPUT
% Ra: extraterrestrial (top of atmosphere) irradiance [W/m2]
% Rso: clear sky global horizontal irradiance (GHI) at the surface [W/m2]
% theta: solar altitude angle [^o]
% USAGE: time = datenum(2009,9,1,0:1:24,0,0)'; [Ra Rso theta] = R(time,32.881269,-117.232867,331*.3048);

dt = median(diff(time));
time_vec = datevec(time);
DoY = ceil( time - datenum(time_vec(1),1,1) );
ToD = time_vec(:,4)+time_vec(:,5)/60+time_vec(:,6)/3600; clear time_vec
Lz = 120; %  longitude of the local time meridian (degrees West), 120° for Pacific time zone 
Lon = - Lon;

% find net longwave radiation
GSC = 0.082; % MJ m -2 min-1 (solar constant)
% Ra: extraterrestrial radiation
% dr = inverse relative distance Earth-Sun (correction for eccentricity of Earth's orbit around the sun) 
% omega1 = Solar time angle 1/2 hour before omega, that is, at the beginning of period (radians) 
% omega2 = Solar time angle 1/2 hour after  omega, at end of period (radians) 
% fe = Station latitude (radians) 
% delta = Declination of the sun above the celestial equator (radians)
% theta = solar altitude 
dr = 1+0.033*cos(2*pi/365 * DoY); 
        b = 2*pi/364*(DoY-81);
delta = 0.409 * sin(2*pi/365*DoY-1.39); clear DoY
    Sc = 0.1645 * sin(2*b) - 0.1255 * cos(b) - 0.025*sin(b);
omega = pi/12 * ( ToD + 4/60*(Lz-Lon)+Sc -12  ) ; % CIMIS: t must be center of time interval; 
omega = pi/12 * ( ToD + 1/15*(Lz-Lon)+Sc -12  ) ; % ASCE
clear ToD

omega1 = omega - (dt/2*24)*pi/12;
omega2 = omega + (dt/2*24)*pi/12;
    phi = pi*Lat/180;


%     omegas = acos(-tan(phi).*tan(delta));
% Ra = 24*60/pi*GSC*dr .* ( omegas *sin(phi).*sin(delta) +cos(phi)*cos(delta).*sin(omegas) ); % FAO daily
theta = asin(sin(phi)*sin(delta) + cos(phi)*cos(delta).*cos(omega) )*180/pi; % CIMIS & ASCE
sin_theta = (omega2-omega1).*sin(phi).*sin(delta) + cos(phi)*cos(delta).*(sin(omega2)-sin(omega1)); % ASCE

Ra = 6*60/pi*GSC*dr .* ( (omega2-omega1) *sin(phi).*sin(delta) + cos(phi)*cos(delta).*(sin(omega2)-sin(omega1)) ); % CIMIS
Ra =12*60/pi*GSC*dr .* sin_theta; % ASCE
Ra (Ra<0) = 0;
% Rso  = clear sky total global radiation at the earth surface (MJ m-2 h-1) 
% Rspm = solar radiation (MJm-2 h-1) during the last hour for which ?>10° 
% Rsopm= clear sky solar radiation (MJ m-2h-1) during the last hour for which ?>10 ° 
Rso = Ra * (0.75+2e-5*Z);
Ra = Ra / 3600 * 1e6 / (dt*24); % convert MJ m-2h-1 to W/m2
Rso = Rso / 3600 * 1e6 / (dt*24);