Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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: enter image description here

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"]
share|improve this question
    
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 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. –  user7388 Jan 26 at 19:07

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

enter image description here

ok, I've found one old map ;p enter image description here

share|improve this answer
    
Great, thank you Kuba, works like a charm! –  user7388 Jan 26 at 19:46

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.