Using values from a for loop in another for loop

I am currently stuck with a problem with my code. The idea is using a range of values for the temperature of a planet due to the flux from two stars, that is calculated in a for loop and then using those values to plot them against the orbital period of one of the stars that has a longer peroid.

The problem I am having is when I try to plot the orbital peroid against time, it only uses the last value of temperature giving me a straight line. (see image below)

Any help with this would be appreciated!

I have attached the code and my image below;

%Written by Adam Kelly
clear all;
clc;

%constants
r=1.5e+10;
rs=6.261e+8;        %Radius of Star
ro=1.49e11;         %Sun orbital radius (1 AU)
rRD=1.25e8 ;        %Radius of Red Dwarf
stefconst=5.67e-8;  %stefan-boltzmann constant
T=5200.2;           %Temp of Sun
TRD=3773.15;        %Temp of Red Dwarf
al=0.3;             %albedo
a=3;                %Red Dwarf's semi major axis

%luminosity of star
Ls=4*pi*rs^2*stefconst*T^4;
Lrd=4*pi*rRD^2*stefconst*TRD^4;

fluxS=(Ls./(4*pi*ro^2));

for roRD=0.25:+0.01:4.5
    %Flux
    fluxRD=(Lrd./(4*pi*(roRD*1.496e11)^2));        %flux of red dwarf 
    fluxAverage=(fluxS+fluxRD);
   
    %Temp of planet
    Tplanet=((fluxAverage*(1-al))./(4*stefconst))^0.25;
    TP=Tplanet-273.15;
    disp(TP)
    
    subplot(2,1,1);
    plot(roRD,TP,'.r')
    hold on
    xlim([0 5])
    ylim([-58 -45])
    xlabel('Orbit(AU)') 
    ylabel('Temperature of Planet(C)') 
end

%Red Dwarf's Period
%Using Kepler's 3rd law

P=a.^1.5;                   %Gives answer in years
PRD=P*(365);
disp(PRD)

for PRD=0:+100:2000
    P_RD=PRD;
    subplot(2,1,2);
    plot(P_RD,TP,'.b'); hold on 
    xlim([0 2e+3])
    ylim([-58 -45])
    xlabel('Time(days)') 
    ylabel('Temperature of Planet(C)')
end

[Arthur Roué]

You need to store the values of TP in an array, because each time it's evaluated in the first loop the value is overwritten.

Without preallocation, you can write :

vTP = [];
% First loop
for roRD=0.25:+0.01:4.5
    ...
    TP=Tplanet-273.15;
    vTP(end+1) = TP;
    ...
end

Then use vTP in the second loop.

[Adam Kelly]

This is the output image I get, it is plotting every temp, where I need singular points of temperture for each day. So at zero since the orbit starts furthest away it should be at the lowest temp -55 and should peak at -44 about half way in time and then decrease back down to -55.

Do you have an idea how to achieve this? If not don't worry, you've helped a lot!

Thank you!

Answer 1

[Divya Yerraguntla]

Hi Adam,

After making the modification suggested by Arthur you are plotting the data for each iteration in your 2nd for loop. So try removing the second for loop and just use the below code in place of it.

PRD = 0:100:42500; % To have PRD and vTP of the same size I changed the upper limit of PRD
subplot(2,1,2);
plot(PRD,vTP,'.b');  
xlim([0 43e+3]);
ylim([-58 -45]);
xlabel('Time(days)'); 
ylabel('Temperature of Planet(C)');

This helps you plot the data at once and you wouldnt need hold on as well.

Hope it helps!