Creating a simulation of our Solar System

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"]]],
      (*radius*)metersToAU[AstronomicalData[object, "Radius"]],
      (*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] :=
  Module[{orbit, radius, color, exponent, position},
   {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]}],
    (*object*)Sphere[position, scale^exponent*radius]}];

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"]],
   metersToAU[AstronomicalData[object, "Radius"]],
   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] := 
  Module[{orbit, position, radius, color, exponent},
   {orbit, position, radius, color, exponent} = objectData[object];

   {If[showDisk, {Opacity[.02], Polygon[orbit]}],
    Tooltip[{Opacity[.1], color, Line[orbit], Opacity[1],
      Sphere[position, scale^exponent*radius]}, tooltipData[object]]}];

(* 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;

(* preload tooltipData *)
tooltipData /@ planets;

Manipulate[
 Graphics3D[
  graphics[#, scale, showDisk] & /@ visiblePlanets,
  Boxed -> False, Background -> Black, SphericalRegion -> True, 
  PlotRangePadding -> .5,
  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"}]

Regarding your first point:

The planetary orbital planes are indeed inclined to that of the Earth and you have the inclinations. However, you also need to know the azimuthal locations of the ascending and descending nodes which together define the line about which to pivot the orbital ellipse.

Regarding your second point:

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.

See also elliptic orbit.