%% Neil Lareau
% EFD
% Mini-Project Part 1
% Radiation Subroutine

close all 
clear all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The purpose of this script is to calculate the incoming solar radiation
% at any point in the northern hemisphere, on any slope, and at any time of
% day and with varying cloud cover 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Solar Declination (ds)
%the solar declination should vary between -23.45 and 23.45 degrees
%depending on the time of year

%define the maximum elpictic tilt
tilt=23.45; %degrees (0.409 Radians)
dsolstice=173; % Julian day of the summer solstice
yearl = 365; %Julian days in a year
d=1:365 ; %day of year

ds=tilt.*cosd((360*(d-dsolstice))./yearl);

figure(1)
plot(d,ds);
xlabel('Day of the Year');
ylabel('Solar Declination');
title('Solar Declination as a Function of DOY');


%% Local Solar Elevation Angle

% this is a function of the time of day, longitude,latitude

%start by defining the time of day function

%we'll be using MST
hr=1:24;%hour of the day

hr=hr+6;
%hr(hr>24)=hr(hr>24)-24;
dayl=24;

%for now lets just hard code the lat/long
slclat=40.6518;%
%slclat=40.788;%dugway
slclon=111.9198;%positive to the west in this case
%slclon=111.978;% % dugway

beta=(360 * (hr)./dayl)-slclon;

%enter day of year


phi=asind((sind(slclat).*sind(ds(232))) - cosd(slclat).*cosd(ds(232)).*cosd(beta));                               
phi2=asind((sind(slclat).*sind(ds(264))) - cosd(slclat).*cosd(ds(264)).*cosd(beta)); 
phi3=asind((sind(slclat).*sind(ds(173))) - cosd(slclat).*cosd(ds(173)).*cosd(beta)); 
figure(2)
clf
plot(hr-6,phi);
xlabel('Hour of Day (MDT)')
ylabel('Solar Elevation Angle');
hold on;
plot(hr-6,phi2,'r')
plot(hr-6,phi3,'g')
%axis([10 28 0 90])


%% Local Azimuth Angle (angle of sun relative to true north)


alpha=acosd((sind(ds(232))-sind(slclat).*cosd(90-phi))./(cosd(slclat).*sind(90-phi)));

alpha((hr-7)>12)=360-alpha((hr-7)>12); %afternoon;


figure(3)
clf
plot(hr-6,alpha)

figure(4)
clf
plot(alpha,phi)



%% Slope calculation

betahat=0:1:90;%5; % deg slope
omegahat=[0 90 180 270]; % (degrees from north);
thetahat=zeros(length(betahat),length(alpha),4);

for j=1:4
    for i=1:length(betahat)
        thetahat(i,:,j)=acosd(cosd(betahat(i)).*cosd(90-phi) +sind(betahat(i)).*sind(90-phi).*cosd(alpha-omegahat(j)));
    end
end
%this is now the zenith angle that can be used in S=Snot*cos(theta)

%% Now we need to determine the incident solar radiation which is a
%% function of the amtospheric tramsimissivity
s0=1368; % TOA solar radiation

middle=0;
low=0;
high=0;

tr=(0.6+.2.*sind(phi)).*(1-.45.*high).*(1-.7.*middle).*(1-.4.*low);
tr2=repmat(tr,[91 1 4]);
rsw=s0.*tr2.*cosd(thetahat);
rsw(rsw<0)=0;



figure(5)
clf
subplot(2,2,1)
[cs,h]=contour(betahat,hr-6,rsw(:,:,1)')
clabel(cs,h);
axis([0 90 4 20])
xlabel('Slope Inclination (Deg)');
ylabel('Time of Local Day');
title('North Slope');
subplot(2,2,2)
[cs,h]=contour(betahat,hr-6,rsw(:,:,2)')
clabel(cs,h);
axis([0 90 4 20])
xlabel('Slope Inclination (Deg)');
ylabel('Time of Local Day');
title('West Slope');
subplot(2,2,3)
[cs,h]=contour(betahat,hr-6,rsw(:,:,3)')
clabel(cs,h);
axis([0 90 4 20])
xlabel('Slope Inclination (Deg)');
ylabel('Time of Local Day');
title('South Slope');

subplot(2,2,4)
[cs,h]=contour(betahat,hr-6,rsw(:,:,4)')
clabel(cs,h);
axis([0 90 4 20])
xlabel('Slope Inclination (Deg)');
ylabel('Time of Local Day');
title('East Slope');


%% Now to check against dugway data

dug=load('dugway.txt');

rsin=dug(:,8);

figure(6)
clf
plot(dug(:,7),rsin,'+');
%% modelled dugway RSW in

dugrsw=rsw(1,:,1);
dugrsw=squeeze(squeeze(dugrsw));
hold on;
plot(hr-6,dugrsw,'r')
xlabel('Hour')
ylabel('Incoming Short Wave [W/m^2]')
legend('Observed','Modeled');
title('Dugway Radiation')


%% Salt lake, oquirrhs, wasatch

index=find(hr-6==14);

contourset=rsw(:,index,:)
contourset=squeeze(contourset);

figure(7)
[c,h]=contour(omegahat,betahat,contourset)
clabel(c,h)
xlabel('Azimuth')
ylabel('Slope Angle')
title('Radiation as a Fucntion of Slope and Azimuth near SLC on Sept 21st');

%% VERY SIMPLE LONG WAVE OUT
b=-98.5; % W/m^2
rout= b.*(1-(.1.*high)-(0.3.*middle)-(0.6.*low));

%% NOW ALBEDO AND RNET

%mean albedo
albedo=0.18;
reflected= -albedo.*rsw(1,:,1);
rnet= rsw(1,:,1)+rout+reflected;
%rnet=squeeze(rnet);


%% NOW WE NEED TO PARTION THE SURFACE FLUXES (RESPONSE TERMS)


% we must first model the heat storage term for the boundary layer in order
% to arrive at the other fluxes

%define constant arrays
a1=[.34 .34 .14 .09 .7 .25 .38 .36];
a1(end+1)=mean(a1);
a2=[.31 .31 .33 .18 .33 .92 .56 .23];
a2(end+1)=mean(a2);
a3=[-31 -31 -6 -1 -.8 -22 -27.3 -19.3];
a3(end+1)=mean(a3);
f=[.193159 .103494 .296653 .003305 .00795 .017845 .041297 .066485 .269812];


%rate of change of the net radiation

dQdt=gradient(rnet);
%%

dQs=0;
for i=1:length(a1)
    dQs=dQs+ (f(i).*a1(i).*rnet) +(f(i).*a2(i).*dQdt(i))+(f(i).*a3(i));
end


%%


%start by loading the observations of T and P necessary for the
%psychrometric constant and the slope of the sat. vap. pres. curve


murray1=load('murray.txt');

%now we need to shift the times to correspond to the GMT times used above
%shift the data by 6 hours
 
murray=zeros(size(murray1));
murray(1:end-1,:)=murray1(2:end,:);
murray(end,:)=murray1(1,:);
%%

obtemp=murray(:,4);
obpres=murray(:,end);
obhl=murray(:,end-2);
obhs=murray(:,end-1);
obhr=murray(:,1);
obmin=murray(:,2)./60;
obdectime=obhr+obmin;





%calculate saturation vapor pressure as a function of obs temp 

lv=(2.501-0.00237.*obtemp).*10^6;
rv=461.5; 
cp=1004;
ep=.622;

claus=6.11*exp((5.417*10^3)*((1/273)-(1./(obtemp+273.15))));



%% approximate slope of the saturation vapor pressure curve

dedt=(lv./rv).*(claus./((obtemp+273.15).^2));
dedt=dedt(1:2:end);
% pschrometric constant

gamma=(cp./lv).*((obpres/10)./ep);
gamma=gamma(1:2:end);


%% now we need to select alpha and beta partioning constants for Qe and Qh

alpha=.5;
beta=20; %advection

fudge=2; % it seems like our simple parameterization above isn't storing enough heat into the ground


Qh=(((1-alpha)+(gamma./dedt))./(1+(gamma./dedt))).*(rnet'-(fudge.*dQs)') -beta;

Qe=(alpha).*((1)./(1+(gamma./dedt))).*(rnet'-(fudge.*dQs)') + beta;


figure(100)
clf
plot(rnet,'b');
hold on;
plot(rsw(1,:,1),'y')
plot(Qh,'r');
plot(Qe,'g');
plot(fudge.*dQs,'k');
plot(obhs(1:2:end),'r*')
plot(obhl(1:2:end),'g*')


%% Now load solar measured data from murray to compare

fid=fopen('murray_solar.txt');
mrry=textscan(fid,'%f %f %f %f %f %s %f %f %f %f %f %f %s %f %f %f','delimiter',',');
fclose(fid);
solar=mrry{15};
hrmrry=mrry{4};
minmrry=mrry{5};
time_mrry=hrmrry+minmrry./60;
time_mrry=time_mrry(end:-1:1);
solar=solar(end:-1:1);
figure(101)
clf
plot(time_mrry,solar)
hold on;
plot(hr-6,rsw(1,:,1),'r')
xlabel('MDT');
ylabel('W/m^2');
legend('Observed','Modeled')

%%