for loop with stored variables

Hi i have written a for loop for several variables that i have defined in other coding. I cannot seem to get matlab to recognise these vectors that i have already defined, I am relatively new to matlab and i am unsure if i have set up the for loop correctly. Any help would be great

this defines daisyworld 1

function dydt =daisyworld1(t,y,Lt,L)

    Lint=interp1(Lt,L,t); %interpolate the data set (Lt,L) at time t

    dydt = [0;0];

    A=((1-y(1)-y(2))*0.5)+(y(1)*0.25)+(y(2)*0.75);       %albedo

    S=917;         %constant solar energy

    Z=5.67*10^(-8);           %Stefan-Boltzmann constant;

    Te=((((S*Lint)/Z)*(1-A)).^(1/4))-273;       %plantary temperature

    B= 1-0.003265*((22.5-Te).^2);           %beta value,

    g=0.3;      %death term

    dydt(1)=y(1)*((  (1-y(1)-y(2))  *B)-g);        %black daisy formula

    dydt(2)=y(2)*((  (1-y(1)-y(2))  *B)-g); %white daisy formula

I then solve it using ode45

clear; % Remove stored variables daisyworldode45 

Lt=linspace(0,1000,25);    %generate t for L

L=Lt/500+ 0.4;               %luminosity function want to keep within a range of 0 and 2 so change depending on tspan 

tspan=[0 1000];             %solve for values of t

IC = [0.2 0.3];     %initial conditions of daisy percentage, black then white

[T, Y]=ode45(@(t,y)daisyworld1(t,y,Lt,L),tspan,IC);             % solves equation, need capital T and Y

I am trying to extract these values of T and Y, and use the vector T to find values for L and then plot these against Te

N=7269;
Te= zeros(1,N);

for T=0:N
    for Y=0:N
        Te(T,Y)=((917*(T/500+0.4)/(5.67*10^(-8))*(1-(((1-y(1)-y(2))*0.5)+(y(1)*0.25)+(y(2)*0.75))).^(1/4)))-273;
    end
end  
plot (T,Y);

7269 is the number of y values that ode45 uses i keep getting the error code

Undefined function 'y' for input arguments of type 'double'.

Error in daisyworldode45 (line 23)
        Te(T,Y)=((917*(T/500+0.4)/(5.67*10^(-8))*(1-(((1-y(1)-y(2))*0.5)+(y(1)*0.25)+(y(2)*0.75))).^(1/4)))-273;

[ChristianW]

What are you trying to calculate at line 23 (Te .. = ...)?

What does this mean: "I am trying to extract these values of T and Y, and use the vector T to find values for L and then plot these against Te" ?

[Kirsty James]

I'm taking the temperature function, Te, from the original daisyworld function. The vectors T and Y are calculated in ode45 and i want to use those values to calculate the Te at each point and plot it against luminosity

[Jan]

Remark: 5.67*10^(-8) requires one multiplication and an expensive power function. 5.67e-8 is parsed once during the reading of the M-file and in consequence much more efficient - and nicer.

Answer 1

[ChristianW]

There are several ways. I would make for the Plantary Temperature (Te) a new function, and vectorize Albedo (A). Like this:

function test
clear % Remove stored variables daisyworldode45 
tspan = [0 1000]; %solve for values of t
L = @(t) t/range(tspan)*2 + 0.4; %luminosity function want to keep within a range of 0 and 2 so change depending on tspan 
IC = [0.2 0.3]; %initial conditions of daisy percentage, black then white
[T,Y] = ode45(@(t,y)daisyworld1(t,y,L),tspan,IC); % solves equation, need capital T and Y 

LT = L(T);
Te = Temp(LT,Y);
subplot(211), plot(T,Y), xlabel('T, Time'), legend('Y_1','Y_2')
subplot(212), plot(LT,Te)
xlabel('L(T), Luminosity'), ylabel('Te, Plantary Temperature')

function dydt = daisyworld1(t,y,L)
dydt = [0;0];
Lt = L(t); % calculate L at time t
Te = Temp(Lt,y);
B = 1-0.003265*((22.5-Te).^2); % beta value
g = 0.3; % death term
dydt(1) = y(1)*((  (1-y(1)-y(2))  *B)-g); % black daisy formula
dydt(2) = y(2)*((  (1-y(1)-y(2))  *B)-g); % white daisy formula

function Te = Temp(Lt,y)
if isvector(y), y = y(:)'; end % make row vector
A = ((1-y(:,1)-y(:,2))*0.5)+(y(:,1)*0.25)+(y(:,2)*0.75); % albedo
S = 917; % constant solar energy
Z = 5.67*10^(-8); % Stefan-Boltzmann constant;
Te = ((((S*Lt)/Z).*(1-A)).^(1/4))-273; % plantary temperature

[Kirsty James]

Thank you so much this has really been a great help, thank you for taking the time to help me.