% INSOLATE.m
% Script to calculate daily insolation
% AP download sept 2021 from: 
% https://pages.jh.edu/lhinnov1/paleoguide/tutorial1.html
% and not from the page indicated in the following article
% https://files.eric.ed.gov/fulltext/EJ1164362.pdf
% Article also saved and put as INSOLATE_README
% in https://people.mio.osupytheas.fr/~petrenko/TEACHING/MATLAB/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% AP explanation OUTPUT   sji = solar radiation
% sji    size 366x1440   days of the year *  minutes of the day (60*24 = 1440 mn)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

si=[ ];
sji=[ ];
dlghthrs=zeros(365);

% Factor to transform degrees to radians
degtorad=pi/180;

% Total input of solar radiation is 1370 W/m2
So=1370;

% Set up increment to do in a day
htot=1440;



% Total degrees in a day is 360
totdegs=360.;


% Define increment (dh)
	dh=(totdegs/htot)*degtorad;
	
% Select latitude (in degrees) to evaluate
theta=input('What latitude? (in degrees; negative for South)  ');
theta1=theta*degtorad;

%   AP following part commented to impose the current solar declination 
% in fact it varies with time from 22.2 to 24.5 (Berger et al., 1993)
decline=23.5
% % Set up declination values to do (seasonal range)
% decline=input('Input maximum declination (in degrees)  ');

decpos=decline*degtorad;
decneg=-1.*decpos;

% Calculate total declination range to do
dec=decpos-decneg;

% Evaluate declination every day of the half-year
ddec=dec/(365.256/2);
totdec=(decpos-decneg)/ddec;
ntotdec=round(totdec);
sum=0.0;

% Set up loop for each day of the first half-year
for J=1:ntotdec;
	d=(J-1)*ddec+decneg;

% Set loop over each hour in a day
for I=1:htot;
% Increment the declination throughout the day
    dd=d+ddec/htot;
	h=(I-1)*dh+pi;
	cz=sin(theta1)*sin(d)+cos(theta1)*cos(d)*cos(h);
	si(I)=So*cz;
	if si(I)<0.00000;
		si(I)=0.00000;  
		end
		sum=sum+si(I);
	if si(I)>0;
		dlghthrs(J)=dlghthrs(J)+24/htot;
	end
	end
	
% Store current calculations into day J
for I=1:htot;
sji(J,I)=si(I);
end
end	
%********
% Set up loop for each day of the second half-year
ntotdec2=2*ntotdec;
for J=ntotdec:ntotdec2;
	d=decpos-(J-ntotdec+1)*ddec;

% Set loop over each hour in a day
for I=1:htot;
% Increment the declination throughout the day
	dd=d-ddec/htot;
	h=(I-1)*dh+pi;
	cz=sin(theta1)*sin(dd)+cos(theta1)*cos(dd)*cos(h);
	si(I)=So*cz;
	if si(I)<0.000000;
		si(I)=0.000000;  
		end
		sum=sum+si(I);
	if si(I)>0;
		dlghthrs(J)=dlghthrs(J)+24/htot;
	end	
	end
	
% Store current calculations into day J
for I=1:htot;
sji(J,I)=si(I);
end
end	
smax1=max(sji);
smax=max(smax1);
hSo=max(sji);

% calculate total daylight hours
totdlghthrs=0.0;
for J=1:ntotdec2;
totdlghthrs=totdlghthrs+dlghthrs(J);
end

insoltot=sum;
fracyr=totdlghthrs/(24*365.256);

% Plot the output
ZMAX=So;
v=[0,smax];
%AP commented
% max(hSo)
% totdlghthrs
% fracyr
% insoltot


scrsz = get(0,'ScreenSize');
figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2]);
h = axes('Position',[0 0 1 1],'Visible','off');
axes('Position',[.25 .1 .7 .8])
surf(sji), view([-70,70]), axis([1 htot 1 365 0 smax]), caxis([0 1370]), colormap(jet), 
shading interp, colorbar, colordef black, whitebg('black')
set(get(gca,'xlabel'),'Rotation',65.0)
xlabel('Single Day (minutes) ','FontSize',11) 
ylabel('Day of Year  ','FontSize',11)
zlabel('Insolation (W/m2)  ','FontSize',11)
title('Incoming Solar Radiation','FontSize',20);
text(1480,700,['Max Insolation Value = ', num2str(max(hSo))]);
text(1380,690,['Total Annual Insolation = ', num2str(insoltot),' W/m2']);
text(1280,680,['Total Annual Daylight Hours = ', num2str(totdlghthrs)]);
text(1180,670,['Fraction of the Year With Daylight = ', num2str(fracyr)]);
text(1080,660,['Latitude (degrees) = ', num2str(theta)]);
text(980,650,['Declination (degrees) = ', num2str(decline)]);




