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}]}}]
AstronomicalData["Moon", "Properties"]
for getting the properties of an object $\endgroup$AstronomicalData["Ganymede", "OrbitRules"]
or similar. Or you can use NASA's SPICE libraries and data if this is more than a casual query. $\endgroup$