%
%       sunearth.m
%
%   ....Simple script to illustrate sun-earth geometry.
%
%       Written by Von P. Walden using eqs. from S. Ackerman's notes.
%                  24 February 1998
%

% Sun-Earth distance
jday = 1:365;   % Julian days in non-leap year

r0 = 1;         % Astronomical unit
r = r0./sqrt(1 + 0.033*cos(2*pi*jday/365));

figure
plot(jday,r)
xlabel('Julian Day')
ylabel('Sun-Earth distance (AU)')

perihelion = find(r == min(r))
aphelion   = find(r == max(r))

% Declination angle
delta = 23.45*sin((2*pi/365)*(jday + 284));

figure
plot(jday,delta)
xlabel('Julian Day')
ylabel('Declination Angle (degs)')

delta_vernal_equinox = find(abs(delta(1:180)) == min(abs(delta(1:180))))             % Zero-crossing in first half of year.
delta_autumnal_equinox = find(abs(delta(181:365)) == min(abs(delta(181:365)))) + 181 % Zero-crossing in second half of year.
% delta_summer_solstice = 
% delta_winter_solstice = 

% Hour angle - 
t = 0:24;
w = 15*(t-12);

figure
plot(t,w)
axis([0 24 -180 180])
xlabel('Time (hours)')
ylabel('Hour Angle (degs)')


% Zenith angle
lat = 43;
lon = -89.5;
dec = delta([31+24]);      % Declination on 24 February.
w = 0;                     % Hour angle = 0 at noon.
theta0_Madison_24Feb = acos( sin(dec*pi/180)*sin(lat*pi/180) + cos(dec*pi/180)*cos(lat*pi/180)*cos(w*pi/180) ) *180/pi
theta0_Madison_Year  = acos( sin(delta*pi/180)*sin(lat*pi/180) + cos(delta*pi/180)*cos(lat*pi/180)*cos(w*pi/180) ) *180/pi;

figure
plot(jday,theta0_Madison_Year)
xlabel('Julian Day')
ylabel('Solar Zenith Angle (deg)')
title('For Madison, Wisconsin, 24 February, Noon')

% Sunrise hour angle and length of day.
ws = acos(-tan(lat*pi/180)*tan(dec*pi/180)) * 180/pi;
Nd = (2/15)*ws;