---
title: Pluto's Closest Approach
author:
  - Jed Rembold
date: Jan 31, 2025
lightbox: true
format:
  html:
    theme: cerulean
    embed-resources: true
    title-block-banner: true
---
In this problem we have been tasked with determining the closest approach the minor planet Pluto makes with the Sun. As a starting point, we have been given a CSV which contains X, Y and Z coordinates of Pluto ranging from the year 1800 until nearly 2100. All coordinate positions are given in units of astronomical units, or AU.

We'll be using Python for analysis in this problem, and thus will start things off by importing the libraries that we will be using, as well as applying my preferred default plot style. Here we are using the Pandas library to facilitate reading in and working with data series. We'll be using the Matplotlib library for visualization purposes, and the Numpy library to give access to some mathematical functions that apply to entire lists.

```{python}
import pandas as pd
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-darkgrid')
import numpy as np
```

Importing in our data, we can take a quick look at what we are working with. Here we see two columns corresponding to the date of the observation, as both a Julian Date and a more normal calendar date. Then we see the X, Y, and Z position of the planet in units of AU. 

```{python}
pluto = pd.read_csv('~/Teaching/Data_in_Cosmos/Materials/Lectures/demos/pluto_positions.csv')
pluto.head()
```

Ideally, one of the coordinates would be much smaller than the others at all times. This is because we need to fit a 2D ellipse to this data, and therefore it would be preferable if the ellipse that the data encodes was contained in purely 2 of the 3 axes. In this case, if we look at some descriptive statistics, we see that the Z coordinate is, on average, about an order of magnitude (10x) smaller than the other two coordinates. As such, we will assume that our ellipse lies mostly in the X-Y plane, even if that isn't perfectly true. 
```{python}
pluto[['JDTDB', 'X', 'Y', 'Z']].mean()
```

To give us some visual to work with, and to confirm that the provided information fully describes the arc of an ellipse, we go ahead and plot the X coordinates vs the Y coordinates in @fig-pluto_orbit below.
```{python}
#| code-fold: true
#| fig-align: center
#| fig-cap: The orbit of Pluto projected onto the X and Y axes as seen from above the Z axis. The elliptical nature of the orbit is clear.
#| label: fig-pluto_orbit
plt.subplots(figsize=(8,8))
plt.plot(pluto.X, pluto.Y, 'o', markerfacecolor='white')
plt.gca().set_aspect('equal')
plt.xlabel('X (AU)')
plt.ylabel('Y (AU)')
plt.title('Orbit of Pluto from 1800 to ~2100')
plt.show()
```

That looks pretty elliptical, and it has very full coverage around the entire edge, so our fitting algorithm will hopefully work well. We need to determine, in particular, the semi-major and semi-minor axes, as they will be necessary to compute the final perihelion distance.

Here we will utilize the provided ellipse fitting library to extract the desired parameters, starting with importing the desired library functions.
```{python}
from fitellipse2 import fit_ellipse, get_ellipse
```

And then we'll apply the `fit_ellipse` function to our X and Y data:
```{python}
efit = fit_ellipse(pluto.X, pluto.Y)
efit
```
We should never just blindly trust the results of a fit, and so as a sanity check we'll plot the fit back over the original data. This revealed that our fit angle was 90 degrees off, so we adjusted as necessary, as can be seen in @fig-fitcheck.
```{python}
#| code-fold: true
#| fig-align: center
#| fig-cap: Fit ellipse overlaid atop original data. While the shapes match nicely, the original fit was off in its rotation by 90 degrees, which occassionally happens. Fixing this results in the excellent agreement between the data and the fit ellipse.
#| label: fig-fitcheck
plt.subplots(figsize=(8,8))
plt.plot(pluto.X, pluto.Y, 'o', markerfacecolor='white', label='Original data')
plt.plot(*get_ellipse(efit), '--', label="Original Fit Results")
#efit['angle'] += np.pi / 2
plt.plot(*get_ellipse(efit), '--', color='C3', label="Adjusted Fit Results")
plt.gca().set_aspect('equal')
plt.xlabel('X (AU)')
plt.ylabel('Y (AU)')
plt.title('Orbit of Pluto from 1800 to ~2100')
plt.legend()
plt.show()
```


Confident that we have a good fit to the data, we can extract out the important ellipse properties from the fit parameters, resulting in a semi-major axis of 
```{python}
#| code-fold: true
f"{efit['semimajor']} AU"
```
and a semi-minor axis of
```{python}
#| code-fold: true
f"{efit['semiminor']} AU"
```
We might note that despite the orbit being clearly elliptical, these values aren't all that different! To use these values now to try to determine our perihelion distance, it can help to look at @fig-ellipse_geom:

![Image showing the geometrical relationships between the semi-major axis, the foci distance, and the perihelion distance](./perihelion.svg){width=50% height=50% fig-align='center' #fig-ellipse_geom }

So the semi-major axis minus the foci distance $f$ will get us the distance that we want. We know that the foci distance $f$ is given in terms of the semi-major and minor axes via:
$$ f = \sqrt{a^2 - b^2} $$
where $a$ is the semi-major axis and $b$ the semi-minor axis. We have both of those values now, and can thus calculate the foci distance.
```{python}
f = np.sqrt((efit['semimajor']) ** 2 - (efit['semiminor']) ** 2)
f"Foci distance = {f} AU"
```

By the above picture, we have that the perihelion distance would be the semi-major axis minus the foci distance. Thus we get:
```{python}
perihelion = efit['semimajor'] - f
f"Perihelion distance = {perihelion} AU"
```
So the closest that Pluto ever gets to the Sun is still over 21 times further away than Earth is from the Sun.

:::{.callout-warning}
My gut was that these coordinates are Sun-centered, so that the Sun should be at the origin.
```{python}
np.sqrt(pluto.X**2 + pluto.Y**2).min()
```
As seen below, that value matches a LOT better with the accepted perihelion distance. But how did we get off? Plotting shows us that our predicted position is suspiciously twice as far from the center as expected.
```{python}
#| fig-align: center
#| code-fold: true
plt.plot(pluto.X, pluto.Y, '.', label='Data')
plt.plot(0,0,'.', label='Origin')
uv = np.array([np.cos(efit['angle']), np.sin(efit['angle'])])
loc = efit['center'] - uv * f
plt.plot(*efit['center'], '.', label='Center')
plt.plot(*loc, '.', label='Calc Foci')
plt.gca().set_aspect('equal')
plt.legend(loc='upper right', fancybox=True, frameon=True, framealpha=1)
plt.show()
```
I'm still not sure why this is happening, or if it is just a wild coincidence. I should perhaps test it with the orbit of Mercury, which is much more planar?
:::

### Appendix
As an aside, the perihelion distance reported on Wikipedia is 29.658 AU, which makes for a fairly substantial 28% error in our estimate. Most likely, this is owing to Pluto's more extreme than usual tilt relative to the ecliptic plane. We might be able to see this if we plot our coordinates instead in 3 dimensions:
```{python}
#| code-fold: true
from mpl_toolkits.mplot3d import Axes3D

with plt.style.context(('default')):
  plt.figure(figsize=(10,6))
  ax = plt.subplot(111, projection='3d')
  ax.plot(pluto.X, pluto.Y, pluto.Z)
  ax.set_xlabel('X (AU)')
  ax.set_ylabel('Y (AU)')
  ax.set_zlabel('Z (AU)')
  ax.set_aspect('equal')

  ax.view_init(elev=20, azim=45)
  plt.show()
```
Here we can see that yes, we _definitely_ have a bit of tilt still, which we aren't accounting for in our about fitting. That means our semi-major axis is measured to be slightly smaller than it should be, which then results in the lower than expected perihelion distance.