Опубликую код Используйте Mathematica 11.3 для создания последнего рисунка.
JSP = QuantityMagnitude[
UnitConvert[
PlanetaryMoonData[{"Io", "Europa", "Ganymede",
"Callisto"}, {"OrbitPeriod", "SemimajorAxis", "Radius"}], "SI"]];
a = Pi*{RandomReal[{0, 2}], RandomReal[{0, 2}], RandomReal[{0, 2}],
RandomReal[{0, 2}]};
b = 31557600*
QuantityMagnitude[Entity["Planet", "Jupiter"]["OrbitPeriod"]]*
Table[1/JSP[[i, 1]], {i, 1, 4}];
R = Table[JSP[[i, 2]], {i, 1, 4}];
c = Table[JSP[[i, 3]], {i, 1, 4}];
k = 3;
radius = QuantityMagnitude[
Entity["Planet", "Jupiter"]["EquatorialRadius"], "Kilometers"];
radius1 =
QuantityMagnitude[Entity["Planet", "Jupiter"]["PolarRadius"],
"Kilometers"];
Xoblateness = Entity["Planet", "Jupiter"]["Oblateness"];
Xobliquity = Entity["Planet", "Jupiter"]["Obliquity"];
distanse =
QuantityMagnitude[
Entity["Planet", "Jupiter"]["AverageOrbitDistance"], "Kilometers"];
angularvelocity =
Entity["Planet", "Jupiter"]["EquatorialAngularVelocity"];
texture =
ImageReflect[
EntityValue[Entity["Planet", "Jupiter"],
"CylindricalEquidistantTexture"], Bottom];
planet = ParametricPlot3D[{radius Cos[t] Sin[p],
radius Sin[t] Sin[p], (1 - Xoblateness) radius Cos[p]}, {t, 0,
2 Pi}, {p, 0, \[Pi]}, Mesh -> None, PlotStyle -> Texture[texture],
Lighting -> "Neutral", Boxed -> False, Axes -> False,
PlotPoints -> 100];
M1 = QuantityMagnitude[Entity["Planet", "Jupiter"]["Mass"]];
M2 = QuantityMagnitude[Entity["Star", "Sun"]["Mass"]];
Ra = QuantityMagnitude[
Entity["Planet", "Jupiter"]["EquatorialRadius"]]*10^3;
S = QuantityMagnitude[Entity["Planet", "Jupiter"]["SemimajorAxis"]];
ar = Ra/(S*149597870700);
m = M1/M2;
G = 6.67384*10^(-11); vS =
Sqrt[G*M1/Ra]; u0 = 16000; {ux, uy,
uz} = {0.28309789428056364`, -2.269693660804425`,
0.24765897137821516`}; tm = 2*0.007188414157797473;
{x0, y0, z0} = {0.9987450742768228`, -0.00006988474848081634`, \
-0.000015444353758731164`};
eq = {x''[t] ==
2*y'[t] + x[t] - (
m (-1 + m + x[t]))/((-1 + m + x[t])^2 + y[t]^2 + z[t]^2)^(
3/2) - ((1 - m) (m + x[t]))/((m + x[t])^2 + y[t]^2 + z[t]^2)^(
3/2), y''[t] == -2*x'[t] + y[t] - (
m y[t])/((-1 + m + x[t])^2 + y[t]^2 + z[t]^2)^(
3/2) - ((1 - m) y[t])/((m + x[t])^2 + y[t]^2 + z[t]^2)^(3/2),
z''[t] == -((m z[t])/((-1 + m + x[t])^2 + y[t]^2 + z[t]^2)^(
3/2)) - ((1 - m) z[t])/((m + x[t])^2 + y[t]^2 + z[t]^2)^(3/2),
x[0.] == x0, y[0.] == y0, x'[0.] == ux, y'[0.] == uy, z'[0.] == uz,
z[0.] == z0};
{xfun, yfun, zfun} = NDSolveValue[eq, {x, y, z}, {t, 0, tm}];
rt = RotationTransform[Xobliquity, {1, 0, 0}]; V =
rt[S*149597870.700*{xfun[t] - 1 + m, yfun[t], zfun[t]}];
Export["C:\\Users\\...\\Desktop\\Jupiter.gif",
Table[Show[{Graphics3D[
Rotate[planet[[1]],
2*Pi*62.91486673636752*t/9.841666666666667/tm, {0, 0, 1}],
Boxed -> False, Axes -> False],
Table[Graphics3D[{White,
Sphere[{R[[i]]*Cos[a[[i]] + b[[i]]*t],
R[[i]]*Sin[a[[i]] + b[[i]]*t], 0}, k*c[[i]]]},
Background -> Black, Boxed -> False,
ImageSize -> .4 {1920, 1080}, PlotRange -> All], {i, 1, 4}]},
Background -> Black, ImageSize -> .4 {1920, 1080},
SphericalRegion -> True, ViewAngle -> Pi/4, PlotRange -> All,
Lighting -> {{"Ambient", GrayLevel[0.05]}, {"Point",
White, {20 radius, 0, 0}}}, ViewVector -> {V, {0, 0, 0}}], {t,
0, (1 - .00125)*tm, .00125*tm}], AnimationRepetitions -> Infinity]