% DAISYWORLD

% An imaginary desert world is populated only by white and black daisies.
% The white daisies love the heat of the sun, but, being white, they
% tend to reflect sunlight and cool down. The black daisies love coolness,
% but, being black, they tend to absorb more sunlight and heat up.
% This program simulates the struggle between the white and black daisies
% for space within a finite spatial domain as the solar flux to the planet 
% increases with time and each daisy seeks to equilibrate to the new conditions 
% - one time step behind.

% the constants:
SB_constant = 5.669E-8 %Stefan-Boltzmann constant W/¡K^4
Solar_Flux_Constant = 4000 %W/m2 for reference, our own solar constant=1368 W/m^2
uncovered_albedo = .5
white_albedo = .75
black_albedo = .25
heat_exchange_fact = 2.*10.^9 %this controls how the local temps differ from the global ave
death_rate=0.3

% initialize amount of area occupied by each daisy & desert
% initially there are few of each daisy
Uncovered_Area(1) = 1.0-0.002
White_Area(1) = 0.001
Black_Area(1) = 0.001
Solar_Luminosity(1) = 0.8 + 1.2/200. 

% as time progresses the Sun becomes hotter, the daisies respond
for itime = 2:30
    Solar_Luminosity(itime) = 0.8+(itime*(1.2/200)) ;
    planetary_albedo = (Uncovered_Area(itime-1)*uncovered_albedo)+(Black_Area(itime-1)*black_albedo)+(White_Area(itime-1)*white_albedo)
    % effective mean temperature calculated from global energy balance
    Avg_Planet_Temp(itime) = ((Solar_Luminosity(itime)*Solar_Flux_Constant*(1-planetary_albedo)/(4*SB_constant))^.25)
    % the local temperatures for each daisy population
    Temp_Black_Land(itime) = (heat_exchange_fact*(planetary_albedo-black_albedo)+Avg_Planet_Temp(itime)^4.)^0.25
    Temp_White_Land(itime) = (heat_exchange_fact*(planetary_albedo-white_albedo)+Avg_Planet_Temp(itime)^4.)^0.25
    % the growth rate for each daisy population
    Black_Growth_rate = 1-.003265*(295.5-Temp_Black_Land(itime))^2 
    White_Growth_rate = 1-.003265*(295.5-Temp_White_Land(itime))^2
    black_growth(itime) = Black_Area(itime-1)*(Uncovered_Area(itime-1)*Black_Growth_rate-death_rate)
    white_growth(itime) = White_Area(itime-1)*(Uncovered_Area(itime-1)*White_Growth_rate-death_rate)
    Black_Area(itime) = Black_Area(itime-1) + black_growth(itime)
    White_Area(itime) = White_Area(itime-1) + white_growth(itime)
    if (Black_Area(itime) < 0.0)
        Black_Area(itime)=0.001 ;
    end ;
    if (White_Area(itime) < 0.0)
        White_Area(itime)=0.001 ;
    end ;
    Uncovered_Area(itime) = 1.0 - Black_Area(itime) - White_Area(itime)
end

figure;
plot(Solar_Luminosity,Black_Area);
xlabel('solar luminosity') ; ylabel('Black Area')