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

Update June 2015

Here is an updated version of the program. I've made it compatible with newer Mathematica versions (AstronomicalData returns Quantity structures in newer versions, which wrangled calculations). It should now work on versions 8 through 10. Let me know if it doesn't.

I added animation and simplified the presentation (no tooltips in the Graphics3D) and also added a star field in the background. I was aiming for reality-mapped stars but StarData coughed and sputtered too much to be usable, and the public star API that I tried gave only partial data. It's possible I misused these tools, but either way it's an improvement to be made.

The program loads the planet colors from WolframAlpha into a memoization table, so it takes a few seconds to run the first time.

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

objects = Prepend[AstronomicalData["Planet"], "Sun"];

(*ClearAll[dataMemo];*)
dataMemo[object_] := (
dataMemo[object] = {
(*orbit*)If[object === "Sun", {{0, 0, 0}}, First[AstronomicalData[object, "OrbitPath"]]],
(*color*)First[Cases[WolframAlpha[object <> " color"], _RGBColor, Infinity]],
(*drawing scale*)If[object === "Sun", 1, 2.6]} /. Quantity[a_, _] :> a
);

objectData[object_, dateOffset_] :=
Append[dataMemo[object],
(*position*)metersToAU[AstronomicalData[object, {"Position", DatePlus[dateOffset]}]] /. Quantity[a_, _] :> a
];

dataMemo /@ objects;

createGraphics[object_, scale_, dateOffset_, showOrbit_: True, showDisk_: False] :=
{orbit, radius, color, exponent, position} = objectData[object, dateOffset];

{color, Glow[color],
(*orbit disk*)If[showDisk, {Opacity[.1], EdgeForm[None], Polygon[orbit]}],
(*orbit*)If[showOrbit, {Opacity[.1], Line[orbit]}],

setterBar = # -> Tooltip[ImageResize[AstronomicalData[#, "Image"], {30, 30}], #] & /@ objects;

(*random star distribution*)
ratio = 1.91;(*8Volume[Cuboid[]]/Volume[Ball[]]*)
{numstars, starDistance} = {10000, 50};
stars = Normalize /@ Select[RandomReal[{-1, 1}, {Floor[numstars*ratio], 3}], Norm[#] < 1 &];

daysInFuture = 0;
Print[DynamicWrapper[Dynamic[DatePlus[daysInFuture]],
Refresh[daysInFuture++, UpdateInterval -> .01, TrackedSymbols :> {}]]];

Module[{viewvector = {0, 8, 0}},
Manipulate[
viewvector[[2]] = zoom;

Graphics3D[{
{White, Opacity[.4], AbsolutePointSize[2], Point[starDistance*stars]},
{Specularity[.1], Dynamic[
createGraphics[#, scale, daysInFuture, showOrbits, showDisk] & /@
visiblePlanets]}},

Boxed -> False, Background -> Black, SphericalRegion -> True,
ImageSize -> Large, PlotRangePadding -> .5, ViewVector -> (*Dynamic[*)viewvector(*]*),
Lighting -> {{"Ambient", White}, {"Directional", White, ImageScaled[{0, 0, 1}]}}],

{{visiblePlanets, Take[objects, 5], "planets"}, setterBar, ControlType -> TogglerBar},
{{showDisk, False, "show planes"}, {False, True}},
{{showOrbits, True, "show orbits"}, {True, False}},
{{zoom, 8}, 2, 30},
{{scale, 13.5, "scale"}, 10, 20, Appearance -> "Labeled", ControlType -> None}]]


Previous version

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"}]

-
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.