# Combining Newtonian 3D N-Body Simulation Output and Earth Sphere

I'm trying to figure out how to make a simulation whereby two satellites orbit Earth, but have run into a little bit of a snag. I have a solution to the system of ODEs for the n-body output, and I found a brilliant thread explaining how to create an Earth sphere here (How to make a 3D globe?) but I'm struggling to combine the output of the two. This is my attempt: (Note that "soln" is the output of the n-body NDSolve computations, which I didn't put in here as I thought combining the two entities would be a trivial syntactic edit. If the full code is needed I'll be happy to supply it as well)

Earth = Import["Desktop/earthtruecolor_nasa_big.jpg"];
EarthSphere =
ParametricPlot3D[{Cos[u] Sin[v], Sin[u] Sin[v], Cos[v]}, {u, 0,
2 Pi}, {v, 0, Pi}, Mesh -> None, PlotPoints -> 100,
TextureCoordinateFunction -> ({#4, 1 - #5} &), Boxed -> False,
PlotStyle -> Texture[Show[Earth, ImageSize -> 1000]],
Lighting -> "Neutral", Axes -> False, RotationAction -> "Clip",
ViewPoint -> {-2.026774, 2.07922, 1.73753418}, ImageSize -> 800];

Show[ParametricPlot3D[
Evaluate[{{x[1][t], y[1][t], z[1][t]}, {x[2][t], y[2][t],
z[2][t]}} /. soln], {t, 0, 20000}, AxesLabel -> {x, y, z},
AspectRatio -> 1, BoxRatios -> 1, PlotStyle -> Automatic,
ImageSize -> Large,
PlotRange -> {{-10000000, 10000000}, {-10000000,
10000000}, {-10000000, 10000000}}] {EarthSphere, {0, 0, 0}}]


This is my output so far:

EDIT: I added in the n-body code as I thought it would cause confusion if I did not:

G = 6.672*10^-11
m[0] = AstronomicalData["Earth", "Mass"];
tmax = 20000;
r[0] = AstronomicalData["Earth", "Radius"];
rx[1] = (r[0] + 300000 ) Cos[45 Degree];
ry[1] = (r[0] + 300000 ) Sin[45 Degree];
rz[1] = 0;
rx[2] = (r[0] + 600000 ) Cos[90 Degree];
ry[2] = (r[0] + 600000 ) Sin[90 Degree];
rz[2] = 0;
vx[1] = Sqrt[(G  m[0])/(r[0] + 300000)] Sin[45 \[Degree]]
vy[1] = Sqrt[(G  m[0])/(r[0] + 300000)] Cos[45 \[Degree]]
vz[1] = 0
vx[2] = Sqrt[(G  m[0])/(r[0] + 600000)] Sin[90 \[Degree]]
vy[2] = Sqrt[(G  m[0])/(r[0] + 600000)] Cos[90 \[Degree]]
vz[2] = 0
soln = NDSolve[{
x[1]''[t] == -((
G m[0] x[1][t])/(x[1][t]^2 + y[1][t]^2 + z[1][t]^2)^(3/2)),
y[1]''[t] == -((
G m[0] y[1][t])/(x[1][t]^2 + y[1][t]^2 + z[1][t]^2)^(3/2)),
z[1]''[t] == -((
G m[0] z[1][t])/(x[1][t]^2 + y[1][t]^2 + z[1][t]^2)^(3/2)),
x[2]''[t] == -((
G m[0] x[2][t])/(x[2][t]^2 + y[2][t]^2 + z[2][t]^2)^(3/2)),
y[2]''[t] == -((
G m[0] y[2][t])/(x[2][t]^2 + y[2][t]^2 + z[2][t]^2)^(3/2)),
z[2]''[t] == -((
G m[0] z[2][t])/(x[2][t]^2 + y[2][t]^2 + z[2][t]^2)^(3/2)),

x[1][0] == rx[1], y[1][0] == ry[1], z[1][0] == rz[1],
x[2][0] == rx[2], y[2][0] == ry[2], z[2][0] == rz[2],
x[1]'[0] == -vx[1], y[1]'[0] == vy[1], z[1]'[0] == 0,
x[2]'[0] == -vx[2], y[2]'[0] == vy[2], z[2]'[0] == 0}, {x[1][t],
y[1][t], z[1][t], x[2][t], y[2][t], z[2][t]}, {t, 0, tmax} ,
MaxSteps -> 1000000, Method -> "StiffnessSwitching"]

• Rescale your orbits to Earth radius units, so divide by ~6350000 and use Show. About your code, what {EarthSphere, {0, 0, 0}} is? there is missing , between plots too.
– Kuba
Jan 26, 2014 at 18:59
• I used {EarthSphere, {0,0,0}} as I thought it might centre the EarthSphere object created in the code at the origin, but it unfortunately did not. Jan 26, 2014 at 19:07

## 1 Answer

I have no access to your desktop so let me skip the texture:

Show[
ParametricPlot3D[{Cos[u] Sin[v], Sin[u] Sin[v], Cos[v]}, {u, 0, 2 Pi}, {v, 0, Pi},
PlotPoints -> 100, Boxed -> False, Lighting -> "Neutral",
Axes -> False, Mesh -> {10, 0}],
ParametricPlot3D[Evaluate[{{x[1][t], y[1][t], z[1][t]}, {x[2][t], y[2][t], z[2][t]}
} /. soln]/6350000, {t, 0, 2000},
AxesLabel -> {x, y, z}, AspectRatio -> 1, BoxRatios -> 1,
PlotStyle -> Thick, ImageSize -> Large],
PlotRange -> 1.1
]


ok, I've found one old map ;p

• Great, thank you Kuba, works like a charm! Jan 26, 2014 at 19:46