I am currently trying to run a numerical model of Daisyworld. I am using the values from the original Lovelock paper.
When I run my code to calculate growth of the daisies, it returns complex numbers for my variables rather than real numbers. What should I do in order to get real numbers? Any help in solving this would be appreciated.
Here is my code:
clear % clear the workspace
%set initial conditions and constants
q = 2.06*10^9; %heat flux coefficient k-4
white_albedo = 0.75; %set white albedo
black_albedo = 0.25; %set black albedo
bare_albedo = 0.5; %set bare albedo
SB = 5.67e-8; % Stefann Boltzmann constant
Solar_lumin = 917; %solar luminosity (W m-2)
TempOpt = 295.65 ; %this is the optimum temperature for daisy growth in celsius 22.5/ 298
B = 0.25; %constant used in analytical solution
k = 17.5^-2; %width parametetr of daisy temperature growth 290.65K/17.5 celsius
qt= 79.79; % rescaled q
DRate = 0.3; %death rate (gamma)
initial_coverage = 0.01; %set initial coverage of B/W daises as 1%
black_area = initial_coverage; %initial coverage of black daisies = 0.01
white_area = initial_coverage; %initial coverage of white daisies = 0.01
bare_area = 0.98; %initial coverage of bare ground = 0.98
L = 1;%set luminosity
t_start = 0t_end = 10tstep = 1t=t_start:tstep:t_end; %time series of 0-100
SL = Solar_lumin *L P = 1; %proportion of fertile ground set to unity
%preallocate arrays for changing variables
dfw_dt= zeros(1,length(t)); %preallocate array
dfb_dt = zeros(1,length(t));white_area=zeros(1,length(t)); %preallocate arrayblack_area=zeros(1,length(t)); %preallocate arraybare_area = zeros(1,length(t)); %preallocate arrayT4w = zeros(1,length(t)); %preallocate arrayBrW = zeros(1,length(t)); %preallocate arrayTb = zeros(1,length(t));BrB= zeros(1,length(t));T4e = zeros(1,length(t)); %preallocate arrayAlbedo_planetary = zeros(1,length(t)); %preallocate array%set initial valies of changing variables
%dfw_dt(1) = 0;
%dfb_dt(1) = 0;
white_area(1) = 0.01;black_area(1) = 0.01;bare_area(1) = 0.98;%T4w(1) = 0;
Albedo_planetary(1) = (white_area(1) * white_albedo) + (black_area(1) *black_albedo) + (bare_area(1) * bare_albedo)T4e(1) = ((SL) / (4*SB ) * (( 1 - Albedo_planetary(1)))) ^0.25T4w(1) = (q*(Albedo_planetary(1) - white_albedo) + (T4e(1)^4)) ^(0.25)%Tb(1) = 0
Tb(1) = ((q * (Albedo_planetary(1) - black_albedo) + (T4e(1)^4)))^(0.25)%BrW(1) = 0;
BrW(1) = (1 - k) * ((T4w(1) - TempOpt))^2%BrB(1)= 0;
BrB(1) = ((1- k) * (Tb(1) - TempOpt)) ^2%T4e(1) = 0;
dfw_dt(1) = white_area(1) * (bare_area(1) * BrW(1) * T4w(1) - DRate)dfb_dt(1) = black_area(1) * (bare_area(1) * BrB (1) * Tb(1) - DRate)for i = 2:length(t) bare_area(i) = P - black_area(i-1) - white_area(i-1) Albedo_planetary(i) = (white_albedo * white_area(i-1)) + (black_albedo * black_area(i-1)) + (bare_albedo * bare_area(i-1)) T4e(i) = ((SL) / (4*SB ) * (( 1 - Albedo_planetary(i-1)))) ^(0.25) T4w(i) = ((q *(Albedo_planetary(i-1) - white_albedo) + ((T4e(i-1)).^4))) ^0.25 Tb(i) = ((q * (Albedo_planetary(i-1) - black_albedo) + (T4e(i-1).^4)))^0.25 BrW(i) = (1 - k) * ((T4w(i-1) - TempOpt))^2 if (T4w(i-1) - TempOpt < k^ 0.5) BrW(i) = 0 end BrB(i) = ((1- k) * (Tb(i-1) - TempOpt)) ^2 if (Tb(i-1) - TempOpt < k^0.5) BrB(i) = 0 end dfw_dt(i) = white_area(i-1) * (bare_area(i-1) * BrW(i-1) * T4w(i-1) - DRate) dfb_dt(i) = black_area(i-1) * (bare_area(i-1) * BrB (i-1) * Tb(i-1) - DRate) white_area(i) = white_area(i-1) + dfw_dt(i-1) * tstep black_area(i) = black_area(i-1) + dfb_dt(i-1) * tstep endfigure (1)plot(t,white_area)hold onplot(t,black_area)
Best Answer