22
$\begingroup$

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", 
      "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, π}, Mesh -> None, PlotStyle -> Texture[texture], 
  Lighting -> "Neutral", Boxed -> False, Axes -> False, 
  PlotPoints -> 100]

We have out this nice picture Figure 1

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

    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]}];
{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 Figure 2

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.8174971529112494627750732`3.9999565727231374*^6/3600, 
  4.392608409884539573243062`3.9999565727231374*^6/3600, 
  5.23208468377358490763874348`3.9999565727231374*^7/24/3600, 
  2.24299376347513480729423521`3.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

Figure 3

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

Figure 4

$\endgroup$
13
  • 1
    $\begingroup$ Version 12.1 It seems to work. Start a new kernel and try again. $\endgroup$ Dec 20 '20 at 11:46
  • 1
    $\begingroup$ Awesome example, my thanks and +1! $\endgroup$
    – Hans Olo
    Dec 20 '20 at 13:24
  • 2
    $\begingroup$ 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. $\endgroup$ Dec 20 '20 at 13:27
  • 2
    $\begingroup$ @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. $\endgroup$ Dec 20 '20 at 14:19
  • 2
    $\begingroup$ 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. $\endgroup$
    – Bob Hanlon
    Dec 20 '20 at 21:20
3
$\begingroup$

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", 
       "Oberon"}, {"OrbitPeriod", "SemimajorAxis", "Radius"}], 
     "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;
radius = QuantityMagnitude[
   Entity["Planet", "Uranus"]["EquatorialRadius"], "Kilometers"];
radius1 = 
  QuantityMagnitude[Entity["Planet", "Uranus"]["PolarRadius"], 
   "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];
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]

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)

Figure 1

ringinnerradius = 34840;

ringouterradius = 50023; Sun = 
 Graphics3D[{White, 
   Sphere[{-40*ringouterradius*Cos[Xobliquity], 0, 
     40*ringouterradius*Sin[Xobliquity]}, .05*radius]}, 
  Boxed -> False];
Shadow = Graphics3D[{Black, 
    Cylinder[{{0, 0, 0}, {30*ringouterradius*Cos[Xobliquity], 
       0, -30*ringouterradius*Sin[Xobliquity]}}, radius1]}];
I0 = SetAlphaChannel[Import["https://i.stack.imgur.com/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, 
  ViewVector -> {{5 ringouterradius, 0, 10 radius}, {0, 0, 0}}, 
  PlotRange -> All, 
  Lighting -> {{"Ambient", GrayLevel[0]}, {"Point", 
     White, {0, 0, 30*ringouterradius}}}]

M1 = QuantityMagnitude[Entity["Planet", "Uranus"]["Mass"]];
M2 = QuantityMagnitude[Entity["Star", "Sun"]["Mass"]];
Ra = QuantityMagnitude[
    Entity["Planet", "Uranus"]["EquatorialRadius"]]*10^3;
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]

Figure 2

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; 
radius = QuantityMagnitude[Entity["Planet", "Neptune"]["EquatorialRadius"], "Kilometers"]; 
radius1 = QuantityMagnitude[Entity["Planet", "Neptune"]["PolarRadius"], "Kilometers"]; 
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]; 
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]

Image rim for rings visualization Figure 3

ringinnerradius = 40900;

ringouterradius = 62302; Sun = 
 Graphics3D[{White, 
   Sphere[{-40*ringouterradius*Cos[Xobliquity], 0, 
     40*ringouterradius*Sin[Xobliquity]}, .05*radius]}, 
  Boxed -> False];
Shadow = Graphics3D[{Black, 
    Cylinder[{{0, 0, 0}, {30*ringouterradius*Cos[Xobliquity], 
       0, -30*ringouterradius*Sin[Xobliquity]}}, radius1]}];
I0 = SetAlphaChannel[Import["https://i.stack.imgur.com/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, 
  ViewVector -> {{5 ringouterradius, 0, 10 radius}, {0, 0, 0}}, 
  PlotRange -> All, 
  Lighting -> {{"Ambient", GrayLevel[0]}, {"Point", 
     White, {-40*ringouterradius*Cos[Xobliquity], 0, 
      40*ringouterradius*Sin[Xobliquity]}}}]

M1 = QuantityMagnitude[Entity["Planet", "Neptune"]["Mass"]];
M2 = QuantityMagnitude[Entity["Star", "Sun"]["Mass"]];
Ra = QuantityMagnitude[
    Entity["Planet", "Neptune"]["EquatorialRadius"]]*10^3;
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]

Figure 4

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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