Kepler and the Orbit of Mars

To participate in the interactive aspects of this notebook, select the ‘Open in Colab’ option!

If you’re running this in Colab, make sure to save a copy of the notebook in Google Drive to save your changes. You will need to download the ‘triangulation.csv’ dataset from the ‘data’ folder and change the path where it is called.

In this notebook, we will use modern techniques and programming to emulate, visualise and validate Kepler’s work, demonstrating the similarities between his methods and ‘data science’ as we know it today.

This notebook contains some ‘helper’ functions which are hidden for clarity. The full code is visible in Colab.

Triangulation

Through the knowledge of Mars’ 687-day orbital period, Kepler was able to plot the positions of Mars in its orbit with a triangulation technique. He used pairs of Brahe’s datapoints, each separated by an integer multiple of 687 days, to triangulate the position of Mars and track its orbit by examining the differences in its apparent position in the sky. This is possible due to the fact that Mars would be in the same place relative to the Sun each time.

[Source: http://galileo.phys.virginia.edu/classes/609.ral5q.fall04/LecturePDF/L08-KEPLER.pdf]

Visualising Earth’s Locations

The following code uses a small dataset of Mars’ locations to mimic Kepler’s work and create a visual reconstruction of Earth’s orbit.

[ ]:

# If you're running this notebook, uncomment the code in this cell to install the required packages.
# ! pip install csv
# ! pip install math
# ! pip install matplotlib
# ! pip install numpy
# ! pip install mpl_toolkits

[ ]:

earthLocations = None
marsAngles = None

earthLocations, marsAngles = loadData1()

plotEarthLocations(earthLocations)

Finding the five positions of Mars

Next, pairs of Earth’s positions separated by multiples of 687 days and their corresponding angles relative to Mars will be used to triangulate the locations of Mars, which will be recorded and displayed below.

[ ]:

marsLocations = []

for i in range(5):
    index1 = 2 * i
    index2 = 2 * i + 1
    marsLocation = findMars(earthLocations[index1], marsAngles[index1],
                                        earthLocations[index2], marsAngles[index2])
    marsLocations.append(marsLocation)

print(marsLocations)

Fitting a circle to the Orbit

Having recorded the positions of Mars, the natural progression for Kepler’s method was to attempt to fix these to a circular orbit, based on the Copernican model of the solar system to which he subscribed at the time.

In the following visualisation, it is clear that the data points do not quite follow the trajectory of a circle.

[ ]:

triangulatedRadius = None

triangulatedRadius = computeRadius(marsLocations)

plotTriangulations(marsLocations, triangulatedRadius)

Fitting the Mars Orbital Plane

So far, all of the data has been two-dimensional, involving projections of Mars onto the ecliptic plane, the orbital plane of Earth around the Sun. To find the actual locations of Mars in three-dimensional space, Kepler used ‘perihelic opposition’ (when Earth, the Sun, and Mars lie on a straight line) to deduce Mars’ orbital plane, as the angle to Mars is the same relative to both the Sun and Earth.

The longitudes to Mars with respect to the Sun and latitudes to Mars with respect to the Earth can be used to calculate the coordinates of Mars on the celestial sphere (an abstract sphere concentric to Earth with an arbitrarily large radius): geocentric latitudes are first used to calculate heliocentric latitudes, which together with longitudes can help us to locate Mars via points of intersection.

[ ]:

helioLong, geoLat = loadData2()

helioLat = []

helioLat = findHelioLat(triangulatedRadius, geoLat)

coordinates = []

coordinates = findCoordinates(helioLong, helioLat)

Fitting a Plane to the Orbit of Mars

From the coordinates of Mars on the celestial sphere, we can fit a plane to the orbit of Mars. First, we will find the parameters of this plane and then we can plot our results!

[ ]:

planeParameters = []

planeParameters = fitPlane(coordinates)

plotPlane(coordinates, planeParameters)

Fitting a Circle and Ellipse to the Orbit of Mars

Lifting Mars’ Projections onto the Orbital Plane

The following code takes our initial coordinates for the locations of Mars and ‘lifts’ them onto the orbital plane that we have found.

[ ]:

liftedMarsLocations = None

liftedMarsLocations = liftCoordinates(planeParameters, marsLocations)

print(liftedMarsLocations)

We now have the true locations of Mars on its orbital plane. Now let’s try to fit a circular orbit and an elliptical orbit, and compare the errors.

We should expect to obtain a lower error for an elliptical fit, consistent with Kepler’s results.

Fitting a Circular Orbit

The fllowing code plots a ‘circle of best fit’ to the orbit of Mars using the distances between the locations of Mars and the Sun to calculate a mean radius. The total error is then calculated using the individual errors of each location and the circle.

[ ]:

circleRadius = None
circleLoss = None

circleRadius, circleLoss = fitCircle(liftedMarsLocations)

print("Circle Radius:", circleRadius)
print("Circle Error (absolute):", circleLoss)

circlePercentageError = math.sqrt(circleLoss) / circleRadius * 100.0
print("Circle Error (percentage):", circlePercentageError)

Fitting an Elliptical Orbit

In a similar manner, the following code plots an ‘ellipse of best fit’ instead, calculating the error in the same way as before.

[ ]:

ellipseParameters = None
ellipseLoss = None

ellipseParameters, ellipseLoss = fitEllipse(liftedMarsLocations)

print("Ellipse Parameters:", ellipseParameters)
print("Ellipse Error (absolute):", ellipseLoss)

semimajorAxis = ellipseParameters[2] / 2.0
ellipsePercentageError = math.sqrt(ellipseLoss) / semimajorAxis * 100.0
print("Ellipse Error (percentage):", ellipsePercentageError)

Comparing the Circular and Elliptical Orbital Errors

As expected, we find that fitting a circle has approximately 16 times the error of fitting an ellipse. Therefore, the orbit of Mars must be elliptical! We can plot both models on one figure to visualise this:

[ ]:

plotBoth(liftedMarsLocations, circleRadius, ellipseParameters)

The plot clearly demonstrates that the ellipse passes through all of our datapoints within a reasonable error, whereas the circle is comparatively skewed. We have successfully validated and visualised Kepler’s results!

How could we analyse Brahe’s data today?

The evolution of data science and technology over the centuries since Kepler carried out his work means that, in today’s age, a neural network could be used to analyse Tycho Brahe’s observations, recognise patterns and relationships in the data and model an orbit. The advantage of this technique is that Kepler’s moment of inspiration to consider the possibility of an elliptical orbit would not have been necessary, not to mention that over 20 years of work would have been reduced to just a matter of minutes!

The following animation is a final example of how we can use modern programming languages to further study celestial orbits; we can investigate the initial conditions necessary for periodic or chaotic motion, as well as consider how massive objects can cause perturbations in orbits.

[Source: https://towardsdatascience.com/use-python-to-create-three-body-orbits-329ffb5b2627]

References

Orbits and Kepler’s Laws
Gentner, Dedre, et al. Analogical reasoning and conceptual change: A case study of Johannes Kepler. The journal of the learning sciences 6.1 (1997): 3-40.
Understanding What Kepler Did
Use Python to Create Three-Body Orbits
Kepler was a Data Scientist
Johannes Kepler and Tycho Brahe
Determination of the Orbit of Mars Using Kepler’s Triangulation Technique

The code for this section was adapted from material found at https://github.com/pulkitsingh/Mars-Orbit-Workshop.

If you wish to get an overview of the remaining topics in this course, click the button below.