0 个投票

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

请先登录，再回答此问题。

 采纳的回答

Ameer Hamza
Ameer Hamza 2020-5-19

0 个投票

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

6 个评论
显示 4更早的评论 隐藏 4更早的评论

christian fares
christian fares 2020-5-19
Thank you so much it worked!
Ameer Hamza
Ameer Hamza 2020-5-19
I am glad to be of help.
christian fares
christian fares 2020-5-19
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
Ameer Hamza 2020-5-19
Yes, from the code in your question. Both appear to be the same. You can say that it for a general case.
christian fares
christian fares 2020-5-19
okay thanks again! you saved my life it's a project for the matlab course.
Have a great day :)
Ameer Hamza
Ameer Hamza 2020-5-19
I am glad to be of help.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Earth and Planetary Science 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by