%cloudfactor_Jessica.m
 function [c]=cloudfactor_Jessica(t_met,latitude,Sin,timestep);
% trying to calculate a cloud factor presuming you have measured SW rad
% Radiation.m describes older attempts at cloud factors...
% for timestep in hours, Sin in W/m^2
%% STEP 1 Potential Solar Radiation

[Y, M, D, H, MN, S] = datevec(t_met);
N= greg2julian(Y, M, D, H, MN, S);%day of year, julian day

S0=1360; %solar constant W/m^2

phi=latitude*2*pi/360; %latitude of location of interest, in radians

ds=23.45*sin(2*pi*(284+N)/365); %solar declination, in degrees
 %the above is from my code, based on UEB, which matches wikipedia
 
delta=ds*2*pi/360; %converted to radians.

hs=acos(-tan(phi).*tan(delta));

Q0=86400*S0.*(hs.*sin(phi).*sin(delta)+cos(phi).*cos(delta).*sin(hs))./pi;
 %the above is in J per m^2
 
 %Note, this is extraterrestrial radiation and so has not been attenuated
 %at all, real would depend on atmospheric transmissivity as well as cloudiness
 
 
%% STEP 2: Convert observed W per m^2 to daily J per m^2 and calculate
%% clearness index
L=length(Sin);
xs=24/timestep;
ys=int32(length(Sin)./(24/timestep));
S1=reshape(Sin,xs,ys);
S2=sum(S1)*3600*timestep;
t1=reshape(t_met,xs,ys);
Q1=reshape(Q0,xs,ys);
Q2=mean(Q1);

% so the clearness index, k, is the ratio of the observed solar flux to the
% potential solar flux at the outside of the atmosphere, Q0!
k=S2./Q2; %clearness index, one value per day.

%% STEP 3:  Convert clearness index to cloud factor
% Campbell [1985] suggests that c can be linearly
% interpolated between c = 1.0 at a clearness index of 0.4 for
% complete cloud cover (kcld) to c = 0.0 at a clearness index of
% 0.7 (kclr). Others have used values of 0.35 and 0.6 for kcld
% and kclr, respectively [Flerchinger, 2000; Xiao et al., 2006];
% the optimum values are unknown and may perhaps depend
% on the location. (quoted from Flerchinger et al. 2009)

% I'm going to vote for the more generous 0.35 and 0.6 limits
% lowlimit=0.35;
% highlimit=0.6;

lowlimit=0.2;
highlimit=0.8;

c2=ones(1,length(k)); %first, make everything cloudy
sunny=find(k>=highlimit);
c2(sunny)=0;
sunbreak=find(k>=lowlimit & k<highlimit);
c2(sunbreak)=1-((k(sunbreak)-lowlimit)./(highlimit-lowlimit));

% %% STEP 4: Check everything with a plot
% figure
% subplot(3,1,1)
% plot(t_met,Q0/10^6,'r')
% hold on
% plot(t1(1,:),S2/10^6,'b')
% ylabel('SW (MJ/day)')
% legend('Incident at top of atmosphere','Observed')
% axis([datenum(1997,12,1) datenum(1998,4,1) -inf inf])
% datetick('keeplimits')
% subplot(3,1,2)
% plot(t1(1,:),S2./Q2,'k')
% ylabel('clearness index')
% axis([datenum(1997,12,1) datenum(1998,4,1) -inf inf])
% datetick('keeplimits')
% subplot(3,1,3)
% plot(t1(1,:),c2,'k')
% ylabel('cloud factor')
% axis([datenum(1997,12,1) datenum(1998,4,1) -inf inf])
% datetick('keeplimits')

%% STEP 5: Assume we have the same cloud factor every day, but distribute
%% it according to the original timeseries
% note that right now only certain timesteps are hardcoded in.
if timestep==2
c1=[c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; c2];
elseif timestep==1
    c1=[c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; ...
        c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; c2];
elseif timestep==3
    c1=[c2; c2; c2; c2; c2; c2; c2; c2];
elseif timestep==0.5
    c1=[c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; ...
        c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; ...
        c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; ...
        c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; c2; c2];
else
    error('only 1/2 hour, 1 hour, 2 hour and 3 hour timesteps are programmed at this time')
end
c=reshape(c1,length(Sin),1);
c=c(1:L); %throw away the last ones

