%% budyko_solver.m
% plot y_i and T_p as functions of fractional change in solar forcing 
% for the Budyko-type EBM with linear radiative heat loss and step-function
% albedo in temperature
%
% outline:
% -define parameters
% -define arrays for key variables
% -calculate insolation as a function of sine of latitude
% -calculate fractional change in solar forcing needed to sustain ice edge 
%     at a given latitude for partially ice-covered climate 
% -calculate stability of partially ice-covered climate states 
% -calculate minimum radiative flux to keep ice-free climate and maximum for snowball
% -plot results (making nice plots takes up much of the code) 
%
% Written by Tim Cronin, twcronin@mit.edu, 8/27-8/28/2018
%
% debugged/tested by: {XX}

% parameters
%name           value       units           comments
%--------------------------------------------------------------------------------------
A       =       201.6;      % W/m^2         OLR at Ts=0C; Roughly Budyko's a+n*a1 converted to SI
B       =       1.45;       % W/m^2/K       radiative damping; Roughly Budyko's B+n*b1 converted to SI
beta    =       3.79;       % W/m^2/K       climate heat transport coefficient; Roughly Budyko's \beta converted to SI
alpha_o =       0.32;       % -             albedo when ice-free
alpha_i =       0.62;       % -             albedo when ice-covered
Tm      =       -8;         % deg C         temperature threshold for ice
S0      =       1338;       % W/m^2         solar constant, perpendicular incidence
obliq   =       23.5;       % degrees       obliquity of orbit

yi      = (0:0.001:1);      % sine-latitude of ice edge
del_yi  = zeros(size(yi));  % fractional change in solar constant required to maintain ice line at yi
J_yi    = zeros(size(yi));  % solar absorption with ice-line at yi for reference S0
Tp_yi   = zeros(size(yi));  % global-mean temperature with ice line at yi
insol   = zeros(size(yi));  % annual-mean insolation as a function of latitude

% variable to loop over angular position of Earth in its orbit around the sun
theta   = 2*pi*(0.001:0.001:1);

% integrate over the year for each latitude to determine annual-mean
% insolation, insol(i) -- for formula see McGehee & Lehman, 2014, "A Paleoclimate Model of 
% Ice-Albedo Feedback Forced by Variations in Earth?s Orbit"
for j=1:length(yi)
    insol(j) = S0/(2*pi^2)*trapz(theta, sqrt(1-(sqrt(1-yi(j)^2)*sind(obliq)*cos(theta)-yi(j)*cosd(obliq)).^2));
end

% for each latitude, calculate J(y_i), delta(y_i), and T_p(y_i)
for j=1:length(yi)
    J_yi(j) = yi(j)*mean(insol(yi<=yi(j)))*(1-alpha_o) + ...
        (1-yi(j))*mean(insol(yi>=yi(j)))*(1-alpha_i); 
    % global-mean absorbed solar radiation with ice line at yi and delta=0

    del_yi(j)  = (A*(1+beta/B) + Tm*(B+beta))/...
        ((insol(j)*(1-(alpha_o+alpha_i)/2)) + beta*J_yi(j)/B)  - 1;
    % relative change in solar forcing required to keep partly ice-covered
    % climate with ice line at yi
    
    Tp_yi(j)   = 1/B*((1+del_yi(j))*J_yi(j) - A);
    % global-mean temperature for partly ice-covered
    % climate with ice line at yi
end

% calculate minimum possible solar forcing that will keep climate ice-free
% by setting polar T equal to Tm and using alpha_o everywhere
del_min_icefree = (A*(1+beta/B) + Tm*(B+beta))/...
    (insol(end)*(1-alpha_o) + beta*J_yi(end)/B)  - 1;

% then calculate maximum possible solar forcing that will allow a snowball
% climate by setting equatorial T equal to Tm and using alpha_i everywhere
del_max_snowball = (A*(1+beta/B) + Tm*(B+beta))/...
    (insol(1)*(1-alpha_i) + beta*J_yi(1)/B)  - 1;

% The stability condition of the Budyko model is that the ice line retreat 
% as the solar forcing increases, so filter
% calculations for that criterion
dyi_dlnS = diff(yi)./diff(del_yi);
dyi_dlnS = [-1 dyi_dlnS]; 
% diff() reduces the array size by 1, so pad with a leading negative value:
% model will always be unstable with the ice edge very close to equator;
% see Roe & Baker (2010), "Notes on a Catastrophe"

%% Make figure
figure('Position',[50 50 500 600]); % bracketed array gives size of figure, in pixels 

% first subplot: y_i against delta
subplot(2,2,1); % set up 2x2 array of subfigures, and select the first (top-left) panel
hold on;
set(gca,'Fontsize',12); % set the font size at 12 pt. YMMV depending on screen size
if sum(dyi_dlnS<0)>0
    plot(del_yi(dyi_dlnS<0), yi(dyi_dlnS<0),':','color',[1 0 0],'LineWidth',1); % unstable branch
end
if sum(dyi_dlnS>0)>0
    plot(del_yi(dyi_dlnS>0), yi(dyi_dlnS>0),'color',[0 0 1],'LineWidth',1.5); % stable branch
end
del_min_plot = min([min(del_yi)-0.1 del_min_icefree-0.1]); % calculate minimum value of delta for plot
del_max_plot = max([max(del_yi)+0.1 del_max_snowball+0.1]); % and max value of delta for plot
line([del_min_plot del_max_snowball],[0 0],...
    'color',[0 1 1],'LineWidth',1.5); % stable snowball branch
line([del_min_icefree del_max_plot],[1 1],...
    'color',[1 0 1],'LineWidth',1.5); % stable ice-free branch
line([0 0],[0 1],'color','k'); % line crossing equilibria with no solar perturbation
xlim([del_min_plot del_max_plot])
box on;
xlabel('fractional change in S_0');
ylabel('Ice-line sine latitude, y_i');
title('a) Ice line stability diagram', 'Fontweight','normal')

% second subplot: insolation against y
subplot(2,2,2)
hold on;
set(gca,'Fontsize',12);
plot(insol, yi,'color',[0 0 0],'LineWidth',1.5);
box on;
ylabel('sine latitude, y');
xlabel('insolation (W/m^2)')
title('b) Insolation against latitude', 'Fontweight','normal')

% third subplot: global-mean temperature T_p against delta
subplot(2,2,3)
hold on;
set(gca,'Fontsize',12);
if sum(dyi_dlnS<0)>0
    plot(del_yi(dyi_dlnS<0), Tp_yi(dyi_dlnS<0),':','color',[1 0 0],'LineWidth',1);
end
if sum(dyi_dlnS>0)>0
    plot(del_yi(dyi_dlnS>0), Tp_yi(dyi_dlnS>0),'color',[0 0 1],'LineWidth',1.5);
end
line([del_min_plot del_max_snowball],...
    [((1+del_min_plot)*mean(insol)*(1-alpha_i)-A)/B ...
    ((1+del_max_snowball)*mean(insol)*(1-alpha_i)-A)/B],...
    'color',[0 1 1],'LineWidth',1.5); % snowball climate branch
line([del_min_icefree del_max_plot],...
    [((1+del_min_icefree)*mean(insol)*(1-alpha_o)-A)/B ...
    ((1+del_max_plot)*mean(insol)*(1-alpha_o)-A)/B],...
    'color',[1 0 1],'LineWidth',1.5); % ice-free climate branch
Tmax_plot = ((1+del_max_plot)*mean(insol)*(1-alpha_o)-A)/B;
Tmin_plot = ((1+del_min_plot)*mean(insol)*(1-alpha_i)-A)/B;
box on;
xlabel('fractional change in S_0');
ylabel('global-mean temperature');
xlim([del_min_plot del_max_plot]);
ylim([Tmin_plot Tmax_plot]);
line([0 0],[Tmin_plot Tmax_plot],'color','k'); % line crossing equilibria with no solar perturbation
if sum(dyi_dlnS<0)==0
    hl = legend('Stable, partial ice-cover','Snowball','Ice-free');
elseif sum(dyi_dlnS>0)==0
    hl = legend('Unstable, partial ice-cover','Snowball','Ice-free');
else
    hl = legend('Unstable, partial ice-cover','Stable, partial ice-cover','Snowball','Ice-free');
end

set(hl,'Position',[0.5050    0.3567    0.3360    0.0958]);
title('c) T_p stability diagram', 'Fontweight','normal')