Plot of temperature changing with orbit

2 次查看(过去 30 天)
Hey, I am stuck with this problem I having with my code.
So I have gotten a bunch of values for temperature of a planet that changes when a red dwarf that is on an elliptical orbit gets closer to the planet. Now I want to plot how the temperature of the planet changes during the Red Dwarf's orbit. The idea is to have the temperature of the planet on the y-axis and the Period of the Red Dwarf on the x-axis.
The problem I am having is when plot it, it is plotting all the values for the temperature of the planet for ever point of the Red Dwarf's period I am plotting which is every 5th day of its period which is 1896 and this isn't what I want, I want to show that for the start of the Red Dwarf's orbit the planet's temperature doesn't change until it gets close to the planet so i should get this spike of temperature and then when the Red Dwarf moves further away it falls back down to the normal temperature of the planet.
Below I have attached the ouput image and the code I have so far. Any help with this would be great!!
%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;
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=0.25:+0.01:3.75
%Flux
fluxRD=(Lrd./(4*pi*(roRD*1.496e11)^2)); %flux of red dwarf
fluxAverage=(fluxS+fluxRD);
%u=G*MRD;
%v=sqrt(u*(2/(4.5*1.496e11))-(1/(a*1.496e11)))
%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 5])
ylim([-58 -45])
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 个)

类别

Help CenterFile Exchange 中查找有关 Earth and Planetary Science 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by