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

Question

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"]

Answer 1

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