% This .m-file computes the difference, at a given day of the
% year and latitude, between current insolation and insolation in ~10000
% years. The insolation is computed at the top of the atmosphere. Variables
% ending with '1' are current and the ones ending with '2'
% are future.

close all;
clear all;

% Degrees latitude, setting up the vector (one dimensional array) phi
phi = (-90:1:90);

n = (1:1:365);   % Days of the year

% Make n and phi into two-dimensional arrays
[n,phi] = meshgrid(n,phi);

% The solar constant
S0 = 1367;

% The declination angle on day n of the year (see delta.m)
delta = delta(n);

% The hour angle at sunrise and sunset
h0 = acos(-tand(phi).*tand(delta));

% The squared ratio of mean Earth-sun to Earth-sun on day n See sqres.m
d1 = sqres(n);
d2= sqres(n+(365/2).*ones(size(d1)));

% Insolation at the top of the atmosphere
q1 = (S0/pi)*d1.*(h0.*sind(phi).*sind(delta) + cosd(phi).*cosd(delta).*sin(h0));
q2 = (S0/pi)*d2.*(h0.*sind(phi).*sind(delta) + cosd(phi).*cosd(delta).*sin(h0));

figure
contourf(q1,[0:50:600])
colorbar
title('Incoming radiation (W/m²)');
xlabel('Day of the year'); % 1 equals january 1st
ylabel('Latitude')

% Removing imaginary components obtained from the acos in the computation
% of h0
q1 = real(q1);
q2 = real(q2);
% Mean along the latitude dimension of q, note that this is different from
% Hartmann p. 31 fig. 2.7
m = mean(q1,2);
figure
plot(phi,m);
axis([-90 90 0 600]);
xlabel('Latitude');
ylabel('W/m²');
title('Annual insolation');

% The difference in insolation
dq = q1 - q2;

figure
contourf(dq,15)

title('Net radiation (W/m²)');
xlabel('Day of the year'); %1 equals january 1st
colorbar
