- Thread starter
- #1
Can you post the ODE's and your Mathematica code so far?
Numerical Simulation of the 2-Body Problem
ClearAll["Global`*"];
We begin with of the the necessary data for this problem....
M = 5974*10^21; (* mass of Earth, kg *)
m = 1000; (* mass of spacecraft, kg *)
\[Mu] = 3.986*10^5; (* gravitaional parameter, based on km units of length, \
km/s for velocity *)
Rearth = 6378; (* radius of the Earth, km *)
Simulation Inputs
r0 = {3950.55, 43197.9, 0};(* initial position vector, km *)
v0 = {3.3809, -7.25046, 0}; (* initial velocity vector, km *)
Days = 1/10; (* elapsed time of simulation days *)
\[CapitalDelta]t = Days*24*3600;(* convert elapsed days to seconds *)
s = NDSolve[
{
x1''[t] == -(\[Mu]/(Sqrt[x1[t]^2 + x2[t]^2 + x3[t]^2])^3)*x1[t],
x2''[t] == -(\[Mu]/(Sqrt[x1[t]^2 + x2[t]^2 + x3[t]^2])^3)*x2[t],
x3''[t] == -(\[Mu]/(Sqrt[x1[t]^2 + x2[t]^2 + x3[t]^2])^3)*x3[t],
x1[0] == r0[[1]], (* intial x-position of satellite *)
x2[0] == r0[[2]],(* intial y-position of satellite *)
x3[0] == r0[[3]],(* intial y-position of satellite *)
x1'[0] == v0[[1]],(* intial vx-rel of satellite *)
x2'[0] == v0[[2]],(* intial vy-rel of satellite *)
x3'[0] == v0[[3]](* intial vy-rel of satellite *)
},
{x1, x2, x3},
{t, 0, \[CapitalDelta]t}
];
Plot of the Trajectory Relative to Earth
g1 = ParametricPlot3D[
Evaluate[{x1[t], x2[t], x3[t]} /. s], {t, 0, \[CapitalDelta]t},
PlotStyle -> {Red, Thick}];
g2 = Graphics3D[{Blue, Opacity[0.6], Sphere[{0, 0, 0}, Rearth]}];
Show[g2, g1, Boxed -> False]