# Creating a simulation of our Solar System

As the title says, my objective is to create a simulation of the motion of the planets in our Solar System using Mathematica. All the theoretical background regarding the equations of motion of planets is well know so, the difficult part is to create a functional code. Below I present what I have so far

Clear["Global"];
pSun = {0, 0};
rMercury = a/(1 + e*Cos[θ]) /. {a -> 0.387, e -> 0.2056, i -> 7.005};
rVenus = a/(1 + e*Cos[θ]) /. {a -> 0.723, e -> 0.0068, i -> 3.3947};
rEarth = a/(1 + e*Cos[θ]) /. {a -> 1, e -> 0.0167, i -> 0};
rMars = a/(1 + e*Cos[θ]) /. {a -> 1.524, e -> 0.0934, i -> 1.851};
rJupiter = a/(1 + e*Cos[θ]) /. {a -> 5.203, e -> 0.0484, i -> 1.305};
rSaturn = a/(1 + e*Cos[θ]) /. {a -> 9.537, e -> 0.0542, i -> 2.484};
rUranus = a/(1 + e*Cos[θ]) /. {a -> 19.191, e -> 0.0472, i -> 0.770};
rNeptune = a/(1 + e*Cos[θ]) /. {a -> 30.069, e -> 0.0086, i -> 1.769};
rPluto = a/(1 + e*Cos[θ]) /. {a -> 39.482, e -> 0.2488, i -> 17.142};

p0 = ListPlot[{pSun}, Axes -> False,
PlotStyle -> {RGBColor[1, 0.65, 0], PointSize[0.035]}];
p1 = PolarPlot[rMercury, {θ, 0, 2 π}, PlotStyle -> Gray];
p2 = PolarPlot[rVenus, {θ, 0, 2 π}, PlotStyle -> Orange];
p3 = PolarPlot[rEarth, {θ, 0, 2 π}, PlotStyle -> Blue];
p4 = PolarPlot[rMars, {θ, 0, 2 π}, PlotStyle -> Red];
p5 = PolarPlot[rJupiter, {θ, 0, 2 π}, PlotStyle -> Brown];
p6 = PolarPlot[rSaturn, {θ, 0, 2 π}, PlotStyle -> Magenta];
p7 = PolarPlot[rUranus, {θ, 0, 2 π}, PlotStyle -> Cyan];
p8 = PolarPlot[rNeptune, {θ, 0, 2 π}, PlotStyle -> Darker[Green]];
p9 = PolarPlot[rPluto, {θ, 0, 2 π}, PlotStyle -> Black];

S1 = Show[{p1, p2, p3, p4, p0}, Axes -> False, Frame -> True,
FrameTicks -> None, PlotRange -> All, AspectRatio -> 1,
ImageSize -> 500]
S2 = Show[{p5, p6, p7, p8, p9, p0}, Frame -> True, FrameTicks -> None,
Axes -> False, PlotRange -> All, AspectRatio -> 1,
ImageSize -> 500]


-------------------- INNER SOLAR SYSTEM --------------------

-------------------- OUTER SOLAR SYSTEM --------------------

I have divided the Solar System in two different plots. The first one contains the Earth-type planets (Mercury, Venus, Earth and Mars) and the second one the gas giants. The orbits of the planets are ellipses and are given in polar form using Kepler's theory. Obviously, Sun is stationary at one of the focuses. The motion of all the planets is two-dimensional. However, all orbits are not co-planar. Here comes the first issue:

(1). Somehow, the ellipses must be rotated according to the inclination of each planet. The inclination (i) in degrees is given. So, we should have a 3D box containing all the 2D inclined ellipses. Earth's inclination is zero thus defying the primary plane (ecliptic) from which we measure inclination.

(2). At every orbit, it would be nice if there was a color dot (like the Sun I already have) indicating each planet. Since we speak of a simulation, every dot (planet) should circulate around Sun following the corresponding orbit. Here we have a problem. Every planet has each one rotational velocity according to its mass. However, the polar equation giving the orbit does not include the mass of the planet. Any suggestions here would be greatly appreciated.

There are also few additional minor issues. For the time being, the first two issues are important and should be resolved first.

From my point of view, this task is indeed not only interesting but also very challenging. Everyone should have a nice model of our Solar System!

-
I think you're working too hard. Create a simulation of a single planet, driven by minimal input: planet's major and minor axes, inclination, rotational speed, orbital parameters, and starting location would do. Then apply the simulation function to an array of such data. You might even contemplate creating a simulation from first principles, using a central gravitational force, initial position and velocity, and angular velocity of rotation, and let Mathematica dynamically solve the equations of motion. Such an approach could be extended to accommodate gravitational interactions. –  whuber Feb 7 '13 at 22:58
Maybe this will give some ideas: Demonstrations - you can inspect the source code there. –  Jens Feb 7 '13 at 23:09
you might find this useful. Did you check Solar System Dynamics by Carl D. Murray and Stanley F. Dermott. Mathematica nb, movies and some explanations available. –  s.s.o Feb 8 '13 at 2:57

I'm not sure what your goal is, exactly, but here is a simulation I cooked up. It should give you some ideas:

metersToAU[m_] := m/(1.496*10^11) ;
orbit = First@AstronomicalData["Earth", "OrbitPath"];
earthCurrentPosition = AstronomicalData["Earth", "Position"] // metersToAU;

(* astronomical data for an object. memoized *)
allData[object_] := allData[object] = Grid[Select[
{#, AstronomicalData[object, #]} & /@ AstronomicalData["Properties"],
Head[#[[2]]] =!= Missing && #[[1]] =!= "OrbitPath" && #[[1]] =!= "Image" &], Alignment -> Left];

(* lets us smoothly animate position of Earth *)
interp = Interpolation[Table[{i, orbit[[i]]}, {i, 1, Length[orbit]}]];

(* from AstronomicalData Help page *)
ColorData["BlackBodySpectrum"][
AstronomicalData[star, "EffectiveTemperature"]],

Manipulate[
Graphics3D[{
Tooltip[
Column[{AstronomicalData["Sun", "Image"], allData["Sun"]}]],

{Opacity[.1], Pink, Line[orbit]},

Tooltip[
{Lighter[Gray, .8], Text["Earth's current position", {0, 0, .5} + earthCurrentPosition],
Red, Arrow[{{0, 0, .35} + earthCurrentPosition, {0, 0, .05} + earthCurrentPosition}],
Column[{AstronomicalData["Earth", "Image"], allData["Earth"]}]]}

, Boxed -> False, Background -> Black, SphericalRegion -> True,

(* this lighting specification is from the original starColorPlot *)
Lighting -> {{"Ambient", Gray}, {"Directional", White, ImageScaled[{0, 0, 1}]}}],

{t, Length[orbit], 1},
{{scale, 10, "radius scale"}, 1, 20, 1, Appearance -> "Labeled"}]


The key here, which should not be missed, is that I'm astoundingly lazy, so I just pulled Earth's orbit out of AstronomicalData. You can in fact pull all sorts of things out, including the radiuses of the planets.

The most interesting aspect is that you can query the location of a planet on a specified date. In this little simulation, I'm just spinning the Earth around based on t, but you can parameterize it based on the date. And you can parameterize the locations of all of the planets based on the date, which would let you have a scientifically-accurate "simulation" without having to account for any celestial mechanics yourself. I'm not sure if you want all your work done for you, but I don't mind. :D

Notice that the Sun and Earth have tooltips.

And if you're bored enough, you can use Texture + SphericalPlot3D (look at the "Applications" section) to draw some detail. You would use Show to combine the spherical plots. It makes the animation a lot slower, though.

Regarding your issues with perspectives, you can specify views on a Graphics3D with some options (ViewAngle and others). And you can show the same Graphics3D from different perspectives by doing something like:

g = Graphics3D[...];
Grid[{{Show[g, ViewAngle -> blah], Show[{g, ovals}, ViewAngle -> different blah]}}]


So for example you could have one of the views that shows ovals in the orbit planes of the planets. You can use a flat Cylinder to draw a disk in Graphics3D, and Scale to make an ellipse. Update: Actually you can just do Polygon on the orbit. Durrr. In the above code, use {Opacity[.1], Pink, Line[orbit], Polygon[orbit]} to get:

I should mention that this isn't as complicated/difficult as it might look. This simulation took me about 15 minutes to make. So don't be intimidated if it looks "complicated." It only looks that way because of all the {s and [s.

Side note: Remember the power of =:

## Update

Version that handles all planets.

metersToAU[m_] := m/(1.496*10^11);

(* these colors aren't available from Astronomical data,
but I got them through the free-form input. "=mars color", "=earth color", etc *)
colorTable = {
"Sun" -> Blend["BlackBodySpectrum", AstronomicalData["Sun", "EffectiveTemperature"]],
"Mercury" -> RGBColor[0.598209, 0.577666, 0.576307],
"Venus" -> RGBColor[0.753588, 0.740667, 0.706174],
"Earth" -> RGBColor[0.598209, 0.577666, 0.576307],
"Mars" -> RGBColor[0.591699, 0.37999, 0.19484],
"Jupiter" -> RGBColor[0.757233, 0.697683, 0.666215],
"Saturn" -> RGBColor[0.767425, 0.699073, 0.563738],
"Uranus" -> RGBColor[0.574328, 0.751126, 0.827463],
"Neptune" -> RGBColor[0.556615, 0.747549, 0.88086]
};

(* returns {orbit, current position, radius, color, scale exponent} *)
objectData[object_] := {
If[object === "Sun", {{0, 0, 0}}, First[AstronomicalData[object, "OrbitPath"]]],
metersToAU[AstronomicalData[object, "Position"]],
object /. colorTable,
If[object === "Sun", 1.1, 2.2]
};

(* all astronomical data for an object. memoized *)
Clear[tooltipData];
tooltipData[object_] := tooltipData[object] = Column[{
AstronomicalData[object, "Image"], Style[object, Bold],
Grid[
Select[{#, AstronomicalData[object, #]} & /@
AstronomicalData["Properties"],
Head[#[[2]]] =!= Missing && #[[1]] =!= "OrbitPath" && #[[1]] =!=
"Image" &],
Alignment -> Left]}];

(* create the Graphics3D primitives for a given object *)
graphics[object_, scale_, showDisk_: False] :=
{orbit, position, radius, color, exponent} = objectData[object];

{If[showDisk, {Opacity[.02], Polygon[orbit]}],
Tooltip[{Opacity[.1], color, Line[orbit], Opacity[1],

(* type AstronomicalData["Planet"], press Ctrl+. twice, then press Ctrl+Shift+Enter *)
planets = {"Sun", "Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune"};
setterBar = # -> Tooltip[ImageResize[AstronomicalData[#, "Image"], {30, 30}], #, TooltipDelay -> .1] & /@ planets;

tooltipData /@ planets;

Manipulate[
Graphics3D[
graphics[#, scale, showDisk] & /@ visiblePlanets,
Boxed -> False, Background -> Black, SphericalRegion -> True,
Lighting -> {{"Ambient", White}, {"Directional", White, ImageScaled[{0, 0, 1}]}}],

{{visiblePlanets, Take[planets, 5], "planets"}, setterBar, ControlType -> TogglerBar},
{{showDisk, False, "show planes"}, {False, True}},
{{scale, 10, "scale"}, 1, 20, Appearance -> "Labeled"}]


This version, however, doesn't animate. You have to parameterize on the date when you query the position of the planets, which you should be able to figure out how to do somewhere in the Help docs for AstronomicalData.

-
That is a very good start! Indeed it is not so complicated. However, it's rather "long" code and only for one planet! Should I repeat it eight more times for the rest planets? Since we exploit AstronomicalData wouldn't be easier to insert somehow the well-known inclinations of the planets in order to rotate each plane? –  Vaggelis_Z Feb 8 '13 at 8:29
Also why do you Manipulate Sun's radius scale? Sun has a fixed radius, at least for the next 4 billion years! –  Vaggelis_Z Feb 8 '13 at 8:31
I tried to insert one more planet (Mercury) but the simulation is not working. Could you please add one more planet? Every planet and its orbit should be in different color. –  Vaggelis_Z Feb 8 '13 at 8:32
@Vaggelis_Z 1) The path returned by AstronomicalData is in 3D, so it's already inclined. The inclinations are very small and almost indistinguishable. 2) I Manipulate the radius of all the objects, not just the Sun. Remember that the planets are separated by ridiculously large distances, so if you draw the radiuses to real-world scale, most of the planets are just dots. 3) Yes, that code requires a lot of refactoring. You don't need to do it for each individual planet. You just do it for one planet and parameterize on its name, which is made easy by AstronomicalData. (See the update I made). –  amr Feb 9 '13 at 4:15
Let me know if you still have questions/issues. –  amr Feb 9 '13 at 4:17

If you assume, as you say, that the "sun is stationary at one of the [foci]" then you are assuming its mass is much larger than any planet. That is, the planet and Sun do NOT orbit about a common centre of mass, rather the planet orbits about the centre of the Sun. In this case of much larger solar mass, the orbital velocity is independent of planetary mass and given by the vis viva equation $v^2=G M_{\unicode{x2609}} \left(\tfrac2r - \tfrac1a\right)$, where $v$ is the velocity, $G$ is the gravitational constant, $M_{\unicode{x2609}}$ is the solar mass, $r$ is the current distance of the planet from the Sun, and $a$ is the major semi-axis.