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"]
Show
. About your code, what{EarthSphere, {0, 0, 0}}
is? there is missing,
between plots too. $\endgroup$