Coding Crowell's Method for Orbiting Bodies
2 views (last 30 days)
Show older comments
I'm trying to make a function that takes some initial conditions and models of planetary bodies and it spits out a plot of what happens over time.
To do this I'm using the following differential equation,
Which is the formula given on Wikipedia (Orbit modeling - Wikipedia). The are vectors that represent the location of the bodies, and n is the number of these bodies, G is the gravitational constant 6.672e-11 . Now I'm using the substitutions, , and , to turn this into into a system of first order ODEs. First, before defining the ODE function I'll explain how the inputs work, I use this function here (the plotting part is unfinished, I'm just testing if the caluclations actually work),
function [T,V] = OrbModel(m,x0,v0,t,evalnum)
% INPUT EXPLANATIONS
% OUTPUT EXPLANATIONS
dim = length(x0(:,1));
p = length(m);
T = linspace(t(1),t(2),evalnum);
y0 = [x0(:);v0(:)];
[~,V]=ode45(@(t,y)ODEFunc(y,m),T,y0);
switch dim
case 2
plotres=zeros(1,p);
hold on
for i=1:p:dim*p
plot(V(:,i),V(:,i+1));
end
hold off
case 3
end
m is a vector that contains the masses of each planet. is a n by p matrix, where n is the dimension I'm working in, e.g. in 2 dimensions their positions are described by 2 dimensional vectors and p are the number of planets. Therefore, the first column of is the initial coordinates of the first planet, etc. is the exact same as but initial velocities instead. Now ideally I'd want to pass the matrix into the ode45( . ) function but I can't so I combined and into a single vector, in the code above. Now ODEFunc( . ) is the following,
function dy=ODEFunc(y,m)
G = 6.672e-11;
p = length(m);
dim = length(y)/(2*p);
dimc=dim-1;
dy = zeros(dim*p*2,1);
for i=1:p
dy(i*dim-dimc:i*dim) = y(i*dim-dimc+dim*p:i*dim+dim*p);
for k=[1:i-1, i+1:p]
dy(i*dim-dimc+p*dim:i*dim+p*dim) = dy(i*dim-dimc+p*dim:i*dim+p*dim)...
+G*m(k)*(y(k*dim-dimc:k*dim)-y(i*dim:i*dim))/norm(y(k*dim-dimc:k*dim)-y(i*dim-dimc:i*dim))^3;
end
end
end
which (hopefully) resembles the function used in Crowell's Method with the substitutions above. Now with these two functions I'm testing it by plotting an Earth-Sun system, the inputs I'm using are m=[1.989e30, 5.972e24], x0=[0 0; 1.4756e11 0], v0=[0 0; 0 30000], t and evalnum are unimportant. These inputs I got from a quick search on google, more specifically I'm assuming the sun to be at the center and the Earth , and the velocity of the sun to be 0 and the earth has an upwards velocity of m/s. With these inputs I get the following plot,
The orange line is Earth, and it starts at the bottom right. Now it seems to have relatively appropriate behavior of planetary bodies, arcing back towards the sun at (0,0), but then the earth just leaves, off to the right to infinity (at least it seems like that). Finally, my question is, since it seems to exhibit decent behavior it's really hard to tell if the code itself is incorrect or just my initial conditions are bad. As such, I'm wondering which one it is. If it's the initial conditions how would I go about finding the correct ones?
0 Comments
Accepted Answer
James Tursa
on 29 Nov 2020
Edited: James Tursa
on 29 Nov 2020
I don't have the patience to go through all of your dim and dimc logic. It would seem to be much easier to code and debug if you reshaped things first. E.g.,
function dy=ODEFunc(y,m)
G = 6.672e-11;
p = length(m);
y = reshape(y,[],2*p); % columns of positions and velocities
x = y(:,1:p); % pick off positions
v = y(:,p+1:end); % pick off velocities
dv = zeros(size(v)); % allocate velocity derivative matrix
for i=1:p
for k=[1:i-1, i+1:p]
dv(:,i) = dv(:,i) + _____; % insert code here to calculate acceleration of i'th planet
end
end
dy = [v(:);dv(:)];
end
You just need to fill in the blank, using x(:,i) and x(:,k) as the positions of the i'th and k'th planets respectively.
2 Comments
James Tursa
on 29 Nov 2020
Edited: James Tursa
on 29 Nov 2020
If your problem was solved then it wasn't a waste!
Note that you really don't need dim. You can just do:
dy = zeros(size(y));
More Answers (0)
See Also
Categories
Find more on Numerical Integration and Differential Equations 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!