For loop running n times wont work

%% earth emission increase 
for n = .2:2
    lwabsorbed_atm = n/2;
    [Tatm(n),Tsurf(n)]=modelwithatm(lwabsorbed_atm);
end
plot(.2:.2:2, Tatm)
hold on 
plot(.2:.2:2, Tsurf)
hold off 
legend('ta','ts')
axis([0 2 200 400])

So i have a function saved in another script called modelwithatm and lwabsorbed_atm is a variable in that function (longwave radiation absorbed in the atmospehere) I commented out the variable definition in the function. I have been able to run this seemlessly for other variables but something tells me there is something wrong with the way I wrote the for loop and suggest what to plot. Tsurf and tatm are for temperatures.

[KSSV]

What error you are getting?

[Theodore Anderson]

Array indices must be positive integers or logical values.

Error in modelingexperimentalchanges (line 95)

[Tatm(n),Tsurf(n)]=modelwithatm(lwabsorbed_atm);

[Theodore Anderson]

function [Tatm, Tsurf] = modelwithatm(lwabsorbed_atm)
%variables 
albedo =.31;
Ks = 1361; %W m-2 
Re = 6173; %km 
Sigma = 5.67E-8;%W m-2 k-4 
Epsilon = 1; 
% with atm variables
atm_albedo = .35; 
surface_albedo = .2; 
swradiation_sun = .9; 
lwradiation_sun = .1; 
%lwabsorbed_atm =.95
% atm absorbs .95 ALL lw and .02 ALL sw multiply by these values
Epsilon_atm = .9; 
%symbols 
syms Tsurf
syms Tatm 
%equations
% absorbed by atm = 2swvectors sun/earthreflection, 2lw vectors sun/earth
% absorbed by surface sw sun, lw sun, lw atm 
%sw 
a = Ks*.9*pi*Re^2*0.02; % sun to atm
d = Ks*.9*pi*Re^2*(1-.02 -.35)*.2*.02; %sun to surface to atm
c = Ks*.9*pi*Re^2*(1-.02 -.35)*(1-surface_albedo); %sw sun to surface 
%lw
q =Ks*.1*pi*Re^2*lwabsorbed_atm; %Lw sun to atm 
f = Ks*.1*pi*Re^2*(1-lwabsorbed_atm); % lw sun to surface 
y = 4*pi*Re^2*Epsilon*Sigma*Tsurf^4*lwabsorbed_atm; %lw earth to atm % whats the vibe?
z = 4*pi*Re^2*Epsilon_atm*Sigma*Tatm^4; %lw from atm to surface
p = 4*pi*Re^2*Epsilon*Sigma*Tsurf^4*lwabsorbed_atm;%lw to atm % whats the vibe?
%what could change? what wont? why would they change?
%emission vectors
%lw emitted from surface = p 
%lw emitted from atm or 2*z = atm emission 
%sw emitted from atm = sun reflected to space = a (.98) rather than .02 &
%albedo .35 so basically c 

%absorbption & emission equations pieces
Energy_absorbed_atm = a + d + q + p; %swsuntoatm + swsuntoearthtoatm + lwsuntoatm + lwearthtoatm
Energy_absorbed_surf = c + f + z; %swsuntoearth + lwsuntoearth + lwatmtoearth
Energy_emitted_atm = 2*z; %lwatmbothvectors
Energy_emitted_surface =  p; %lwfromearth
%packed equations solving for temperatures
equationa = (Energy_absorbed_atm  == Energy_emitted_atm); 
equations = (Energy_absorbed_surf == Energy_emitted_surface);
[tempa, temps] = solve([equationa, equations],[Tatm,Tsurf]);
%equation solutions for temperature in absolute terms
Tsurf = abs(real(double(temps(1))));
Tatm = abs(real(double(tempa(1))));
end

Answer 1

[Alan Stevens]

You have

for n = .2:2
...

But then call

[Tatm(n),Tsurf(n)] ...

However, the indices to Tatm(n) etc must be integers, not 0.2.

Change to something like

for n = 1:2  
    i = 0.2*n;  % or whatever
    lwabsorbed_atm = i/2;
    [Tatm(n),Tsurf(n)]=modelwithatm(lwabsorbed_atm);
end

Possibly you want to go from 0 to 2 in steps of 0.2 in which case

for n = 1:11  
    i = 0.2*(n-1);  % or whatever  
    lwabsorbed_atm = i/2;
    [Tatm(n),Tsurf(n)]=modelwithatm(lwabsorbed_atm);  % the indices must be >= 1
end

[Theodore Anderson]

Thank you so much!

[Theodore Anderson]

but i have one more question, when I change the n = to positive integers i get a strange result for my surface temperature. Theoretically if I raisie the longwave absorbtion in the atmosphere it will raise the temperature of the surface as you can see it is a multiplier for several of the equations but when I ran the code i got this result:

[Alan Stevens]

Show the modified code you are using to call modelwithatm (including the for loop and lwabsorbed_atm specification).

[Alan Stevens]

I get a somewhat different result having manipulated your equations, however, Tsurf still decreases with increasing lwabsorbed_atm.

lwabsorbed_atm = zeros(1,10);
Tatm = zeros(1,10);
Tsurf = zeros(1,10);

for n = 1:10 
    i = 0.2*n; 
    lwabsorbed_atm(n) = i/2;
    [Tatm(n),Tsurf(n)]=modelwithatm(lwabsorbed_atm(n));  
end

plot(2*lwabsorbed_atm,Tatm,'b',2*lwabsorbed_atm,Tsurf,'r'),grid
legend('Tatm','Tsurf')

function [Tatm, Tsurf] = modelwithatm(lwabsorbed_atm)
%variables 
albedo =.31;
Ks = 1361; %W m-2 
Re = 6173; %km 
Sigma = 5.67E-8;%W m-2 k-4 
Epsilon = 1; 
% with atm variables
atm_albedo = .35; 
surface_albedo = .2; 
swradiation_sun = .9; 
lwradiation_sun = .1; 
%lwabsorbed_atm =.95
% atm absorbs .95 ALL lw and .02 ALL sw multiply by these values
Epsilon_atm = .9; 
%symbols 
%equations
% absorbed by atm = 2swvectors sun/earthreflection, 2lw vectors sun/earth
% absorbed by surface sw sun, lw sun, lw atm 
%sw 
a = Ks*.9*pi*Re^2*0.02; % sun to atm
d = Ks*.9*pi*Re^2*(1-.02 -.35)*.2*.02; %sun to surface to atm
c = Ks*.9*pi*Re^2*(1-.02 -.35)*(1-surface_albedo); %sw sun to surface 
%lw
q =Ks*.1*pi*Re^2*lwabsorbed_atm; %Lw sun to atm 
f = Ks*.1*pi*Re^2*(1-lwabsorbed_atm); % lw sun to surface 
% y = 4*pi*Re^2*Epsilon*Sigma*Tsurf^4*lwabsorbed_atm; %lw earth to atm % whats the vibe?
% z = 4*pi*Re^2*Epsilon_atm*Sigma*Tatm^4; %lw from atm to surface
%p = 4*pi*Re^2*Epsilon*Sigma*Tsurf^4*lwabsorbed_atm;%lw to atm % whats the vibe?
%what could change? what wont? why would they change?
%emission vectors
%lw emitted from surface = p 
%lw emitted from atm or 2*z = atm emission 
%sw emitted from atm = sun reflected to space = a (.98) rather than .02 &
%albedo .35 so basically c 

%absorbption & emission equations pieces
%Energy_absorbed_atm = a + d + q + p; %swsuntoatm + swsuntoearthtoatm + lwsuntoatm + lwearthtoatm
% Energy_absorbed_surf = c + f + z; %swsuntoearth + lwsuntoearth + lwatmtoearth
% Energy_emitted_atm = 2*z; %lwatmbothvectors
% Energy_emitted_surface =  p; %lwfromearth

% Energy_absorbed_surf = Energy_emitted_surface
% c + f + z = p;               ...(1)
% Energy_absorbed_atm = Energy_emitted_atm
% a + d + q + p = 2*z 
% 2*z - a - d - q = p          ...(2)
% (2) - (1):  z-a-d-q-c-f = 0
%             z = a+d+q+c+f    ...(3)

z = a + d + q + c + f;     %  (3)
Tatm = (z/(4*pi*Re^2*Epsilon_atm*Sigma))^(1/4);

p = c + f + z;             % (1)
Tsurf = (p/(4*pi*Re^2*Epsilon*Sigma*lwabsorbed_atm))^(1/4);

end

This results in