function [Sd] = dailyinsolation(doy, S0, phi)
% DAILY_INSOLATION  Daily-mean insolation at top of atmosphere.
% 
%   DAILY_INSOLATION(DAY, S0, PHI) returns the daily-mean insolation
%   at day of year DAY [1,...,364] for solar constant S0 and latitude
%   PHI [angle in radians].
  
delta   =declination_angle(doy);
phi=phi*pi/180;  
  h0   = real(acos(-tan(phi).*tan(delta)));
  Sd   = S0/pi * (h0.*sin(phi).*sin(delta) ...
		+ cos(phi).*cos(delta).*sin(h0));