10
$\begingroup$

I have seen a lot of solar system models (both animated and not), but I am trying to replicate the same thing with moons. Either just one planet and its lunar system or the our whole solar system with the respective moons. However the "OrbitPath" function does not seem to work with these.

I am able to get the OrbitCircumference and other data like that, but I am not sure how to implement this into a system. Any help would be highly appreciated.

Currently I have the following code to determine the circumference of the moons around their respective planets:

data = PlanetaryMoonData[EntityClass["PlanetaryMoon", "EarthMoon"], 
  "OrbitCircumference", "EntityAssociation"]; data2 = 
 PlanetaryMoonData[EntityClass["PlanetaryMoon", "MarsMoon"], 
  "OrbitCircumference", "EntityAssociation"];
data3 = PlanetaryMoonData[EntityClass["PlanetaryMoon", "JupiterMoon"],
    "OrbitCircumference", "EntityAssociation"];
data4 = PlanetaryMoonData[EntityClass["PlanetaryMoon", "SaturnMoon"], 
   "OrbitCircumference", "EntityAssociation"];
data5 = PlanetaryMoonData[EntityClass["PlanetaryMoon", "UranusMoon"], 
   "OrbitCircumference", "EntityAssociation"];
data6 = PlanetaryMoonData[EntityClass["PlanetaryMoon", "NeptuneMoon"],
    "OrbitCircumference", "EntityAssociation"];
UnitConvert[{data, data2, data3, data4, data5, 
  data6}, "astronomical units"]

but I don't know how I would get something like this, but with moons rather than planets:

Graphics3D[PlanetData[PlanetData[], "OrbitPath"]]

I also tried this

Graphics3D[MinorPlanetData[MinorPlanetData["Mars"], "OrbitPath"]]

but all i got was errors and it took ages to run.

$\endgroup$
8
  • 1
    $\begingroup$ Could you be more specific in your question? Please show your due diligence, share the code you are working on. Only good questions are likely to get great answers. Please edit your question to improve it making it more specificly about Mathematica programming. Otherwise it may be closed. $\endgroup$
    – rhermans
    Feb 4, 2016 at 16:49
  • 5
    $\begingroup$ AstronomicalData["Moon", "Properties"] for getting the properties of an object $\endgroup$ Feb 4, 2016 at 16:50
  • 3
    $\begingroup$ You'll have to compute the ellipse yourself based on AstronomicalData["Ganymede", "OrbitRules"] or similar. Or you can use NASA's SPICE libraries and data if this is more than a casual query. $\endgroup$
    – user1722
    Feb 4, 2016 at 20:57
  • 1
    $\begingroup$ I'm guessing you would use some combination of ParametricPlot3D and Animate? It would really depend on how you want the model to look. $\endgroup$
    – user1722
    Feb 6, 2016 at 15:53
  • 1
    $\begingroup$ There's a fair bit of interest in Jupiter's 4 largest moons, so I'm somewhat surprised no one has done this (at least for Jupiter + the 4 big moons) yet. BTW, feel free to contact me (see my profile) if you want help in real time, and you can post the results here. $\endgroup$
    – user1722
    Feb 6, 2016 at 17:25

1 Answer 1

5
$\begingroup$

Please provide credit if you use (parts of) this code anywhere, I spent a lot of time on it.

The Ephemeris data comes from The NASA HORIZONS Web-Interface, and has been generated for 00:00-1-1-2000, with respect to the Barycenters of the Solar system, Jupiter System, and Saturn System, using the VECTORS Ephemeris type.

My answer here doesn't rely PlanetaryMoonData[], the reason is that it is built around an all-purpose physical interaction system which I have also used for fictional physical systems.

Where possible I use AU units with Standard gravitational parameters instead of Mass, point is to be consistent. AU is more accurate than SI. For clarity's sake I'll put the actual init data in pastebin, please paste the result here into an input cell:

LINK TO FIRST MMA CELL

As J.M. pointed out, you can import this directly into MMA with

Import["http://pastebin.com/raw/8PpfwbXe", "Package"]

This is the actual engine: it is a function which inputs the data in the form of the pastebin file: As you may want to generate a system for many different particles/ planets/ moons, there is a progress monitor:

generate[dataset_, rep_, nn_, {xxm_, yym_, zzm_}, mm_, tmax_] := 
 Module[{n, positions, velocities, m, xt, xz, yt, yz, zt, zz, xm, ym, 
   zm, rm, xf1, xf, yf1, yf, zf1, zf, pos, vel},
  (*Setting up basic data*)
  n = Length[dataset];
  nn = n;
  positions = Transpose[Table[dataset[[i, 2]], {i, n}]];
  velocities = Transpose[Table[dataset[[i, 3]], {i, n}]];
  m = Table[dataset[[i, 4]], {i, n}];
  mm = m;
  (*Setting up stuff for NDSolve[]*)
  {xt, yt, 
    zt} = {ToExpression[Table["x" <> ToString[i], {i, n}]], 
    ToExpression[Table["y" <> ToString[i], {i, n}]], 
    ToExpression[Table["z" <> ToString[i], {i, n}]]};
  {xm, ym, zm} = {Through[xt[t]], Through[yt[t]], Through[zt[t]]};
  {xz, yz, zz} = {Through[xt[0]], Through[yt[0]], Through[zt[0]]};
  yt = ToExpression[Table["y" <> ToString[i], {i, n}]];
  zt = ToExpression[Table["z" <> ToString[i], {i, n}]];
  rm = Flatten[
     Table[If[i != j, 
       Sqrt[(xm[[j]] - xm[[i]])^2 + (ym[[j]] - ym[[i]])^2 + (zm[[j]] -
             zm[[i]])^2]], {i, n}, {j, n}]] /. Null -> Sequence[];
  (*The final equations for the differential equation*)

  xf1 = (Flatten[
       Table[If[i != j, (xm[[j]] - xm[[i]]) m[[j]]], {i, n}, {j, 
         n}]] /. Null -> Sequence[])/rm^3;
  xf = Thread[
    D[D[xm, t], t] == 
     Table[Total[Take[xf1, {(n - 1) i - n + 2, (n - 1) i}]], {i, 
       n}]];
  yf1 = (Flatten[
       Table[If[i != j, (ym[[j]] - ym[[i]]) m[[j]]], {i, n}, {j, 
         n}]] /. Null -> Sequence[])/rm^3;
  yf = Thread[
    D[D[ym, t], t] == 
     Table[Total[Take[yf1, {(n - 1) i - n + 2, (n - 1) i}]], {i, 
       n}]];
  zf1 = (Flatten[
       Table[If[i != j, (zm[[j]] - zm[[i]]) m[[j]]], {i, n}, {j, 
         n}]] /. Null -> Sequence[])/rm^3;
  {xxm, yym, zzm} = {xm, ym, zm};
  zf = Thread[
    D[D[zm, t], t] == 
     Table[Total[Take[zf1, {(n - 1) i - n + 2, (n - 1) i}]], {i, n}] ];
  pos = {Thread[xz == positions[[1]]], Thread[yz == positions[[2]]], 
    Thread[zz == positions[[3]]]};
  vel = {Thread[D[xm, t] == velocities[[1]]] /. t -> 0, 
    Thread[D[ym, t] == velocities[[2]]] /. t -> 0, 
    Thread[D[zm, t] == velocities[[3]]] /. t -> 0};
  Print["Now generating Interpolation functions for ", 
   dataset[[All, 1]], "..."];
  rep = Monitor[NDSolve[Flatten[{xf, yf, zf,
       vel, pos}], Flatten[{xt, yt, zt}], {t, 0, tmax}, 
     MaxSteps -> \[Infinity], MaxStepSize -> 0.5, 
     EvaluationMonitor :> (time = t)], ProgressIndicator[time/tmax]];
  Print["Done."]]

Now we input the data. dataSol is the data generate of the pastebin cell, the solution gets stored in s, n stores the amount of particles, {xm, ym, zm} gets the InterpolationFunctions[], m gets the relative masses, for plotting purposes, and you need to set tmax the time in days. In the first case a sets the amount of years the solar ystem is computed.

a = 256;
Clear[s, n, xm, ym, zm, m]
generate[dataSol, s, n, {xm, ym, zm}, m, 
 tmax1 = 365*a + Ceiling[a/4]]
Clear[s2, n2, xm2, ym2, zm2, m2]
generate[dataJupiter, s2, n2, {xm2, ym2, zm2}, m2, tmax2 = 20]
Clear[s3, n3, xm3, ym3, zm3, m3]
generate[dataSaturn, s3, n3, {xm3, ym3, zm3}, m3, tmax3 = 20]

This is an example of how to plot it: Top-left: Solar System upto Ceres, top-right: Solar System starting at Jupiter. Bottom-left: Jupiter system, bottom-right: Saturn system.

Grid[{{Manipulate[
    Show[ParametricPlot3D[
      Evaluate[
       Flatten[Thread[Transpose[{xm, ym, zm}][[1 ;; 7]] /. s], 
        1]], {t, t0, a}, 
      PlotRange -> {{-2.5, 2.5}, {-2.5, 2.5}, {-0.75, 0.75}}, 
      PlotStyle -> ps[[1 ;; 7]], Background -> Black, 
      ImageSize -> Large], 
     Graphics3D[
      Transpose[{colours[[1 ;; 7]], 
        Thread[Sphere[
          Evaluate[
           Flatten[
            Thread[Transpose[{xm, ym, zm} /. t -> a][[1 ;; 7]] /. s], 
            1]], sizes1a[[1 ;; 7]]]]}]]], {t0, 0, Min[tmax1, 700] - 1,
      1}, {{a, Min[tmax1, 700], "t"}, t0 + 1, Min[tmax1, 700], 1}], 
   Manipulate[
    Show[ParametricPlot3D[
      Evaluate[
       Flatten[Thread[
         Transpose[{xm, ym, zm}][[{1, 8, 9, 10, 11, 12}]] /. s], 
        1]], {t, t0, a}, 
      PlotRange -> {{-45, 45}, {-45, 45}, {-10, 10}}, 
      PlotStyle -> ps[[{1, 8, 9, 10, 11, 12}]], Background -> Black, 
      ImageSize -> Large], 
     Graphics3D[
      Transpose[{colours[[{1, 8, 9, 10, 11, 12}]], 
        Thread[Sphere[
          Evaluate[
           Flatten[
            Thread[Transpose[{xm, ym, zm} /. t -> a][[{1, 8, 9, 10, 
                11, 12}]] /. s], 1]], 
          sizes1b[[{1, 8, 9, 10, 11, 12}]]]]}]]], {t0, 0, tmax1 - 1, 
     1}, {{a, tmax1, "t"}, t0 + 1, tmax1, 1}]}, {Manipulate[
    Show[ParametricPlot3D[
      Evaluate[
       Flatten[Thread[Transpose[{xm2, ym2, zm2}] /. s2], 1]], {t, t0, 
       a}, PlotRange -> {{-0.02, 0.02}, {-0.02, 0.02}, {-0.005, 
         0.005}}, PlotStyle -> ps2, Background -> Black, 
      ImageSize -> Large], 
     Graphics3D[
      Transpose[{colours2, 
        Thread[Sphere[
          Evaluate[
           Flatten[Thread[Transpose[{xm2, ym2, zm2} /. t -> a] /. s2],
             1]], sizes2]]}]]], {t0, 0, tmax2 - 1, 
     1}, {{a, tmax2, "t"}, t0 + 1, tmax2, 1}], 
   Manipulate[
    Show[ParametricPlot3D[
      Evaluate[
       Flatten[Thread[Transpose[{xm3, ym3, zm3}] /. s3], 1]], {t, t0, 
       a}, PlotRange -> {{-0.02, 0.02}, {-0.02, 0.02}, {-0.005, 
         0.005}}, PlotStyle -> ps3, Background -> Black, 
      ImageSize -> Large], 
     Graphics3D[
      Transpose[{colours3, 
        Thread[Sphere[
          Evaluate[
           Flatten[Thread[Transpose[{xm3, ym3, zm3} /. t -> a] /. s3],
             1]], sizes3]]}]]], {t0, 0, tmax3 - 1, 
     1}, {{a, tmax3, "t"}, t0 + 1, tmax3, 1}]}}]

enter image description here

$\endgroup$
4
  • $\begingroup$ Where do the numbers in dataSol, dataJupiter, etc. come from? $\endgroup$
    – shrx
    Feb 3, 2017 at 11:45
  • $\begingroup$ @shrx Good point, I've added that to the top of the answer. $\endgroup$
    – Feyre
    Feb 3, 2017 at 11:53
  • 1
    $\begingroup$ If you want to load the data without needing to copy-pasta: Import["http://pastebin.com/raw/8PpfwbXe", "Package"]. $\endgroup$ Feb 3, 2017 at 14:41
  • $\begingroup$ @J.M. Thank you, I've added that to the answer. $\endgroup$
    – Feyre
    Feb 3, 2017 at 15:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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