function [MeanT, Tday, MeanTday, MeanTnight, d, n]=meandailytemp(doy,lat,Tmax,Tmin, daily)
%Determines mnean daily temperature including daylength, Gaussian curve
%during daytime and exponential decrease during night. See Gaudriaan's book
%
%Usage
%[MeanT, Tday, MeanTday, MeanTnight, d, n]=meandailytemp(doy,lat,Tmax,Tmin, daily)
%
%INPUT
%doy        [m]   day of year for m days
%lat        [m]   latitude
%Tmax       [m]   Maximum temperature for the day
%Tmin       [m]   Minimum temperature for the day
%daily            calculate Tday (1) or not (0). Slows process dramatically
%
%OUTPUT
%MeanT      [m]  Mean daiyly temperature
%Tday       [m*n]Temperature pattern from sunset to sunset
%MeanTday   [m]  Mean day time temperature
%MeanTnight [m]  Mean night time temperature
%d          [m]  Hrs of day
%n          [m]  Hrs of night
%
%Tom Schut, 2007

if nargin<5 daily=1; end

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

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

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

    %a day is defined as from sunset to sunset for conveniance
    d=12.0 .* (1 + 2 .* asin(aob) / pi() );
    if ~isreal(d) & Tmax(i) ~= Tmin(i)
        error('No daylight. Impossible to calculate diurnal pattern for given temp values');
    end
    
    n=24.00-d;
    sunr=12.00 - 0.5 .* d; 
    suns=12.00 + 0.5 .* d; 
    
    d=round(d./interval).*interval;
    n=round(n./interval).*interval;
    sunr=round(sunr./interval).*interval;
    suns=round(suns./interval).*interval;
    
    
    dtime=[interval : interval : d];
    ntime=[interval : interval : n];

    Tad=Tmin(i)+(Tmax(i)-Tmin(i)) .* sin( pi() .* ((dtime+sunr-12) + d ./ 2) ./ (d + 2 .* 1.5));
    Tsset=Tad(end);
    Tan=( Tmin(i) - Tsset .* exp(-n./TC) + (Tsset-Tmin(i)) .* exp(-ntime./TC) ) ./ (1 - exp(-n/TC));
    
    
    MeanTnight(i,1)=mean(Tan);
    MeanTday(i,1)=mean(Tad);
    if daily
        Tday(i,:)=[Tad Tan];
    end
    MeanT(i,1)= (d .* MeanTday(i,1) + n .* MeanTnight(i,1)) ./ (n+d);          
end

    