MATLAB Answers

Hello can someone show me a simple code to show the earths orbit around the sun using ordinary differential equations.

272 views (last 30 days)
Matthew Worker
Matthew Worker on 25 Jan 2018
Edited: James Tursa on 13 Mar 2018
function Earthsorbit
e=2;
y0=[1 2*pi 0 -pi];
tspan=[0, 10];
options = odeset('RelTol',1e-4,'AbsTol',1e-4);
[t,y]=ode45(@earthorbit,tspan,y0,options);
for i = 1:length (t)
fprintf ('%2i %7.5f %7.5f \n',i,t(i),y(i));
end
plot(t,y(:,1),'-');
title('The solution will give us the orbit of a body in two dimensions.');
xlabel('time t');
ylabel('solution y');
end
function f=earthorbit(t,y)
f=zeros(4,1);
GM=4*pi^2;
r=sqrt ( y(1)^2 + y(3)^2 ); %distance between objects
f(1) = y(2);
f(2) =- GM*y(1) / r^(3/2);
f(3) = y(4);
f(4) = - GM*y(3) / r^(3/2);
end
This isnt showing me the circle I desire.
  5 Comments

Sign in to comment.

Answers (2)

Álvaro Romero Calvo
Álvaro Romero Calvo on 25 Jan 2018
Hi!
If you want a simple 2 body approach in 2D without considering perturbations, take a look at this (parameters extracted from NASA):
function [r,v,T] = simpleorbit
% Executes the 2D Earth orbit integration around Sun with 2 bodies
% equations
% Orbit parameters:
e = 0.0167; % Orbit eccentricity
mu_S = 132712439935.5; % Sun parameter (Km^3/s^2)
R_S = 6.957*10^5; % Sun radius (Km)
a = 149.6*1e6; %Semimajor axis (Km)
theta0 = 0;
% Time
T = 2*pi*sqrt(a^3/mu_S);
tfin = T;
% IC:
rx0 = a*(1-e^2)/(1+e*cos(theta0));
vy0 = sqrt(2*mu_S/rx0-mu_S/a);
X0 = [rx0,0,0,0,vy0,0];
options = odeset('Reltol',1e-15,'AbsTol',1e-15);
% Integration:
[T, X] = ode113(@dynamics, [0,tfin], X0, options, mu_S);
r = [X(:,1),X(:,2),X(:,3)];
v = [X(:,4),X(:,5),X(:,6)];
% Additional results:
h = cross(r,v);
ev = cross(v(1,:),h(1,:))/mu_S-r(1,:)/norm(r(1,:));
% Draw planet
figure
[sx,sy,sz] = sphere(100);
surf(R_S*sx,R_S*sy,R_S*sz,'EdgeColor','none','FaceColor','y')
% Aspect
hold on
axis equal
grid on
axis([-2*a,2*a,-2*a,2*a,-2*R_S,2*R_S])
% Plot:
% quiver3(0,0,0,h(1,1),h(1,2),h(1,3))
quiver3(0,0,0,ev(1),ev(2),ev(3))
comet3(r(:,1),r(:,2),r(:,3))
quiver3(r(1:10:end,1),r(1:10:end,2),r(1:10:end,3),v(1:10:end,1),v(1:10:end,2),v(1:10:end,3))
end
function [dX] = dynamics(~,X,mu)
rx = X(1);
ry = X(2);
rz = X(3);
vx = X(4);
vy = X(5);
vz = X(6);
r = sqrt(rx.^2+ry.^2+rz.^2);
dX(1,1) = vx;
dX(2,1) = vy;
dX(3,1) = vz;
dX(4,1) = -mu*rx/r.^3;
dX(5,1) = -mu*ry/r.^3;
dX(6,1) = -mu*rz/r.^3;
end
For advanced orbit calculations, you will need to use a propagator (like Gauss equations).
Hope it helps!
  8 Comments

Sign in to comment.


James Tursa
James Tursa on 2 Feb 2018
Edited: James Tursa on 2 Feb 2018
Your biggest problem is that you have a 4-element state vector defined with y0 (indicating only x-y plane motion), but your derivative function assumes a 6-element state vector for motion in 3D. That's a mismatch. Pick one or the other and be consistent. I would advise going with the 3D version ... you can always set up planar motion if you want.
Next problem is your initial state vector y0. Why would pi be involved as a state vector element? I might expect to see sin(pi) or cos(pi) here if you were setting things up to initially be somewhere on a circle, but not pi directly as you have.
You need to define a mu value for your gravity calculations and pass that into your derivative code. E.g., something like this:
mu = _____; % some appropriate value for the state vector units you are using
[t,y]=ode45(@(t,y)earthtwo(t,y,mu),tspan,y0,options);
and then change the signature of your derivative function to include mu:
function f= earthtwo(t,y,mu)
So, you need to fix up your initial state vector, define a value of mu, and then change your derivative code to use y and mu properly. E.g.,
rx = y(1);
ry = y(2);
rz = y(3);
vx = y(4);
vy = y(5);
vz = y(6);
You've been struggling with this for about a week now, so I went ahead and modified your code to run in 3D making the above corrections. I went ahead and picked the units to be meters and seconds since that is what the units of the sun gravity was that I got off the Wiki. But you can change the units to be something else and the program will still run properly (the derivative function doesn't care what the units are ... it just needs mu to be consistent units with y). Finally, the signature of the derivative function that ode45 calls must be (t,y). Since we also want mu to be passed in (I think this is nicer than hard-coding a mu value into your derivative routine), the clean way to do this is to create a function handle on the fly that has a (t,y) signature but calls your actual derivative function with any extra parameters that you want such as mu: @(t,y)earthtwo(t,y,mu)
To change the code to run satellites orbiting the Earth, all you need to do is change mu to Earth gravity and set up an appropriate initial state vector and time span.
function earthmain
mu = 1.32712440018e20; % m^3/s^2
r = 149597870700; % m (= 1 AU)
vcirc = sqrt(mu/r); % m/s (circular velocity for the given value of r)
y0 = [r; 0; 0; 0; vcirc; 0]; % m and m/s (start on x-axis with +y velocity direction)
tspan=[0, 365.25*86400]; % One year time span
options = odeset('RelTol',1e-4,'AbsTol',1e-4);
[t,y]=ode45(@(t,y)earthtwo(t,y,mu),tspan,y0,options);
% for i = 1:length (t)
% fprintf ('%2i %7.5f %7.5f \n',i,t(i),y(i));
% end
figure;
plot(y(:,1),y(:,2),'-'); % Plot only the x-y plane since that is how we set things up
axis square
grid on
title('The solution will give us the orbit of a body in two dimensions.');
xlabel('x (m)');
ylabel('y (m)');
end
function f= earthtwo(t,y,mu)
rx = y(1);
ry = y(2);
rz = y(3);
vx = y(4);
vy = y(5);
vz = y(6);
r = sqrt(rx.^2+ry.^2+rz.^2);
f=zeros(6,1);
f(1) = vx;
f(2) = vy;
f(3) = vz;
f(4) = -mu*rx/r.^3;
f(5) = -mu*ry/r.^3;
f(6) = -mu*rz/r.^3;
end
  9 Comments
James Tursa
James Tursa on 13 Mar 2018
It was surprising to me how difficult it was to locate an online reference for the equations you are using. There are a lot of references for the orbital mechanics equations (technically, the integrals of the relative differential equations of motion), and lots of references for the inertial differential equations of motion, but I couldn't find any references to the relative equations of motion for the two-body problem ... until today. The reference is here:
What is the difference between the inertial equations of motion and the relative equations of motion? Simply put, the inertial equations of motion for the two body problem are for an inertial frame (i.e. Newtonian frame) of reference and comprise a 12th order system of ODE (3 position elements and 3 velocity elements for each body with respect to an inertial frame). The relative equations of motion for the two-body problem are where you subtract their inertial equations of motion from each other, resulting in the position and velocity of one body relative to the other. This reduces what started out as a 12th order ODE system to a 6th order ODE system. And it is the solution of this 6th order "relative" ODE system that is always presented in the literature (6 integrals called orbital mechanics equations, and 6 constants of integration called orbital elements). You can find this relative ODE system discussion in section 1.3 of the above link. In particular, you will note equation 1.33 (this is the ODE that you are solving in your MATLAB code):
(Equation 1.33) r2 dotdot - r1 dotdot = r dotdot = -G*(m1+m2) r / r^3
On the left of this equation you can see the position derivatives of the two bodies subtracted which is the big clue that these are the relative equations of motion and not the inertial equations of motion. Also, you see the term with the masses added as G*(m1+m2) instead of the masses multiplied as G*m1*m2. This is the other clue that these are the relative equations of motion instead of the inertial equations of motion. It is the vectors r and r dot that you are solving for.
This relates to your code as follows:
The r vector in the relative ODE equations (1.33) is the first three elements y(1:3) of your y state vector, and the velocity r dot is the last three elements y(4:6) of your y state vector. The mu in your code is simply the G*(m1+m2) term in the ODE (Equation 1.34). From this, you should easily be able to correlate the relative ODE given in the 1.33 line with your earthtwo( ) function code.
Side Note: For satellite work where the mass of the planet far exceeds the mass of the satellite, only the mu of the planet is needed. But for planets orbiting the sun the planet mass, while much smaller than the sun, probably shouldn't be considered negligible. In that case, it would be technically more appropriate to combine the two mu's. E.g., in your case you would use a mu = mu_sun + mu_earth in your equations. E.g.,
mu_sun = 1.32712440018e20; % m^3/s^2
mu_earth = 3.986004418e14; % m^3/s^2
mu = mu_sun + mu_earth;

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by