Planetary Orbits

Mathematica has some basic planetary data built in.  We can use it to re-discover Kepler’s third law, and, from that, infer the form of Newton’s law of Universal Gravitation.  Since the orbits are not quite circular, the appropriate quantity to use is the semi-major axis (i.e. half of the major axis of the orbital ellipse.)

Clear["Global`*"];

planets = {"Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune"};

data = Table[{PlanetData[p, "MajorAxis"]/2.0, 
PlanetData[p, "OrbitPeriod"]}, {p, planets}];

Pluto is technically classified as a “Dwarf Planet”, so we have to add it in separately.

AppendTo[planets, "Pluto"];AppendTo[data, {MinorPlanetData["Pluto", "MajorAxis"]/2.0, MinorPlanetData["Pluto", "OrbitPeriod"]}];

TableForm[data,
TableHeadings->{planets, {"Orbital Radius", "Orbital Period"}}]

For plotting and computation, we should switch all data to uniform units of meters and seconds.

data = UnitConvert[data, "SIBase"]; 
TableForm[data,
TableHeadings->{planets, {"Orbital Radius", "Orbital Period"}}]

Animation (inner planets only)

For animation, it is convenient to scale all lengths by Earth’s orbital radius, and all times by Earth’s orbital period (1 year):

scaled = Table[{data[[i, 1]] / data[[3, 1]], data[[i, 2]]/data[[3, 2]]}, {i, 1, Length[data]}];

x[r_, T_, t_] := r Cos[2 π t/T]
y[r_, T_, t_] := r Sin[2 π t/T]
pos[r_, T_, t_] := {x[r, T, t], y[r, T, t]}

planets = {1, 2, 3, 4};   (* Just show the 4 inner planets *)
Animate[Show[{
Graphics[{Yellow, Disk[{0, 0}, 0.1]}] (* Sun *),
ParametricPlot[ Table[pos[scaled[[i, 1]], 1, ϕ], {i, planets}], {ϕ, 0, 1}, PlotRange -> All],
ListPlot[Table[{pos[scaled[[i, 1]], scaled[[i, 2]], t]}, {i, planets}], 
PlotStyle->{{Brown, PointSize[0.04]}, {Gray, PointSize[0.05]}, {Blue, PointSize[0.06]}, {Red, PointSize[0.04]}
(*, {Brown, PointSize[0.10]}, {Yellow, PointSize[0.08]} *)}]
}, ImageSize->Scaled[0.6]
], {t, 0,5, 0.01, DefaultDuration -> 60, AnimationRunning-> False}]

opts = {AspectRatio -> 1/2,LabelStyle -> Larger, Frame -> True, FrameLabel -> {"Mean Radius (m)", "Orbital Period (s)"}, ImageSize -> Scaled[0.75]}; (* Common plot options *)

ListPlot[data, opts]

This plot is not terribly useful for visualization, since all of the inner planets are squished down in the lower left corner, and it doesn’t follow any obvious curve.

Power-Laws and Log-Log Plots

The T vs. R  plot is not terribly useful for visualization, since all of the inner planets are squished down in the lower left corner.  Using a Log-Log plot spreads things out  nicely.

ListLogLogPlot[data, opts]

This plot shows a suggestive linear trend.  Try fitting it to find the slope:

logdata = Log[10, QuantityMagnitude[data]]

llfit[x_] = Fit[logdata, {1, x}, x]

The slope is very nearly 1.5.

Show[{ListPlot[logdata], Plot[llfit[x], {x, logdata[[1, 1]], logdata[[-1, 1]]}, PlotStyle -> Red, PlotLegends -> "Slope=1.49983"]}, LabelStyle -> Larger, Frame -> True, FrameLabel -> {"Log(Mean Radius)", "Log(Orbital Period)"}, ImageSize -> Scaled[0.75], AspectRatio -> 1/2
]

Log-Log Plots and Power Laws

Suppose the period T is related to the orbital radius R by a power law   T = c (R^m).  How would that show up in a plot?  Try taking logs of both sides:
	T = c (R^m)
	log(T) = log(c) + log((R^m)) = log(c) + m log(R)
This suggests that if you plot log(T) vs. log(R), you will get a straight line with slope ‘m’.  In this case, we have m = 1.5.  This means
	T = c (R^1.5) = c (R^(3/2))
squaring both sides gives
	(T^2) = c (R^3)
which is Kepler’s 3rd law.

Kepler’s Third Law -- Log Log Plot

The Log-Log analysis suggests that another way to plot it is to plot (T^2) vs. (R^3).  It shows a linear trend, so we can fit it with a straight line.  We will see that the slope is given by ((4 (π^2))/(G  (Subscript[M, sun])))

data32 = Table[{((QuantityMagnitude[data[[i, 1]]])^3), ((QuantityMagnitude[data[[i, 2]]])^2)}, {i, 1, Length[data]}];

fit =LinearModelFit[data32, { x}, x,IncludeConstantBasis->False ];

Show[{ListPlot[data32], Plot[fit[x], {x, data32[[1, 1]], data32[[-1, 1]]}, PlotStyle -> Red]}, LabelStyle -> Larger, Frame -> True, FrameLabel -> {"(Mean Radius (m))^3", "(Orbital Period (s))^2"}, ImageSize -> Scaled[0.75], AspectRatio -> 1/2
]

slope = fit'[x]

G = Quantity[1, "GravitationalConstant"];
Msun = Quantity[1, "SolarMass"];
UnitConvert[4 (π^2)/(G * Msun),"SIBase"]

Kepler’s Third Law

Ttheory[R_, GM_] := √(((4 π^2 R^3)/GM))

fit = NonlinearModelFit[QuantityMagnitude[data], Ttheory[R, GM], {{GM,10^(20)}}, R]

G = Quantity[1, "GravitationalConstant"];
Msun = Quantity[1, "SolarMass"];
UnitConvert[4 (π^2)/(G * Msun),"SIBase"]

Show[{
ListPlot[data, PlotStyle->Red, LabelStyle->Larger, AxesLabel -> {"Mean Radius (m)", "Period (s)"}, ImageSize -> Scaled[0.6]],
Plot[fit[r], {r, 0, Max[QuantityMagnitude[data[[All, 1]]]]}]
}]

fit["BestFitParameters"]

UnitConvert[G*Msun, "SIBase"]