I'm trying to solve a 3 body problem consisting of the earth moon and sun. I'm wanting to show the earth orbiting the sun and the moon orbiting the earth as a result of the gravitational forces.
Where Force due to gravity is given by:
$$\mathbf F_{ij}=\frac{Gm_i m_j}{\|\mathbf q_j-\mathbf q_i\|^2}\cdot\frac{(\mathbf q_j-\mathbf q_i)}{\|\mathbf q_j-\mathbf q_i\|}=\frac{Gm_i m_j(\mathbf q_j-\mathbf q_i)}{\|\mathbf q_j-\mathbf q_i\|^3}$$
And the acceleration of the bodies can be given as
$$\vec{a}_j=\sum_{i\ne j}^n G\frac{M_i}{\left|\vec{r}_i-\vec{r}_j\right|^3}\left(\vec{r}_i-\vec{r}_j\right)$$
INITIAL CONDITIONS
The initial conditions of the bodies are as follows:
Ecc = {0.01671022, 0.00549}(*eccentricity of earth and moon orbits*);
Inc = {0.00005, 5.145};
RAAN = {-11.260640, 0};
ArgPer = {102.94719, 0};
G = 0.00029589743849552926`;(*gravitational constant in AU*)
m = {1, 3.004*10^-6,
3.694*10^-8}; (*mass of sun, earth and moon in solar masses*)
μ = Table[
G*m[[i]], {i, 1, 3}];(*standard gravitational parameter of *)
(*μ={8.8878*10^-10, 1.093*10^-11};*)(*standard gravitational \
constant\[Rule] Sunm, earth, moon*)
ψ = 0;
a = {1.00000011, 0.99743};(*semi major axis in AU*);
r = a (1 - Ecc^2)/(1 +
Ecc Cos[ψ])(* disatances for sun to earth and sun to moon
in AU*)
rx = r Cos[ψ](* x component of distance for: sun to earth and \
sun to moon*);
ry = r Sin[ψ](* y component of distance forsun to earth and sun \
to moon*);
v = Table[
Sqrt[μ[[i]] (2/r[[i]] - 1/a[[i]])], {i,
2}](*earth velocity due to sun , moon velocity due to earth units of AU/day*)
T = {365.2, 27}(*earth moon period in days*);
EQUATIONS OF MOTION
Nbody = Table[
NDSolve[{(x''[t] + (
G*m[[i]]*x[t])/((x[t])^2 + (y[t])^2)^(3/2) + (
G*m[[j]]* (x[t]))/((x[t])^2 + (y[t])^2)^(3/2) ==
0, (y^''[
t] == -((G*m[[i]] (y[t]))/((x[t])^2 + (y[t])^2)^(3/2)) - (
G*m[[j]] (y[t]))/((x[t])^2 + (y[t])^2)^(3/2), x[0] == rx[[i]],
y[0] == ry[[i]], [x]'[0] == 0,
[y]'[0] == v[[i]]}, {x, y}, {t, 0, 365.2}], {i,
2}, {j, 2}]
I can plot the results using `
plot2D = Show[
Table[ParametricPlot[
Evaluate[{x[t], y[t]} /. Nbody[[i]]], {t, 0, T[[i]]}(*,
PlotStyle\[Rule]None*), PlotRange -> 1.1], {i, 2}]];
Animate[Show[plot2D,
Graphics[Table[{Red, PointSize[0.02],
Point[{x[t], y[t]} /. Nbody[[i]]]}, {i, 2}]]], {t, 0, 365.2},
AnimationRate -> 1, AnimationRunning -> False]
However I get this graph Along with the error message
ReplaceAll::reps: {{{{x->InterpolatingFunction[{<<1>>},{<<13>>},{<<1>>},{<<3>>},{<<1>>}],y->InterpolatingFunction[{<<1>>},{<<13>>},{<<1>>},{<<3>>},{<<1>>}]}}}[[2]]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.
Coordinate ReplaceAll[{$CellContext`x[0], $CellContext`y[0]}, Part[{{{$CellContext`x -> InterpolatingFunction[{{0., 365.2}}, {5, 7, 2, {84}, {4}, 0, 0, 0, 0, Automatic, {}, {}, False}, {{0., 0.011199207293206922`, 0.022398414586413843`, 0.27915745480600035`, 0.5359 should be a pair of numbers, or a Scaled or Offset form.
What I think is wrong
From the graph, the orbit of the earth is correct. The orbit of the moon should be around the earth which obviously it's not!
From section of Code for Nbody;
The masses for the sun and the moon are correct However I think the respective vectors for each body are not.
What I think I need to do is find a way for apply the initial conditions to both bodies individually, which I don't think I've done using this code.
Would anyone be able to help in correcting this?
Thank you for your time.
planet
is used in theGraphics[ ... ]
expression, butplanet
is not defined in the code above. $\endgroup$(x^\[Prime]\[Prime])[t]
is wrong, output only looks likeDerivative[2][x][t]
! That's whyNBody
isn't evaluated! $\endgroup$NBodySimulation
instead? $\endgroup$