Plot of temperature changing with orbit

Question

Adam Kelly 2020-4-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/515928-plot-of-temperature-changing-with-orbit

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