% ZSUN_M Program
% Calculates Solar Zenith Angle vs time of day
% Calculates Solar Irradiance and Total Daily Sunshine above Atmosphere
pi=atan(1)*4;
drc=pi/180;
sday=86400;%                Seconds per 24 hours                  
dt=300;
Sol=1370;%                  Solar Constant
dSEavg=149.5;%              Mean Sun-Earth Distance
ddSE=2.5;%                  +/- Variation of Sun-Earth Distance
tilt=23.45;%                Tilt of Earth's Axis to its Orbit
d_aphel = 13;%              Days from Solstice to Aphelion
% CAN MAKE LOOPS FOR LATITUDE(phi) AND DAY NUMBER (dsol) INSTEAD OF INPUTS
phi=input('Latitude = ').*drc;
dsol=(input('Month Number =')-6).*pi/6;
dsol=dsol+(input('Day Number of Month =')-21).*drc;
phisun=tilt.*drc.*cos(dsol);
daph=dsol-d_aphel.*drc;
dSE=dSEavg+ddSE.*cos(daph);
Fsun0=Sol.*(dSEavg/dSE).^2;
Qsun=0;
hold on
grid on
set(gca,'XLim',[0 24])
set(gca,'XTick',[0 6 12 18 24])
%Without this command, x axis tick marks occur at intervals of 5 hours
set(gca,'YLim',[0 1500])
set(gca,'YTick',[0 300 600 900 1200 1500])
for t =-sday/2:dt:sday/2
    csZ=sin(phi).*sin(phisun)+cos(phi).*cos(phisun).*cos(2.*pi.*t/sday);
    Fsun=Fsun0.*max(0,csZ);
    Qsun=Qsun+Fsun;
    plot(t/3600 + 12,Fsun,'o')
end
disp('Units for Qsun are Joules/m^2/day and for Favg are W/m^2')
Qsun=Qsun.*dt
Favg=Qsun/sday