# 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[t], y[t], z[t]}, {x[t], y[t],
z[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 = AstronomicalData["Earth", "Mass"];
tmax = 20000;
rx = (r + 300000 ) Cos[45 Degree];
ry = (r + 300000 ) Sin[45 Degree];
rz = 0;
rx = (r + 600000 ) Cos[90 Degree];
ry = (r + 600000 ) Sin[90 Degree];
rz = 0;
vx = Sqrt[(G  m)/(r + 300000)] Sin[45 \[Degree]]
vy = Sqrt[(G  m)/(r + 300000)] Cos[45 \[Degree]]
vz = 0
vx = Sqrt[(G  m)/(r + 600000)] Sin[90 \[Degree]]
vy = Sqrt[(G  m)/(r + 600000)] Cos[90 \[Degree]]
vz = 0
soln = NDSolve[{
x''[t] == -((
G m x[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2)),
y''[t] == -((
G m y[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2)),
z''[t] == -((
G m z[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2)),
x''[t] == -((
G m x[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2)),
y''[t] == -((
G m y[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2)),
z''[t] == -((
G m z[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2)),

x == rx, y == ry, z == rz,
x == rx, y == ry, z == rz,
x' == -vx, y' == vy, z' == 0,
x' == -vx, y' == vy, z' == 0}, {x[t],
y[t], z[t], x[t], y[t], z[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 '14 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. – InquisitiveInquirer Jan 26 '14 at 19:07

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[t], y[t], z[t]}, {x[t], y[t], z[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! – InquisitiveInquirer Jan 26 '14 at 19:46