%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function do_all=energy_balance_0d;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% a zero dimensional global energy balance model including ice albedo
% feedback, showing an animation of the bifurcation diagram.  for APM
% 115.  Eli tziperman.

% to produce the quicktime animation of the bifurcation diagram:
% $ export DISPLAY=""; matlab -nodesktop -nosplash -nojvm -r energy_balance_0d_plus_bif_diagram > out.txt

close all; clear all;

% albedo specification:
a_min=0.1;
a_max=0.9;
T_min=-40;
T_max=20;
% longwave specification:
emissivity=0.6;
sigma=5.67e-8; % J/(sec m^2 K^4)
S0=1366/4;

if 0
  % integrate the time dependent equation:
  figure
  S0=1366/4; % watts/m^2
  T_init=-10;
  [time,T]=integrate_energy_balance_eqn(T_init,S0,T_min,T_max,a_min,a_max,emissivity,sigma);
  hold on
  %pause
  T_init=-20;
  [time,T]=integrate_energy_balance_eqn(T_init,S0,T_min,T_max,a_min,a_max,emissivity,sigma);
  %pause
end

% plot albedo, terms in the energy balance, and RHS of energy balance equation:
if 0
  S0=1366/4; % watts/m^2
  steady_states=calc_steady_states(S0,T_min,T_max,a_min,a_max,emissivity,sigma);
  plot_energy_balance_terms_and_rhs(S0,T_min,T_max,a_min,a_max,emissivity,sigma,steady_states)
  pause(8)
end

if 0
  S0=1000/4; % watts/m^2
  steady_states=calc_steady_states(S0,T_min,T_max,a_min,a_max,emissivity,sigma);
  plot_energy_balance_terms_and_rhs(S0,T_min,T_max,a_min,a_max,emissivity,sigma,steady_states)
  pause(8)
end

if 0
  S0=6000/4; % watts/m^2
  steady_states=calc_steady_states(S0,T_min,T_max,a_min,a_max,emissivity,sigma);
  plot_energy_balance_terms_and_rhs(S0,T_min,T_max,a_min,a_max,emissivity,sigma,steady_states)
  pause(5)
end

animate_bif_diagram(S0,T_min,T_max,a_min,a_max,emissivity,sigma)
quit

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot_energy_balance_terms_and_rhs(S0,T_min,T_max,a_min,a_max,emissivity,sigma,steady_states)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ScreenSize is a four-element vector: [left, bottom, width, height]:
scrsz = get(0,'ScreenSize'); 
figure('Position',[scrsz(4)*0.1 scrsz(4)*0.35 scrsz(3)*0.8 scrsz(4)*0.5])

T_min_plot=-150; T_max_plot=100; dT=1;

i=0;
for T=T_min_plot:dT:T_max_plot
  i=i+1;
  alpha(i)=albedo(T,T_min,T_max,a_min,a_max);
  temperature(i)=T;
end
subplot(2,3,1)
plot(temperature,alpha);
% AXIS([XMIN XMAX YMIN YMAX]):
V=axis; axis([T_min_plot,T_max_plot,0, 1]);
h1=ylabel('albedo');
h2=xlabel('temperature');
h3=title(sprintf('S0=%gW/m^2',4*S0));
set([h1,h2,h3,gca],'fontsize',14)

% plot incoming*(1-albedo) vs outgoing radiation as function of temperature:
i=0;
for T=T_min_plot:dT:T_max_plot
  i=i+1;
  incoming(i)=S0*(1-albedo(T,T_min,T_max,a_min,a_max));
  outgoing(i)=longwave(T,emissivity,sigma);
  temperature(i)=T;
end
subplot(2,3,2)
plot(temperature,incoming,'r');
hold on
plot(temperature,outgoing,'b');
% AXIS([XMIN XMAX YMIN YMAX]):
V=axis; axis([T_min_plot,T_max_plot,V(3),V(4)]);
h1=ylabel('radiation (watts/m^2)');
h2=xlabel('temperature');
h3=title('red: incoming; blue: outgoing');
set([h1,h2,h3,gca],'fontsize',14)


% plot incoming*(1-albedo)-outgoing radiation as function of temperature:
subplot(2,3,3)
plot(temperature,incoming-outgoing,'r');
% AXIS([XMIN XMAX YMIN YMAX]):
V=axis; axis([T_min_plot,T_max_plot,V(3),V(4)]);
hold on
plot(temperature,zeros(length(temperature),1),'k')
plot(steady_states-273,[0 0 0],'ro');
h1=ylabel('radiation (watts/m^2)');
h2=xlabel('temperature');
h3=title('dT/dt=incoming-outgoing');
set([h1,h2,h3,gca],'fontsize',14)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function a=albedo(T,T_min,T_max,a_min,a_max)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% T in degree celsius

if T>T_max 
  a=a_min;
elseif T<T_min
  a=a_max;
else
  a=a_max-(a_max-a_min)*(T-T_min)/(T_max-T_min);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function LW=longwave(T,emissivity,sigma)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% outgoing long wave radiation
% T in celsius

LW=emissivity*sigma*(T+273)^4;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [time,T]=integrate_energy_balance_eqn(T_init,S0,T_min,T_max,a_min,a_max,emissivity,sigma);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tspan=[0:0.011:10];
T0=T_init;  % -10 or -20 to get the two equilibria
[time,T] = ode45(@(t,X) myrhs(t,X,S0,T_min,T_max,a_min,a_max,emissivity,sigma),tspan,T0);
plot(time,T)
h1=title('zero dimensional energy balance model');
h2=xlabel('time');
h3=ylabel('T');
set([h1,h2,h3,gca],'fontsize',14)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rhs=myrhs(t,T,S0,T_min,T_max,a_min,a_max,emissivity,sigma)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
incoming=S0*(1-albedo(T,T_min,T_max,a_min,a_max));
outgoing=longwave(T,emissivity,sigma);
rhs=incoming-outgoing;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function steady_states=calc_steady_states(S0,T_min,T_max,a_min,a_max,emissivity,sigma);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initialize:
steady_states=zeros(3,1);

% T<T_min:
T1=(S0*(1-a_max)/(emissivity*sigma))^0.25;
if imag(T1)~=0 || real(T1) > T_min+273; T1=NaN; end

% T>T_max:
T3=(S0*(1-a_min)/(emissivity*sigma))^0.25;
if imag(T3)~=0 || real(T3) < T_max+273; T3=NaN; end

% T_min<T<T_max: search numerically
T0=(T_min+T_max)/2;
options = optimset('Display', 'off'); % Turn off Display
[x,fval,exitflag]=fzero(@(T) rhs_at_mid_range_Temps(T,S0,T_min,T_max,a_min,a_max,emissivity,sigma),T0,options);
T2=x+273; if exitflag<0; T2=NaN; end

steady_states=[T1 T2 T3];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rhs=rhs_at_mid_range_Temps(T,S0,T_min,T_max,a_min,a_max,emissivity,sigma)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t=0;
rhs=myrhs(t,T,S0,T_min,T_max,a_min,a_max,emissivity,sigma);
if T>T_max || T<T_min; rhs=NaN; end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function animate_bif_diagram(S0,T_min,T_max,a_min,a_max,emissivity,sigma)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% calculate the bifurcation diagram first:
S_s=400/4;dS=0.2;S_e=7500/5; % watts/m^2
N=length(S_s:dS:S_e);
bif_plot=zeros(3,N); S_plot=zeros(N,1);
i=0;
for S0=S_s:dS:S_e
  i=i+1;
  steady_states=calc_steady_states(S0,T_min,T_max,a_min,a_max,emissivity,sigma);
  bif_plot(:,i)=steady_states;
  S_plot(i)=S0;
end
if 0
  figure
  plot(S_plot,bif_plot(1,:),'r.')
  hold on
  plot(S_plot,bif_plot(2,:),'g.')
  plot(S_plot,bif_plot(3,:),'b.')
  xlabel('S_0')
  ylabel('T')
end

MakeQTMovie('start','budyko-bifurcation-animation.qt')

figure
% animate next:
for S0=S_s:dS*100:S_e
  clf
  subplot(2,2,1)
  % plot bifurcation diagram, marking location of current bifurcation
  % parameter (S0) with a vertical line:
  axis([0 1500 100-273 450-273]);
  hold on
  plot(S_plot,bif_plot(1,:)-273,'b.')
  plot(S_plot,bif_plot(2,:)-273,'g.')
  plot(S_plot,bif_plot(3,:)-273,'r.')
  h1=xlabel('S_0 (watts/m^2)');
  h2=ylabel('Temperature (C)');
  h3=title('bifurcation diagram, (green=unstable)');
  set([h1,h2,h3,gca],'fontsize',14)
  V=axis;
  plot([S0,S0],[V(3) V(4)])
  hold off


  subplot(2,2,2)
  % plot incoming*(1-albedo) vs outgoing radiation as function of temperature:
  T_min_plot=-150; T_max_plot=180; dT=1;
  i=0;
  for T=T_min_plot:dT:T_max_plot
    i=i+1;
    incoming(i)=S0*(1-albedo(T,T_min,T_max,a_min,a_max));
    outgoing(i)=longwave(T,emissivity,sigma);
    temperature(i)=T;
  end
  plot(temperature,incoming,'r');
  hold on
  plot(temperature,outgoing,'b');
  % AXIS([XMIN XMAX YMIN YMAX]):
  axis([T_min_plot,T_max_plot,0,1500]);
  h1=ylabel('radiation (watts/m^2)');
  h2=xlabel('Temperature (C)');
  h3=title('red: incoming; blue: outgoing');
  set([h1,h2,h3,gca],'fontsize',14)
  hold off
  
  print -djpeg tmp.jpg
  MakeQTMovie('addimage','tmp.jpg')

end

MakeQTMovie finish
