% function daisyworld(lower_L,upper_L)
% Tian, Yang <yangtian@fas.harvard.edu>, 20150429
% model that shows life can influence global temperature, based on initial work by
% Watson and Lovelock (1983), Lovelock (2001)
%
% lower_L fraction of current insolation energy at the start of daisyworld.
% upper_L fraction of current insolation energy at the end of daisyworld.
%
%
% -Black daisy albedo (black).
% -White daisy albedo (white).
% -Bare planet albedo (bare).
% -Horizontal heat flux (q).
% -Daisy death rate (death).
% -Present day solar luminosity (S).
% -Lower temperature limit for daisy growth (lowgrowthK).
% -Upper temperature limit for daisy growth (highgrowthK).
% -The minimum daisy area to start growth (startarea).
% -Size of each increment determined by (numba).
%
%
% e.g. daisyworld(0.5,17)

function [A_b,A_w,A_g,T_p]=daisyworld(lower_L,upper_L)
S=1000;
q=15;
black=.25;
bare=.5;
white=.75;
cloud =.8;
death=.3;
lowgrowthK=278;
highgrowthK=313;
startarea=0.01; %start area as a fraction of 1
numba=100; %number of time steps relative to insolation increase
%calculate number of time steps (numba per 3 billion years)
length=(upper_L-lower_L)*numba;
if length~=round(length)
length=round(length);
boof=abs(length);
else boof=abs(length);
end
if length>0
jump=(upper_L*S-lower_L*S)/(length-1);
E(1,1)=lower_L*S; %initialises the lower limit relative to todays insolation
for a=2:length
E(a,1)=E((a-1),1)+jump; %creates the column matrix of solar Energy
end
end
%allows the insolation to be decreasing throughout daisyworld
if length<0
jump=(lower_L*S-upper_L*S)/(boof+1);
E(1,1)=lower_L*S; %initialises the lower limit relative to todays insolation
for a=2:boof
E(a,1)=E((a-1),1)-jump; %creates the column matrix of solar Energy
end
end
if length<0
length=abs(length);
end
SB=5.669*10^-8; %sets the Stefan-Boltzman constant
%initialises the daisy area to 1% to allow growth to start
Ablack=startarea;
Awhite=startarea;
%begin iterations for each time step
for a=1:length
%calculates area of bare ground and planet albedo
Abare=1-Ablack-Awhite;
Ap=Abare*bare+Awhite*white+Ablack*cloud;
%calculate emission and local temperatures
Temission=((E(a,1)/SB)*(1-Ap))^0.25;
Trock=((E(a,1)/SB)*(1-bare))^0.25;
Tblack=q*(Ap-black)+Temission;
Twhite=q*(Ap-white)+Temission;
%calculate daisy growth- define growth and add previous daisy cover
growthb=Ablack*(beta(lowgrowthK,highgrowthK,Tblack)*Abare-death);
growthw=Awhite*(beta(lowgrowthK,highgrowthK,Twhite)*Abare-death);
Ablack=Ablack+growthb;
Awhite=Awhite+growthw;
%ensure that there is no negative/too much area (numerical procedure)
if Ablack<=startarea
Ablack=startarea;
end
if Ablack>1
Ablack=1;
end
if Awhite<=startarea
Awhite=startarea;
end
if Awhite>1
Awhite=1;
end
%store results
T_p(1,a)=Temission;
A_b(1,a)=Ablack;
A_w(1,a)=Awhite;
A_g(1,a)=1-Ablack-Awhite;
T_b(1,a)=Tblack;
T_w(1,a)=Twhite;
T_r(1,a)=Trock;
A_p(1,a)=Ap;
growth_b(1,a)=growthb;
growth_a(1,a)=growthw;
end
%plot the results
f_S=[lower_L+0.01:0.01:upper_L];
%figure(1)
%set(gca,'fontsize',13)
%plot(f_S,A_b,'k-')
%hold on;
%plot(f_S,A_w,'c-')
%hold on;
%plot(f_S,A_w+A_b,'g--')
%xlabel('Solar luminosity')
%ylabel('Daisy Area')
%legend('Black daisies','White daisies','Total area',2)
figure
subplot(2,1,1)
%hold on
set(gca,'fontsize',13)
plot(f_S,A_w,f_S,A_b,f_S,A_g,'linewidth',1.5)
legend('a_w','a_b','a_g')
ylabel('fractional area')
title('Area and temperature of a daisy planet')
xlim([0.5 1.6])
subplot(2,1,2)
%hold on
set(gca,'fontsize',13)
plot(f_S,T_p-273,'linewidth',1.5)
ylim([-10 70])
xlim([0.5 1.6])
hold on;
plot(f_S,T_r-273,'k--')
subplot(2,1,2)
xlabel('Solar luminosity')
ylabel('Temperature (^0C)')
%legend('Daisyworld','Lifeless world',2)

