For loop with variables that are dependent upon other variables (Daisyworld numerical model)

I am attempting to produce a numerical model to simulate the growth of daisies on Daisyworld.

I have a for loop that should calculate the variables as time changes. These variables are used to define otehr variables through equations, however, one of my initial variables does not change as it should with each iteration of the loop.

My Albedo_planetary variable is used to define Te and Tw,Tb, but it does not change with each iteration.

Albedo_planetary = (white_albedo * white_area) + (black_albedo * black_area) + (bare_albedo * bare_area)

The initial variables are defined before the loop however the next iteration of them is defined at the end of the loop and as such Albedo_planetary does not change. I think I may need to use nested for loops but I am not sure how to implement that.

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% 
 %initial coverage of black daisies = 0.01
 %initial coverage of white daisies = 0.01
bare_area = 0.98; %initial coverage of bare ground = 0.98 
L = 1;%set  luminosity 
t_start = 0
t_end = 50
tstep = 1
t=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 array
black_area=zeros(1,length(t)); %preallocate array
bare_area = zeros(1,length(t)); %preallocate array
T4w = zeros(1,length(t)); %preallocate array
BrW = zeros(1,length(t)); %preallocate array
Tb = zeros(1,length(t));
BrB= zeros(1,length(t));
T4e = zeros(1,length(t)); %preallocate array
Albedo_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;


bare_area(1) = P - black_area(1) - white_area(1)

Albedo_planetary(1) = (white_area(1) * white_albedo) + (black_area(1) *black_albedo) + (bare_area(1) * bare_albedo)
T4e(1) = ((SL) / (SB ) * (( 1 - Albedo_planetary(1)))) ^0.25
T4w(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
if (T4w(1) - TempOpt < k^ -0.5) BrW(1) = 0;
    end
%BrB(1)= 0;
 BrB(1) = ((1- k) * (Tb(1) - TempOpt)) ^2
 if (Tb(1) - TempOpt < k^-0.5) BrB(1) = 0;
    end
 
    
%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)
    white_area(2) = white_area(1) + dfw_dt(1) * tstep
    black_area(2) = black_area(1) + dfb_dt(1) * tstep
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) / (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
   bare_areanew(i) = P - black_area(i) - white_area(i)
   Albedo_planetnew = (white_albedo * white_area(i)) + (black_albedo + black_area(i)) +(bare_albedo * bare_areanew(i))
         
   
end
Real_whitearea = real(white_area)
Real_blackarea = real(black_area)
figure (1)
plot(t,white_area)
hold on
plot(t,black_area)

Any help would be greatly appreciated

Answer 1

[SaiDileep Kola]

check the discussion about nested loops here