%daisyworld.m gka 29 March 2008%
%PARAMETERS%
T_opt_b=295.0; %Optimal growing temperature of black daisies in K
T_opt_w=295.0; %Optimal growing temperature of black daisies in K
k_b=17.5^(-2); %growth rate bracketed between 30 K of T_opt
k_w=17.5^(-2); %growth rate bracketed between 30 K of T_opt
p=0.8; %Fraction of land in Daisy World that can support daisies
a_g=0.5; %Albedo of bare ground
a_w=0.9; %Albedo of white daisies
a_b=0.1; %Albedo of black daisies
S=1368/4; %Solar radiation in W/m2 received by a planet located 1 au from the Sun
L=3; %Luminosity
q=10^9; %Heat transfer coefficient to ensure thermal equilibrium
sigma=5.67*10^-8;% Stefan Boltzman constant J/s/m2/K4
gamma=.0001; %Death rate

%INITIALIZATION%
alpha_w=0.5; %Fraction of land covered by white daisies
alpha_b=0.0; %Fraction of land covered by black daisies
T=300; %Initial temperature of daisy world in K

%MAIN LOOP%
for itime=1:0.5*10^5
%	S=1368/4*(1+0.3*sin(2*pi*itime/10^4)); %insolation for problem 3
	A=alpha_w*a_w+alpha_b*a_b+a_g*(1-alpha_w-alpha_b); %mean planetary albedo
	
	alpha_w_store(itime)=alpha_w; %store variables for plotting
	alpha_b_store(itime)=alpha_b;
	temp_store(itime)=T;
	Astore(itime)=A;

	alpha_g=1-alpha_b-alpha_w; %amount of bare ground

	A=alpha_w*a_w+alpha_b*a_b+a_g*(1-alpha_w-alpha_b); %mean planetary albedo

	T=(S*L*(1-A)/sigma)^(1/4); %compute mean planetary temperature in radiative equilibrium
	T_w=(q*(A-a_w)+T^4)^(1/4); %compute temperature of patch of white daisies
	T_b=(q*(A-a_b)+T^4)^(1/4); %compute temperature of patch of black daisies
	T_g=(q*(A-a_g)+T^4)^(1/4); %compute temperature of bare ground
	%if mod(itime,10^3)==1     %next three lines are for problem 4
	%T_opt_b=T_opt_b+0.2;
	%end

	alpha_w=alpha_w+alpha_w*(alpha_g*betafn(T_w,T_opt_w,k_w)-gamma);
	alpha_b=alpha_b+alpha_b*(alpha_g*betafn(T_b,T_opt_b,k_b)-gamma);
end