MATLAB Answers

Plotting period using values from a for loop

2 ビュー (過去 30 日間)
Adam Kelly
Adam Kelly 2020 年 4 月 10 日
コメント済み: Adam Kelly 2020 年 4 月 10 日
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

  0 件のコメント

サインインしてコメントする。

採用された回答

James Browne
James Browne 2020 年 4 月 10 日
編集済み: James Browne 2020 年 4 月 10 日
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])

  1 件のコメント

Adam Kelly
Adam Kelly 2020 年 4 月 10 日
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.

サインインしてコメントする。

その他の回答 (0 件)

タグ

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by