The plot is empty

The purpous of this code is for gravity attraction. I found a code that was perfectly working but i wanted to add more planets to it, so i added mercury planet.The plot was working before adding a new planet but now it's not. Actually, the code is running and there's no errors

Gravity:

%% Inputs

stepsize=1; %step size in days, should not exceed 1

finaltime=100; %days from t0

%% System Conditions

global sunm G

sunm=1.9891e30; %kg

G=1.4879*10^(-34); %AU^3/(kg*Day^2)

%% Sun Initial Conditions

x=0;

y=0;

z=0;

Vx=0;

Vy=0;

Vz=0;

SunInitial=[x,y,z,Vx,Vy,Vz]; %summary of initial conditions of the sun

%Note: not used in rest of code.

%% Earth Initial Conditions

x=-7.829690890802676E-01; %AU

y=6.009223815219861E-01; %AU

z=-1.377986550047871E-05; %AU

Vx=-1.075646922487205E-02;%AU/day

Vy=-1.370584048238274E-02;%AU/day

Vz=3.974096024543801E-08; %AU/day

EarthInitial=[x,y,z,Vx,Vy,Vz]; %summary of initial conditions of the earth

%% New planet

Mx=-2.829690890802676E-01; %AU

My=8.009223815219861E-01; %AU

Mz=-4.377986550047871E-05; %AU

MVx=-0.075646922487205E-02;%AU/day

MVy=-2.370584048238274E-02;%AU/day

MVz=5.974096024543801E-08; %AU/day

MercuryInitial=[Mx,My,Mz,MVx,MVy,MVz];

%% Call the Function

%outpu=call the function, time to run, initial conditions

[attime,finalposition]=ode45(@Propagator,[0:stepsize:finaltime],[EarthInitial,MercuryInitial ]);

%% Plot the Result

hold on

plot3(finalposition(:,4),finalposition(:,5),finalposition(:,6),'g'); %plot of position of earth in green

plot3(finalposition(:,10),finalposition(:,11),finalposition(:,12),'g');

xlabel('Distance (AU)');ylabel('Distance (AU)');zlabel('Distance (AU)'); %label the axis

legend('Earth'); %tell me what the green line is

axis equal

hold off

%fix the axis so it has the same scale in x, y and z

%note: sun is at 0,0,0

Propagator function:

function [out]=Propagator(t,PlanetPositions);

global sunm G

%state space representation idot=A1*PlanetPositions+A2(PlanetPositions)

%where idot is the derivative of the PlanetPositions

%this represents the linear part in A1 added to the non-linear part in A2

% |Vx| |0 0 0 1 0 0| |x | | 0 |

% |Vy| |0 0 0 0 1 0| |y | | 0 |

% |Vz| = |0 0 0 0 0 1| * |z | + | 0 |

% |Ax| |0 0 0 0 0 0| |Vx| | -G*sunm/R^2*|Rx| |

% |Ay| |0 0 0 0 0 0| |Vy| | -G*sunm/R^2*|Ry| |

% |Az| |0 0 0 0 0 0| |Vz| | -G*sunm/R^2*|Rz| |

% idot = A1 * PlanetPositions + A2(PlanetPositions)

% where |Rx| is the unit vector of R in the x direction

% where R=sqrt(x^2+y^2+z^2);

PlanetPositions= zeros(6,2);

dx=zeros(6,1); %initialize the xdot matrix

A1=zeros(6,6);

A2=zeros(6,1);

A1(1,4)=1;

A1(2,5)=1;

A1(3,6)=1;

%% Earth

x=PlanetPositions(1); %x position

y=PlanetPositions(2); %y position

z=PlanetPositions(3); %z position

R = norm([x,y,z]); %find the radius from the earth to the sun

unitvector=[x,y,z]/R;

% do these calculations now so we only have to calculate it once

A2(4)=-G*sunm/R^2*unitvector(1); %acceleration in the x

A2(5)=-G*sunm/R^2*unitvector(2); %acceleration in the y

A2(6)=-G*sunm/R^2*unitvector(3); %acceleration in the z

TotAccelEarth=A1*PlanetPositions+A2; %sum up the linear + nonlinear accelerations

out=TotAccelEarth; %output the accelerations

%% Mercury

Mx=PlanetPositions(7); %x position

My=PlanetPositions(8); %y position

Mz=PlanetPositions(9); %z position

R = norm([Mx,My,Mz]); %find the radius from mercury to the sun

unitvector=[Mx,My,Mz]/R;

A2(4)=-G*sunm/R^2*unitvector(1); %acceleration in the x

A2(5)=-G*sunm/R^2*unitvector(2); %acceleration in the y

A2(6)=-G*sunm/R^2*unitvector(3); %acceleration in the z

TotAccelMercury=A1*PlanetPositions+A2; %sum up the linear + nonlinear accelerations

out=TotAccelMercury; %output the accelerations

out = out(:);

end

Answer 1

[Ameer Hamza]

I am not sure why you have written the line

PlanetPositions = zeros(6,2);

so ode45 is not able to do anything. Also, it seems that you are trying to solve two independent systems of equations. Therefore, I suggest two calls to ode45. Run this code

%% Inputs
stepsize=1;   %step size in days, should not exceed 1
finaltime=100; %days from t0
%% System Conditions
global sunm G
sunm=1.9891e30; %kg
G=1.4879*10^(-34); %AU^3/(kg*Day^2)
%% Sun Initial Conditions
x=0;
y=0;
z=0;
Vx=0;
Vy=0;
Vz=0;
SunInitial=[x,y,z,Vx,Vy,Vz]; %summary of initial conditions of the sun
%Note: not used in rest of code.
%% Earth Initial Conditions
x=-7.829690890802676E-01; %AU
y=6.009223815219861E-01;  %AU
z=-1.377986550047871E-05; %AU
Vx=-1.075646922487205E-02;%AU/day
Vy=-1.370584048238274E-02;%AU/day
Vz=3.974096024543801E-08; %AU/day
EarthInitial=[x,y,z,Vx,Vy,Vz]; %summary of initial conditions of the earth
%% New planet 
Mx=-2.829690890802676E-01; %AU
My=8.009223815219861E-01;  %AU
Mz=-4.377986550047871E-05; %AU
MVx=-0.075646922487205E-02;%AU/day
MVy=-2.370584048238274E-02;%AU/day
MVz=5.974096024543801E-08; %AU/day
MercuryInitial=[Mx,My,Mz,MVx,MVy,MVz];
%% Call the Function
%outpu=call the function,      time to run,      initial conditions
[attime_earth,finalposition_earth]=ode45(@Propagator,0:stepsize:finaltime,EarthInitial);
[attime_mercu,finalposition_mercu]=ode45(@Propagator,0:stepsize:finaltime,MercuryInitial);
%% Plot the Result
hold on
plot3(finalposition_earth(:,4),finalposition_earth(:,5),finalposition_earth(:,6),'g'); %plot of position of earth in green
plot3(finalposition_mercu(:,4),finalposition_mercu(:,5),finalposition_mercu(:,6),'g--');
xlabel('Distance (AU)');ylabel('Distance (AU)');zlabel('Distance (AU)'); %label the axis
legend('Earth'); %tell me what the green line is
axis equal
hold off
%fix the axis so it has the same scale in x, y and z

function out = Propagator(t, PlanetPositions)
global sunm G
%state space representation idot=A1*PlanetPositions+A2(PlanetPositions) 
%where idot is the derivative of the PlanetPositions
%this represents the linear part in A1 added to the non-linear part in A2
% |Vx|   |0 0 0 1 0 0|   |x |   | 0                |
% |Vy|   |0 0 0 0 1 0|   |y |   | 0                |
% |Vz| = |0 0 0 0 0 1| * |z | + | 0                |
% |Ax|   |0 0 0 0 0 0|   |Vx|   | -G*sunm/R^2*|Rx| |
% |Ay|   |0 0 0 0 0 0|   |Vy|   | -G*sunm/R^2*|Ry| |
% |Az|   |0 0 0 0 0 0|   |Vz|   | -G*sunm/R^2*|Rz| |
% idot =      A1       * PlanetPositions   +     A2(PlanetPositions)
% where |Rx| is the unit vector of R in the x direction
% where R=sqrt(x^2+y^2+z^2);
% PlanetPositions= zeros(6,2);
dx=zeros(6,1); %initialize the xdot matrix
A1=zeros(6,6);
A2=zeros(6,1);
A1(1,4)=1;
A1(2,5)=1;
A1(3,6)=1;
%% Earth
x=PlanetPositions(1); %x position
y=PlanetPositions(2); %y position
z=PlanetPositions(3); %z position
R = norm([x,y,z]); %find the radius from the earth to the sun 
unitvector=[x,y,z]/R;
% do these calculations now so we only have to calculate it once
A2(4)=-G*sunm/R^2*unitvector(1); %acceleration in the x
A2(5)=-G*sunm/R^2*unitvector(2); %acceleration in the y
A2(6)=-G*sunm/R^2*unitvector(3); %acceleration in the z
TotAccelEarth=A1*PlanetPositions+A2; %sum up the linear + nonlinear accelerations
out=TotAccelEarth(:); %output the accelerations
end

[christian fares]

Thank you so much it worked!

[Ameer Hamza]

I am glad to be of help.

[christian fares]

But i have a question , the propagator function is for any planet and we call it to make it work for any planet? because i need to remove %% earth in the propagator function since its for any planet

[Ameer Hamza]

Yes, from the code in your question. Both appear to be the same. You can say that it for a general case.

[christian fares]

okay thanks again! you saved my life it's a project for the matlab course.

Have a great day :)

[Ameer Hamza]

I am glad to be of help.