orbit equation with ode45
6 views (last 30 days)
Show older comments
Hi, I am trying to solve an equation with ode45 but it does not actually give me what it's expect to. The equation is the classical mechanics one for orbits. a(phi)'' = a(phi) + mi*gamma/l^2 , a(phi) = 1/r(phi)
I have included F = -gm1m2/r^2 already. I guess my initial values might be wrong but I do not see another alternative. ( I dont want to use the analitical solution for this). I put the code down below the first part is just a adjustment of scales plus variables. Could anyone help me see what's wrong? What I get is a vector y with y(1) = -y(2) which doesnt make sense.
Since a = 1/r and I get y1 and y2 ... and I have used y1=da/dphi and y2=a then, what I got is y1 and y2 and it's simple to convert them to give me radius and velocity but they dont return reasonable values.
function SunsOrbit
% Simple 2 body problem with Jupiter and the Sun
% CIRCULAR ORBIT AND JUPITER's VELOCITY CONSTANT WERE SUPPOSED
%
%Scales
Scale.d=7.785e11;
Scale.M=1.989e30;
Scale.mi=Scale.M;
Scale.V=1e3;
Scale.G=6.67384e-11;
Scale.gamma=Scale.G*Scale.M^2;
Scale.P=Scale.M*Scale.V;
Scale.l=Scale.d*Scale.P;
%Setting up Constants
K.d=7.785e11/Scale.d ; %m distance Jupiter-Sun
K.Mj=1.898e27/Scale.M; %kg
K.Msun= 1.989e30/Scale.M; %kg
K.mi=K.Mj*K.Msun/(K.Mj+K.Msun);
K.Vj=1.707e4/Scale.V; %m/s Jupiter's velocity
K.Pj=K.Mj*K.Vj; % Jupiter's linear momentum
K.l=K.d*K.Pj; % Angular Momentum Modulus
K.G=6.67384e-11/Scale.G; % m3 kg-1 s-2
K.gamma = K.G*K.Mj*K.Msun;
%
Inits=[1e-5,1/K.d];
ThetaRange=linspace(-pi,pi,500);
[r,y]=ode45(@orbit,ThetaRange,Inits,[],K); % y(1) is 1/r(theta)
r=1./y(:,2);
Vj=-(r.^2).*y(:,1);
Ueff=-K.gamma./r + (K.l.^2)./2*K.mi.*r.^2;
[u,i]=pol2cart(Theta,r);
plot(r*Scale.d,Vj*Scale.V,'*')
axis image;
% Orbit's equation
function da=orbit(Theta,y,K)
da=[-y(2)+K.mi*K.gamma/K.l^2;-y(1)]; % u'' = -u + mi*gamma/l
%
% TRY TO SET IT ALL UP SO THAT THE
% ENERGY IS CONSTANT ALSO USE ENERGY AS
% A FUNCTION OF POTENTIAL ENERGY AND
% KINETIC ENERGY
0 Comments
Answers (1)
Torsten
on 30 Oct 2014
y(2) contains the solution of the ODE
-u'' = -u + K.mi*K.gamma/K.l^2
I guess you either want to solve
u'' = -u + mi*gamma/l
as written as a comment behind your orbit-function, or
u'' = u + mi*gamma/l
as written in your introductionary text.
Best wishes
Torsten.
0 Comments
See Also
Categories
Find more on Gravitation, Cosmology & Astrophysics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!