Plotting period using values from a for loop

I have this code I am working on that is plotting the temperature of a planet that is orbit a star but there is a red dwarf star on an elliptical orbit around the main star and I wanted to calculate the temperature of the planet based its distance from the two stars. That part I have done, next part I want to do is plot the change of temperature against the orbital period of the red dwarf.

Some info on the Red Dwarf, its orbital period is 1896.59 days and its distance from the planet varies from 3.75 AU to 0.25 AU. The idea for the plot is to say that at day 1 the red dwarf starts its orbit at 3.75 AU, then takes half of its orbit (948.3 days) to get to the distance of 0.25 AU then the rest of its orbital period to return back to 3.75 AU. The problem I was running into when trying to plot this was that it kept plotting the whole range of values for the planets temperature every day. What I want is to say at day 1 when it is at 3.75 AU the planets temperature is this, and then do the same for every day during the red dwarf's orbit.

Any help would be great, I have attached my code below and an output image of what I am getting currently.

%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
G=6.67408e-11;
MS=2e+30;
MP=5.972e+24;
MRD=1.2e+30;

%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));

%Red Dwarf's Period
%Using Kepler's 3rd law
P=a.^1.5;                   %Gives answer in years
PRD=P*(365)

vTP=[];
for roRD=3.75:-0.01:0.25
    %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
    vTP(end+1) = TP;
    
    subplot(2,1,1);
    plot(roRD,TP,'.r')
    hold on
    
    xlim([0 4])
    ylim([-58 -44])
    xlabel('Orbit(AU)') 
    ylabel('Temperature of Planet(C)') 
    
 
end
    
   
for P=1:+5:1896
    
    subplot(2,1,2);
    plot(P,vTP,'.b'); 
    hold on 
    xlim([0 1896])
    ylim ([-58 -45])
end

Answer 1

[James Browne]

Hello,

I did some modifications to your code and I think I got it most of the way there but it is hard to finish because I am not quite sure how to structure the data for subplot 2 since I do not know how the variables relate to each other.

One thing I did was vectorize the results of your calculations and run the plots only once, that sped things up and seemed to clean up the second plot. Keep in mind that when you have a plot command with "hold on" in a loop, not only does MATLAB call the plot function on every loop cycle, which can cause high processing times, it can also be very difficult to troubleshoot when you don't get the results you expect.

I think I got you most of the way there but I was not sure how to structure the x-axis data for the second subplot so I just created a generic index vector so I could look at the shape of the plot and see if it was still doing crazy things.

Hope this at least helps:

%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
G=6.67408e-11;
MS=2e+30;
MP=5.972e+24;
MRD=1.2e+30;
%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));
%Red Dwarf's Period
%Using Kepler's 3rd law
P=a.^1.5;                   %Gives answer in years
PRD=P*(365)
vTP=[];

roRD = 3.75:-0.01:0.25;

for i = 1:length(roRD)
    %Flux
    fluxRD=(Lrd/(4*pi*(roRD(i)*1.496e11)^2));        %flux of red dwarf
    fluxAverage=(fluxS+fluxRD);
    
    %Temp of planet
    Tplanet=((fluxAverage*(1-al))/(4*stefconst))^0.25;
    TP(i)= Tplanet-273.15;
    vTP(i) = TP(i);
end

subplot(2,1,1);
plot(roRD,TP,'.r')

xlim([0 4])
ylim([-58 -44])
xlabel('Orbit(AU)')
ylabel('Temperature of Planet(C)')

plotVectorX = 1:length(vTP);

subplot(2,1,2);
plot(plotVectorX,vTP,'.b');
xlim([0 1896])
ylim ([-58 -45])

[Adam Kelly]

This is the output image and it is on the right tracks of what I need, so for the x axis I need to plot for every day so a step of +1 and the peak of the temperature should be at half the orbital period which is 948 days since that is when the red dwarf is closests and then it should decrease back down to normal as it goes from 0.25 AU back to 3.75 AU. So the desired output I need is a realtive straight line till that peak at 948 days (on the x-axis) and then it falling back down to that straight line till 1896 days.