Travelling around Jupiter

Since we have soon historical conjunction of Jupiter and Saturn I have prepared code with visualization planet Jupiter and Saturn with moons. This code first been prepared for version 11.3 and with this version it can be used as it is. Visualization of planet

JSP = QuantityMagnitude[
UnitConvert[
PlanetaryMoonData[{"Io", "Europa", "Ganymede",
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;
"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];

2 Pi}, {p, 0, π}, Mesh -> None, PlotStyle -> Texture[texture],
Lighting -> "Neutral", Boxed -> False, Axes -> False,
PlotPoints -> 100]


We have out this nice picture

Visualization of travelling orbit (we solve restricted 3 body problem)

    M1 = QuantityMagnitude[Entity["Planet", "Jupiter"]["Mass"]];
M2 = QuantityMagnitude[Entity["Star", "Sun"]["Mass"]];
Ra = QuantityMagnitude[
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]}];
{Orbit = ParametricPlot3D[V, {t, 0, tm}], Plot[Norm[V], {t, 0, tm}]}

Show[{Graphics3D[Orbit[[1]], Boxed -> False], planet}]


We have this picture with Jupiter and orbit around

Now we combine orbit, planet and moons in one scene. And here we have first problem with "OrbitPeriod" for moons. In version 12.0.0.0 it gives out period in h for "Io", "Europa", and in d for "Ganymede", "Callisto". Therefore we need to redefine b as

b1 = {8.81749715291124946277507323.9999565727231374*^6/3600,
4.3926084098845395732430623.9999565727231374*^6/3600,
5.232084683773584907638743483.9999565727231374*^7/24/3600,
2.242993763475134807294235213.9999565727231374*^7/24/3600}; Th =
tm*31557600*
QuantityMagnitude[
Entity["Planet", "Jupiter"]["OrbitPeriod"]]/(3600)/(2*Pi);


Now we cam prepare scene and animation with traveler, planet and moons rotation

frames=Table[Show[{Graphics3D[
Rotate[planet[[1]], 2*Pi*Th*t/9.841666666666667/tm, {0, 0, 1}],
Boxed -> False, Axes -> False, ImageSize -> .25 {1920, 1080}],
Table[Graphics3D[{White,
Sphere[{R[[i]]*Cos[a[[i]] + b1[[i]]*t],
R[[i]]*Sin[a[[i]] + b1[[i]]*t], 0}, k*c[[i]]]},
Background -> Black, Boxed -> False, PlotRange -> All], {i, 1,
4}]}, Background -> Black, 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 - .0025)*tm, .0025*tm}];


Finally we make animation with ListAnimate[frames] or export it as a .gif

The problem is that this code not working with versions 12.1 and 12.2. How we can adapted code for these versions?

Update1. Second code with visualization of Saturn and moons "Mimas", "Enceladus", "Tethys", "Dione", "Rhea", "Titan" is differ from above since we need to show rings and also orbit around Saturn is not like orbit around Jupiter. The code is in notebook attached here, and animation looks like this one

• Version 12.1 It seems to work. Start a new kernel and try again. Commented Dec 20, 2020 at 11:46
• Awesome example, my thanks and +1! Commented Dec 20, 2020 at 13:24
• Thanks for doing this! I tried the code out with Mathematica Version 12.2 Mac OS 10.15.7 -- the code and the animation making work fine. Please consider making a dedicated post at community.wolfram.com. Commented Dec 20, 2020 at 13:27
• @AntonAntonov I will do at community.wolfram.com also for Saturn, Uranus and Neptune. But there is still question about b. In v.11.3 "OrbitPeriod" been in second in SI, but in 12.0, 12.1, 12.2 it is in hours for some moons and in days for other. Why we need to define b1. It looks like UnitConvert[] not working any more. Commented Dec 20, 2020 at 14:19
• From the documentation for UnitConvert, "Conversion to a unit system will return a unit that exists in that unit system, not necessarily the base unit of that system". However, "UnitConvert[q, "SIBase"] typically converts all units to their SI base units". Consequently, UnitConvert[PlanetaryMoonData[{"Io","Europa","Ganymede","Callisto"},{"OrbitPeriod","SemimajorAxis","Radius"}], "SIBase"] would provide the periods in seconds; however, the distances would then be in meters rather than kilometers. Commented Dec 20, 2020 at 21:20

Thanks to this discussion I been able to develop some general approach to the problem. In this answer I have demonstrated code to simulate journey around Uranus and Neptune. Code for Uranus

JSP = QuantityMagnitude[
UnitConvert[
PlanetaryMoonData[{"Miranda", "Ariel", "Umbriel", "Titania",
"SI"] /.
q : Quantity[_, "Days" | "Hours"] :>
UnitConvert[q,
"Seconds"]];(*Miranda,Ariel,Umbriel,Titania,and Oberon.*)
js = Length[JSP];
a = Pi*RandomReal[{0, 2}, js];
b = 31557600*
QuantityMagnitude[Entity["Planet", "Uranus"]["OrbitPeriod"]]*
Table[1/JSP[[i, 1]], {i, js}];

R = Table[JSP[[i, 2]], {i, js}];
c = Table[JSP[[i, 3]], {i, js}];
k = 3;
"Kilometers"];
Xoblateness = Entity["Planet", "Uranus"]["Oblateness"];
Xobliquity = Entity["Planet", "Uranus"]["Obliquity"];
distanse =
QuantityMagnitude[
Entity["Planet", "Uranus"]["AverageOrbitDistance"], "Kilometers"];
angularvelocity =
Entity["Planet", "Uranus"]["EquatorialAngularVelocity"];
texture =
ImageReflect[
EntityValue[Entity["Planet", "Uranus"],
"CylindricalEquidistantTexture"], Bottom];
2 Pi}, {p, 0, \[Pi]}, Mesh -> None, PlotStyle -> Texture[texture],
Lighting -> "Neutral", Boxed -> False, Axes -> False,
PlotPoints -> 100]


Next we need to define rim as image for rings (just copy and past image on the line I0 = SetAlphaChannel[rim, .5] instead of rim)

ringinnerradius = 34840;

Graphics3D[{White,
Boxed -> False];
I0 = SetAlphaChannel[Import["https://i.sstatic.net/KoALx.jpg"], .5];
T = Table[{ImageValue[I0, {i, j}][[1]], ImageValue[I0, {i, j}][[2]],
ImageValue[I0, {i, j}][[
3]], (ImageValue[I0, {i, j}][[1]] + ImageValue[I0, {i, j}][[2]] +
ImageValue[I0, {i, j}][[3]])/1.5}, {j, 5, 30}, {i, 1, 748}];
rings = ParametricPlot3D[{r Cos[t], r Sin[t], 0}, {r, ringinnerradius,
ringouterradius}, {t, 0, 2 Pi}, Mesh -> None,
PlotStyle -> Texture[T], PlotPoints -> 100, MaxRecursion -> 2];
Th = 31557600*
QuantityMagnitude[
Entity["Planet", "Uranus"]["OrbitPeriod"]]/(3600); Trp =
QuantityMagnitude[Entity["Planet", "Uranus"]["RotationPeriod"]];
gr[t_] :=
Show[{Graphics3D[Rotate[planet[[1]], 2*Pi*Th*t/Trp, {0, 0, 1}],
Boxed -> False, Axes -> False, ImageSize -> .3 {1920, 1080}],
Graphics3D[{Lighting -> {{"Ambient", GrayLevel[.53]}},
rings[[1]]}], Sun}, Background -> Black,
ImageSize -> .4 {1920, 1080}, SphericalRegion -> True,
ViewAngle -> Pi/10,
PlotRange -> All,
Lighting -> {{"Ambient", GrayLevel[0]}, {"Point",

M1 = QuantityMagnitude[Entity["Planet", "Uranus"]["Mass"]];
M2 = QuantityMagnitude[Entity["Star", "Sun"]["Mass"]];
Ra = QuantityMagnitude[
S = QuantityMagnitude[Entity["Planet", "Uranus"]["SemimajorAxis"]];

ar = 3 Ra/(S*149597870700);
m = M1/M2;
G = 6.67384*10^(-11); vS =
Sqrt[G*M1/(3 Ra)]; u0 = 6800; {ux, uy,
uz} = -vS/u0 {-.1, 0, 1.2}; tm = 0.00065;
{x0, y0, z0} = {1 - m - ar, 0., 0.};
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]}];

{Orbit = ParametricPlot3D[V, {t, 0, tm}], Plot[Norm[V], {t, 0, tm}]};
Show[{Graphics3D[Orbit[[1]], Boxed -> False], planet}]


And finally we have animation

Export["C:\\...\\Uranus.gif",
Table[Show[{gr[t],
Table[Graphics3D[{White,
Sphere[{R[[i]]*Cos[a[[i]] + b[[i]]*t],
R[[i]]*Sin[a[[i]] + b[[i]]*t], 0}, 2.5*c[[i]]]},
Background -> Black, Boxed -> False,
ImageSize -> .3 {1920, 1080}, PlotRange -> All], {i, 1, 5}]},
Background -> Black, ImageSize -> .3 {1920, 1080},
SphericalRegion -> True, ViewAngle -> Pi/4, PlotRange -> All,
Lighting -> {{"Ambient", GrayLevel[0.05]}, {"Point",
White, {0, 0, 20 radius}}}, ViewVector -> {V, {0, 0, 0}}], {t,
0, (1 - .00501)*tm, .005*tm}], AnimationRepetitions -> Infinity]


Code for Neptune

JSP = QuantityMagnitude[UnitConvert[PlanetaryMoonData[{Entity["PlanetaryMoon", "Naiad"], Entity["PlanetaryMoon", "Thalassa"], Entity["PlanetaryMoon", "Despina"],
Entity["PlanetaryMoon", "Galatea"], Entity["PlanetaryMoon", "Larissa"], Entity["PlanetaryMoon", "S2004N1"], Entity["PlanetaryMoon", "Proteus"],
Entity["PlanetaryMoon", "Triton"], Entity["PlanetaryMoon", "Nereid"], Entity["PlanetaryMoon", "Halimede"], Entity["PlanetaryMoon", "Sao"],
Entity["PlanetaryMoon", "Laomedeia"], Entity["PlanetaryMoon", "Psamathe"], Entity["PlanetaryMoon", "Neso"]}, {"OrbitPeriod", "SemimajorAxis", "Radius"}],
"SI"] /. q:Quantity[_, "Days" | "Hours"] :> UnitConvert[q, "Seconds"]]; js = 8;
a = Pi*RandomReal[{0, 2}, js];
b = 31557600*QuantityMagnitude[Entity["Planet", "Neptune"]["OrbitPeriod"]]*Table[1/JSP[[i,1]], {i, 1, js}];
R = Table[JSP[[i,2]], {i, 1, js}];
c = Table[JSP[[i,3]], {i, 1, js}];
k = 3;
Xoblateness = Entity["Planet", "Neptune"]["Oblateness"];
Xobliquity = Entity["Planet", "Neptune"]["Obliquity"];
distanse = QuantityMagnitude[Entity["Planet", "Neptune"]["AverageOrbitDistance"], "Kilometers"];
angularvelocity = Entity["Planet", "Neptune"]["EquatorialAngularVelocity"];
texture = ImageReflect[EntityValue[Entity["Planet", "Neptune"], "CylindricalEquidistantTexture"], Bottom];
PlotStyle -> Texture[texture], Lighting -> "Neutral", Boxed -> False, Axes -> False, PlotPoints -> 100]


Image rim for rings visualization

ringinnerradius = 40900;

Graphics3D[{White,
Boxed -> False];
I0 = SetAlphaChannel[Import["https://i.sstatic.net/YUp2p.jpg"], .5];
T = Table[{ImageValue[I0, {i, j}][[1]], ImageValue[I0, {i, j}][[2]],
ImageValue[I0, {i, j}][[
3]], (ImageValue[I0, {i, j}][[1]] + ImageValue[I0, {i, j}][[2]] +
ImageValue[I0, {i, j}][[3]])/1.5}, {j, 2, 28}, {i, 1, 753}];
rings = ParametricPlot3D[{r Cos[t], r Sin[t], 0}, {r, ringinnerradius,
ringouterradius}, {t, 0, 2 Pi}, Mesh -> None,
PlotStyle -> Texture[T], PlotPoints -> 100, MaxRecursion -> 2];

grn[t_] :=
Show[{Graphics3D[Rotate[planet[[1]], 2*Pi*Th*t/Trp, {0, 0, 1}],
Boxed -> False, Axes -> False, ImageSize -> .3 {1920, 1080}],
Graphics3D[{Lighting -> {{"Ambient", GrayLevel[.53]}},
rings[[1]]}], Sun, Shadow}, Background -> Black,
ImageSize -> .3 {1920, 1080}, SphericalRegion -> True,
ViewAngle -> Pi/10,
PlotRange -> All,
Lighting -> {{"Ambient", GrayLevel[0]}, {"Point",

M1 = QuantityMagnitude[Entity["Planet", "Neptune"]["Mass"]];
M2 = QuantityMagnitude[Entity["Star", "Sun"]["Mass"]];
Ra = QuantityMagnitude[
S = QuantityMagnitude[Entity["Planet", "Neptune"]["SemimajorAxis"]];

ar = 8 Ra/(S*149597870700);
m = M1/M2;
G = 6.67384*10^(-11); vS =
Sqrt[G*M1/(8 Ra)]; u0 = 6800; {ux, uy, uz} =
vS/u0 {0, 1., .01}; tm = 0.000323;
{x0, y0, z0} = {1 - m - ar, 0., 0.};
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[Pi/24, {1, 0, 0}]; V =
rt[S*149597870.700*{xfun[t] - 1 + m, yfun[t], zfun[t]}];


Animation

Th = 31557600*
QuantityMagnitude[
Entity["Planet", "Neptune"]["OrbitPeriod"]]/(3600); Trp =
QuantityMagnitude[Entity["Planet", "Neptune"]["RotationPeriod"]];

Export["C:\\Users\\...\\NeptuneM8.gif", Table[
Show[{grn[t],
Table[Graphics3D[{White,
Sphere[{R[[i]]*Cos[a[[i]] + b[[i]]*t],
R[[i]]*Sin[a[[i]] + b[[i]]*t], 0}, 3*c[[i]]]},
Background -> Black, Boxed -> False,
ImageSize -> .3 {1920, 1080}, PlotRange -> All], {i, 1, js}]},
Background -> Black, ImageSize -> .3 {1920, 1080},
SphericalRegion -> True, ViewAngle -> Pi/4, PlotRange -> All,
ViewVector -> {V, {0, 0, 0}}], {t, 0, (1 - .005)*tm, .005*tm}],
AnimationRepetitions -> Infinity]