%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function doall=gaia
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
% Eli Tziperman: example for great papers course.
% Solve Lovelock's Daisy-world Gaia example using ode45
% see LaTeX notes for the equations.

% define parameters:
p=set_parameters;

% initial conditions:
X0=[0.3 0.3];
T0=calc_temperatures(X0,p)';

%% run models for different values of solar luminosity, controlled by
%% nondimensional parameter f_S:
f_S=0.5:0.01:1.7;
N=length(f_S);
Areas=zeros(N,3);
Temperatures=zeros(N,3);

for i=1:N
  %% set luminosity:
  p.S=p.S0*f_S(i);
  %% integrate:
  [t,X] = ode45(@rhs,[0 p.T],X0,p.options,p);

  %% plot one example of time evolution:
  if i==floor(N/2)
    figure(1); clf;
    b=1;
    e=length(t(:));
    plot(t(b:e),X(b:e,1),'b');
    hold on
    plot(t(b:e),X(b:e,2),'r');
    plot(t(b:e),1-X(b:e,1)-X(b:e,2),'g');
    h1=legend('a_w','a_b','a_g');
    h2=xlabel('time');
    h3=ylabel('areas');
    set([gca,h1,h2,h3],'fontsize',14);
  end
  
  %% store values of steady states:
  Areas(i,:)=[X(end,1) X(end,2) 1-X(end,1)-X(end,2)];
  Temperatures(i,:)=calc_temperatures(X(end,:),p)';
 
  
  %% set initial condition for next run:
  X0=X(end,:);
end

%% plot results as function of luminosity:
fig2=figure(2); clf
subplot(1,2,1)
plot(f_S,Areas(:,1),'b');
hold on
plot(f_S,Areas(:,2),'r');
plot(f_S,1-Areas(:,1)-Areas(:,2),'g');
h1=legend('a_w','a_b','a_g');
h2=xlabel('luminosity');
h3=ylabel('areas');
set([gca,h1,h2,h3],'fontsize',14);
xlim([min(f_S) max(f_S)]);

subplot(1,2,2)
plot(f_S,Temperatures(:,1),'b');
hold on
plot(f_S,Temperatures(:,2),'r');
plot(f_S,Temperatures(:,3),'g');
plot(f_S,Areas(:,1).*Temperatures(:,1)+Areas(:,2).*Temperatures(:,2)+...
    (1-Areas(:,1)-Areas(:,2)).*Temperatures(:,3),'k');
h1=legend('T_w','T_b','T_g','T_p','location','northwest');
h2=xlabel('luminosity');
h3=ylabel('Temperatures (K)');
set([gca,h1,h2,h3],'fontsize',14);
xlim([min(f_S) max(f_S)]);



%% save figure as pdf
set(gcf, 'PaperUnits', 'inches'); set(gcf, 'PaperSize', [6 4]);
set(gcf, 'PaperPosition', [0 0 6 4]); % [left, bottom, width, height];
saveas(gcf,sprintf('dasyworld-model-results.pdf'));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = rhs(t,X,p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get state variables (area fractions) from X:
a_w=X(1);
a_b=X(2);
a_g=1-a_w-a_b;

% Calculate temperatures by solving linearized equations A*T=Gamma:
T=calc_temperatures(X,p);
T_w=T(1); T_b=T(2); T_g=T(3);

% RHS of area fraction equations:
out=[a_w*(a_g*beta(T_w)-p.gamma) ...
     a_b*(a_g*beta(T_b)-p.gamma)]';


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function b=beta(T);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b=max(0,1-(T-273-22.5)^2/17.5^2);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function T=calc_temperatures(X,p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get state variables (area fractions) from X:
a_w=X(1);
a_b=X(2);
a_g=1-a_w-a_b;
A=[p.q*(a_w-1)-a_w*p.B, p.q*a_b, p.q*a_g; ...
p.q*a_w, p.q*(a_b-1)-a_b*p.B, p.q*a_g; ...
 p.q*a_w, p.q*a_b,p.q*(a_g-1)-a_g*p.B];
Gamma=[-a_w*p.S*(1-p.alpha_w)+a_w*(p.A-p.B*p.T_0), ...
       -a_b*p.S*(1-p.alpha_b)+a_b*(p.A-p.B*p.T_0), ...
       -a_g*p.S*(1-p.alpha_g)+a_g*(p.A-p.B*p.T_0)]';
T=A\Gamma;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p=set_parameters;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Integration time:
p.T=400.0;
% set error bounds for time integration:
p.options=odeset('RelTol',10^(-6),'AbsTol',[10^(-6)*ones(1,2)]);

% Stefan-Boltzmann constant:
p.sigma=5.670373e-8; % mks
% parameters for linearizing black body radiation law:
p.T_0=273;
p.A=p.sigma*p.T_0^4;
p.B=p.sigma*4*p.T_0^3;
% albedos:
p.alpha_w=0.75;
p.alpha_b=0.25;
p.alpha_g=0.5;
% solar constant:
p.S0=1340/4;
p.S=p.S0;
% transport coeff:
p.q=1.5*(4*p.sigma*p.T_0^3);
% death rate:
p.gamma=0.5;
