function c = reedcf(yd,lat,dsw);
% REEDCF: computes daily mean cloud cover following Reed (1977).
% c = REEDCF(yd,lat,dsw) computes daily averaged cloud cover c from
% yearday, latitude, and measured insolation following Reed (1977), J. Phys. 
% Oceanog., 7, 482-485. Assumes hourly input series are either both 
% column or both row vectors of equal length. c is output as a boxcar 
% function with c constant over a 24 hr period.
%
%   INPUT:  yd  - yearday (e.g., Jan 10 is yd=10)
%           lat - latitude  [deg]
%           dsw - (measured) insolation  [W/m^2] 
%
%   OUTPUT: c - daily averaged cloud cover

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3/8/97: version 1.0
% 8/5/99: version 2.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% determine if input is either row or column vectors
[m,n] = size(dsw);

% first daily average swr and yd
davesw = binave(dsw,24);
daveyd = binave(yd,24);

% estimate daily averaged cloud factor
c=cloudfac(daveyd,lat,davesw);

% cloudfac always outputs a column vector, if initial input is
% row vector, convert cf to row vector
if m==1
  c  = c';
end

% convert daily averaged cloudfactor to boxcar function
if n== 1
  c = c*ones(1,24);
  c = reshape(c',prod(size(c)),1);

elseif m==1
  c = ones(24,1)*c;
  c = reshape(c,1,prod(size(c)));
end
c(length(dsw)+1:end)=[];


function cf=cloudfac(yd,lat,davesw)
% CLOUDFAC: computes daily mean cloud factor following Reed (1977).
% cf=CLOUDFAC(yd,lat,davesw) computes the daily average cloud factor 
% based on Reed (1977), J. Phys. Ocean., 7, 482-485.  Assumes year is 
% not a leap year. If yd is negative, then yd=yd+365. 
%
%   INPUT:  yd - yearday (e.g., Jan 10th is 10)
%           lat - latitude [deg]
%           davesw - average measured insolation  [W/m^2]
%
%   OUTPUT: cf - daily average cloud factor 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3/8/97: version 1.0
% 8/5/99: version 2.0 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% convert to column vectors
[N,M]=size(davesw);
if N==1
  davesw=davesw';
  yd=yd';
end

% check if yd is negative
ind=find(yd<0);
yd(ind)=yd(ind)+365;
% truncate to integer yearday for use with formula
yd=fix(yd);

% determine noon solar altitude
nsa=nsunang(yd,lat);

% compute ratio of observed to clear sky insolation
cssw=clskswr(yd,lat);
QdQ=davesw./cssw;

% correct for measurement error

ind1=find(QdQ>1);
QdQ(ind1)=ones(length(ind1),1);

% compute cloud factor
cf=(1-QdQ+.0019.*nsa)./.62;

% truncate limits

ind2=find(cf>1);
cf(ind2)=ones(length(ind2),1);
ind3=find(cf<0.3);
cf(ind3)=zeros(length(ind3),1);

% reconvert if necessary
if N==1
  cf=cf';
end


function nsa=nsunang(yd,lat)
% NSUNANG: computes noon solar altitude angle.
% nsa=NSUNANG(yd,lat) computes the noonday solar altitude angle as a 
% function of yearday and latitude using Smithsonian table (see 
% Reed (1976), NOAA Tech Memo., ERL PMEL-8, 20 pp., or Reed (1977), 
% J. Phys. Ocean., 7, 482-485).
%
%  INPUT:  yd  - yearday (e.g., Jan 10 is 10)
%          lat - latitude  [deg]
%
%  OUTPUT: nsa - noonday solar altitude  [deg]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3/8/97: version 1.0
% 8/5/99: version 2.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ensure that yd is positive
ind=find(yd<0);
yd(ind)=yd(ind)+365;
yd=fix(yd);

% compute noon solar altitude
k=pi./180;
t=2.*pi.*(yd./365);
d=.397+3.630*sin(t)-22.98.*cos(t)+.040.*sin(2.*t)-.388.*cos(2.*t)...
             +.075.*sin(3.*t)-.160.*cos(3.*t);
dl=k.*d;
ll=k.*lat;
z=sin(ll).*sin(dl)+cos(ll).*cos(dl);
nsa=asin(z)./k;