%   EarthSWR_Model0a.m

%   Version:  30 September 2025 (25273)
%   Modified EarthSWR_Model0.m as follows:
%     - clarified the header section and added section breaks for debugging
%     - modified display of results, including figure numbering & titles
%     - added annual calculation of time and Earth-Sun distance calculation
%     - corrected TOASWR to include only the latitude of interest

%   Version:  01 May 2010  (0800)
%   Calculation Top of Atmos. (TOA) SWR Irradiance and orbital parameters
%   Code features:  Manual input of timeStart/timeEnd/SamplesPerHour  WORKING
%   Subroutines:  none

%% #########################################################################

%   Header Section
clear variables;
close all;

format long;
home;

%% #########################################################################

%   Input Parameters
TitleDateRange_Model = '17 March - 01 April';
year = 2010;

timeStart = 76;
timeEnd = 91;
SamplePerHour = 6;
t_offset_Model = 0;

% #########################################################################
% #########################################################################

%       Data Synthesis

%     Time:  Julian Day for the data window
NtModel = SamplePerHour*24*(timeEnd - timeStart) + 1;   %  Time with SamplePerHour resolution
timeModel = linspace(timeStart,timeEnd,NtModel);
fracdayModel = timeModel - floor(timeModel);

timeModelAnnual = linspace(1,365,SamplePerHour*24*365);

% -------------------------------------------------------------------------

Lat_NSB = 41.402*pi/180;   %  Latitude of New Science Building (NSB)

% #########################################################################

%       Calculation of TOA SWR Irradiance
%     Earth Orbital Parameters
SemiMajorAxis = 149.457e9;
eccentricity = 0.0167;
Torbit = 365.2564;

MeanAnomaly = 2*pi*(timeModel-3.0)/Torbit;
EarthSunDistance = SemiMajorAxis*(1-eccentricity*eccentricity)./(1+eccentricity*cos(MeanAnomaly));

SolarConstant = 1368;   %  TOA Solar Irradiance (annual mean value)
MeanEarthSunDistance = 149.600e9;   % Mean Earth Sun Distance (m);

% -------------------------------------------------------------------------

%     Special code for creating assignment Figure 1
MeanAnomalyAnnual = 2*pi*(timeModelAnnual-3.0)/Torbit;
EarthSunDistanceAnnual = SemiMajorAxis*(1-eccentricity*eccentricity)./(1+eccentricity*cos(MeanAnomalyAnnual));
MeanAnomalyAnnualBirthday = 2*pi*(22.5-3.0)/Torbit;
EarthSunDistanceAnnualBirthday = SemiMajorAxis*(1-eccentricity*eccentricity)./(1+eccentricity*cos(MeanAnomalyAnnualBirthday));

% #########################################################################

%     Calculation of cosine of the solar zenith angle for a single location.
phi_max = 23.45*pi/180;
SDA_Model = phi_max*cos(2*pi*(timeModel-173)/365);

CSZA_Model = sin(Lat_NSB)*sin(SDA_Model)-cos(Lat_NSB)*cos(SDA_Model).*cos(2*pi*(fracdayModel+t_offset_Model));
mu_oModel = zeros(1,NtModel);   %  Output quantity of interest, Preallocate for speed
for n = 1:NtModel;
    if CSZA_Model(n) < 0;
        mu_oModel(n) = 0.0;
    else
        mu_oModel(n) = CSZA_Model(n);
    end
end

%   Calculation of the Top of Atmosphere Solar Irradiance at a specified latitude
TOASWR = SolarConstant*mu_oModel.*(MeanEarthSunDistance./EarthSunDistance).^2;

%% ########################################################################

%   Plot Model data

figure(100)
hold on;
plot(timeModelAnnual,EarthSunDistanceAnnual,'k')
plot(22.5,EarthSunDistanceAnnualBirthday,'ro','MarkerFaceColor','r')
hold off;
box on;
xlabel(['Time (Julian Day, ',num2str(year),') UTC'])
ylabel('Earth/Moon - Sun Distance (m)')
title('Earth - Elliptical Orbit')

figure(101)
subplot(2,1,1)
hold on;
plot(timeModel,(180/pi)*SDA_Model,'r-')
refline(0,0);
hold off;
box on;
grid on;
ylabel('Solar Declination (degrees)')
title(['Orbital Mechanics & Solar Radiation:  ',char(TitleDateRange_Model),', Latitude = ',num2str((180/pi)*Lat_NSB),'^o'])

subplot(2,1,2)
plot(timeModel,mu_oModel,'b')
xlabel(['Time (Julian Day, ',num2str(year),') UTC'])
ylabel('cos(Zenith Angle)')

figure(102)
plot(timeModel,TOASWR,'r')
xlabel(['Time (Julian Day, ',num2str(year),') UTC'])
ylabel('Top of Atmosphere Irradiance (W m^{-2})')
title(['Orbital Mechanics & Solar Radiation:  ',char(TitleDateRange_Model),', Latitude = ',num2str((180/pi)*Lat_NSB),'^o'])
