---
title: "Kepler's Ellipses"
author: Jed Rembold
date: January 21, 2025
slideNumber: true
theme: tokyo-night-light
highlightjs-theme: tokyo-night-light
width: 1920
height: 1080
transition: slide
---


## Announcements
- Homework 1 due a week from Friday
	- Find time to meet with your partner this week to make progress!
	- You should be prepped to do all but the last part of problem 1 and probably problem 3 after today
- I still need about 10 of you to fill out the info and reflection form [here](https://docs.google.com/forms/d/e/1FAIpQLSdaUjeiwdLBnon0uamSm0lKLgxC88S24TsB70UcPs0_8yrzFQ/viewform?usp=sharing)


## A Note on Packages
- Given that all of you generally fall into 1 of 2 camps (R or Python), where possible I'm going to try to not rely on package specific implementations
- That said, there are some common ones that will repeatedly come up:
	- Python
		- _Pandas_ introduces dataframes, which makes ideas much more similar to R's dataframes
		- _Matplotlib_ is still the workhorse of data visualization in Python
		- _Numpy_ will probably make occasional appearances for pure number crunching
	- R
		- _Tidyverse_ brings so many nice quality of life improvements and added functionality
		- _ggplot_ is one of the best data visualization libraries I've seen
- Specific applied libraries may also spring up, but as much as possible I'll try to give both camps equivalent options

## Today's Plan
- Angular Differences and Distances
- What can we determine with positions?
	- Kepler's 1st Law
- Ellipse Math
- Ellipse Fitting

# Differences and Distances
## Angular Separation
::::::cols
::::col
- Frequently it is useful to describe the angular separation, or distance, between two points (as seen from Earth)
	- Pythagorean theorem will be close for small separations, but not quite accurate. Better to use:
	  $$ \theta \approx \sqrt{\left(\left(\alpha_A - \alpha_B\right)\cos\delta_A\right)^2 + \left(\delta_A - \delta_B\right)^2} $$
	  where $\alpha$ indicates Right Ascension and $\delta$ the Declination

::::

::::col
![](../images/ang_sep_diagram.svg)

::::
::::::


## Relating to Physical Distance
::::::cols
::::col
![](../images/physical_separation.svg)

::::

::::{.col style='font-size:.9em'}
- If you know both things to be approximately the same distance away, you can convert angular separation to a physical distance:
	$$ \frac{\theta}{360^\circ} = \frac{\text{separation}}{2\pi(\text{distance})} $$
	- $\text{distance}$ is the distance from the observer to the object in question, 
	- $\theta$ is the angular separation (in degrees)
	- $\text{separation}$ is the physical separation
- The physical separation will have the same units as the distance you use

::::
::::::


# Orbits
## Tycho Brahe
::::::cols
::::col
- Last real "naked-eye" astronomer
- Took very precise angle measurements of stars and planets
- Did not believe in a heliocentric universe, as he saw no stellar parallax
- Hired Johannes Kepler to analyze all his data
::::

::::col
![Brahe's observations](../images/brahe.jpg)
::::
::::::

## Johannes Kepler
- Kepler then spent **8 years** trying to reconcile Tycho's observations of Mars with the Ptolemaic (Earth centered) model

::::{.quote style='font-size:.9em'}
> Who would have thought it possible? This hypothesis, which so closely agrees with the observed oppositions, is nevertheless false? If I had believed that we could ignore those 8 minutes, I would have patched up my hypothesis accordingly. But since it was not possible to ignore them, those 8 minutes point the road to a complete reform of astronomy...

> Thou seest now, diligent reader, that the hypothesis based on this method not only satisfies the four positions on which it was based, but also correctly represents within 2 minutes all the other observations.

:::attribution
_Nova Astronomica_ -- Johannes Kepler
:::
::::


## Kepler's 1st Law
- The orbits of the planets are ellipses
- The Sun sits at one focus, with nothing at the other


\begin{tikzpicture}%%width=50%
[every node/.style={font=\sf, color=Black}]
  \def\a{4}
  \def\b{3}
  \draw (0,0) ellipse (\a cm and \b cm);
  \node[draw,fill,circle, inner sep=1pt,label={below:Nothing}] (f1) at ({-sqrt(\a*\a-\b*\b)},0) {};
  \node[fill,circle, inner sep=3pt, inner color=yellow, outer color=orange,label={below:Sun}] (f2) at ({sqrt(\a*\a-\b*\b)},0) {};
  \node at ($(0,0)+(45:{\a} and {\b})$) {\includegraphics[width=3mm]{world.png}};
\end{tikzpicture}


# Ellipse Math
## An Ellipses Aside
::::::{.cols style='align-items: flex-start;'}
::::col
- Look like squished circles

- Described by:
	- Major Axis (longest width)
		- Half lengths are "semi-"
	- Minor Axis (shortest width)
	- Foci distances
	- Eccentricity (squished-ness)
::::

:::::col

<!--![](../images/test_ellipses.svg)-->

::::rstack

:::{.fragment .current-visible .only-fragment}
\begin{tikzpicture}%%width=100%
[every node/.style={font=\sf, color=Black}]
  \draw[very thick] (0,0) ellipse (4cm and 3cm);
\end{tikzpicture}
:::

:::{.fragment .current-visible .only-fragment}
\begin{tikzpicture}%%width=100%
[every node/.style={font=\sf, color=Black}]
  \draw[very thick] (0,0) ellipse (4cm and 3cm);
  \draw[latex-latex,Cyan, thick] (-4,0) -- node[midway,above] {Major Axis}(4,0);
\end{tikzpicture}
:::

:::{.fragment .current-visible .only-fragment}
\begin{tikzpicture}%%width=100%
[every node/.style={font=\sf, color=Black}]
  \draw[very thick] (0,0) ellipse (4cm and 3cm);
  \draw[latex-latex, Red, thick] (0,0) -- node[midway,below] {Semi-Major Axis} (4,0);
\end{tikzpicture}
:::

:::{.fragment .current-visible .only-fragment}
\begin{tikzpicture}%%width=100%
[every node/.style={font=\sf, color=Black}]
  \draw[very thick] (0,0) ellipse (4cm and 3cm);
  \draw[latex-latex, Purple, thick] (0,0) -- node[midway,right, align=left] {Semi-Minor\\Axis} (0,3);
\end{tikzpicture}
:::

:::{.fragment .current-visible .only-fragment}
\begin{tikzpicture}%%width=100%
[every node/.style={font=\sf, color=Black}]
  \draw[very thick] (0,0) ellipse (4cm and 3cm);
  \fill ({sqrt(4*4-3*3)},0) node (f1) {} circle (3pt);
  \fill ({-sqrt(4*4-3*3)},0) node (f2) {} circle (3pt);
  \node (flab) at (0,0) {Foci};
  \draw[-latex] (flab) -- (f1);
  \draw[-latex] (flab) -- (f2);
  \coordinate (a) at ($(0,0)+(45:4cm and 3cm)$);
  \coordinate (b) at ($(0,0)+(195:4cm and 3cm)$);
  \draw[Green, thick] (f1) -- (a) -- (f2);
\end{tikzpicture}
:::

:::{.fragment}
\begin{tikzpicture}%%width=100%
[every node/.style={font=\sf, color=Black}]
  \coordinate (a) at ($(0,0)+(45:4cm and 3cm)$);
  \coordinate (b) at ($(0,0)+(195:4cm and 3cm)$);
  \draw[very thick] (0,0) ellipse (4cm and 3cm);
  \node[above right] at (a) {$\varepsilon=0.661$};
  \draw[thick, Yellow] (0,0) ellipse (4cm and 2cm);
  \node[above,Yellow] at (0,2) {$\varepsilon=0.866$};
  \draw[thick, Blue] (0,0) ellipse (4cm and 1cm);
  \node[above,Blue] at (0,1) {$\varepsilon=0.968$};
\end{tikzpicture}
:::

::::

:::::
::::::

## Some Ellipse Math
You can always get the foci locations and eccentricity from the semi-major and minor axes.

::::::cols
::::col
- Foci:
  $$ f = \sqrt{a^2 - b^2} $$
- Eccentricity:
  $$ \varepsilon = \sqrt{1 - \frac{b^2}{a^2}} $$
- Area:
  $$ A = \pi a b $$
::::

::::col
\begin{tikzpicture}%%width=100%
[every node/.style={font=\sf, color=Black}]
\draw[very thick] (0,0) ellipse (3cm and 2cm);
\draw[latex-latex, Red, thick] (0,0) -- node[midway,below] {$a$} (3,0);
\draw[latex-latex, Purple, thick] (0,0) -- node[midway,left] {$b$} (0,2);
\fill[Cyan] (2.23,0) circle (3pt);
\draw[latex-latex, Cyan] (0,0.2) -- node[midway,above] {$f$} (2.23,0.2);
\fill[Cyan] (-2.23,0) circle (3pt);
\draw[latex-latex, Cyan] (0,0.2) -- node[midway,above] {$f$} (-2.23,0.2);
\end{tikzpicture}

::::
::::::


## Elliptical Orbits
- Most planets in the Solar System have orbits that are not very eccentric

:::{.fragment .current-visible .only-fragment}
\begin{tikzpicture}%%width=50%
[every node/.style={font=\sf, color=Black}]
\draw[thick,Cyan] (0,0) circle (3cm);
\node[left,align=center,Cyan,xshift=-5mm] at ($(0,0)+(180:3)$) {$\varepsilon=0$\\\small(circle)};
\end{tikzpicture}
:::

:::{.fragment .current-visible .only-fragment}
\begin{tikzpicture}%%width=55%
[every node/.style={font=\sf, color=Black}]
\draw[thick,Red] (0,0) ellipse (3cm and 2.91cm);
\node[right,align=center,Red] at ($(0,0)+(0:3)$) {$\varepsilon=0.24$\\\small(most eccentric\\ planet)};
\end{tikzpicture}
:::

:::{.fragment .current-visible .only-fragment}
\begin{tikzpicture}%%width=60%
[every node/.style={font=\sf, color=Black}]
\draw[thick,Cyan] (0,0) circle (3cm);
\node[left,align=center,Cyan,xshift=-5mm] at ($(0,0)+(180:3)$) {$\varepsilon=0$\\\small(circle)};
\draw[thick,Red] (0,0) ellipse (3cm and 2.91cm);
\node[right,align=center,Red] at ($(0,0)+(0:3)$) {$\varepsilon=0.24$\\\small(most eccentric\\ planet)};
\end{tikzpicture}
:::

:::{.fragment}
\begin{tikzpicture}%%width=60%
[every node/.style={font=\sf, color=Black}]
\draw[thick,Cyan] (0,0) circle (3cm);
\node[left,align=center,Cyan,xshift=-5mm] at ($(0,0)+(180:3)$) {$\varepsilon=0$\\\small(circle)};
\draw[thick,Red] (0,0) ellipse (3cm and 2.91cm);
\node[right,align=center,Red] at ($(0,0)+(0:3)$) {$\varepsilon=0.24$\\\small(most eccentric\\ planet)};
\draw[fill,Cyan] (0,0) circle (3pt);
\draw[fill,Red] (0.719,0) circle (3pt);
\end{tikzpicture}
:::

. . .

- The Sun is clearly not at the center however


<!--
## Kepler's 2nd Law
- A line joining a planet to the Sun sweeps out equal areas in equal times
	- Results in planets moving faster when close to the Sun and slower when far away

\begin{tikzpicture}%%width=60%
[every node/.style={font=\sf, color=Black}]
  \coordinate (sun) at (1.65,0);
  \draw[very thick] (0,0) ellipse (3cm and 2.5cm);
  \fill[opacity=0.5,Red] (sun) -- (-25:3cm and 2.5cm) arc (-25:25:3cm and 2.5cm);
  \fill[opacity=0.5,Red] (sun) -- (172:3cm and 2.5cm) arc (172:188:3cm and 2.5cm);
  \fill[inner color=yellow, outer color=orange] (sun) circle (2mm);
  \node[right] at (3,0) {Fastest};
  \node[left] at (-3,0) {Slowest};
\end{tikzpicture}

## Kepler's 3rd Law
:::incremental
- The square of a planet's orbital period is proportional to the cube of its semi-major axis.
	- A planet's _period_ is the time it takes to complete an entire orbit
	- An orbit's semi-major axis we defined earlier to be half the major axis
- Mathematically, this look like:
  $$ \frac{p^2}{a^3} = \text{constant} $$
  where $p$ is the period and $a$ is the semi-major axis
- For objects orbiting the Sun, if you choose units of years and AU, then
  $$ p^2 \approx a^3 $$
  - 1 _AU_ = 1 astronomic unit = average distance from Earth to Sun
:::


## Projections
- When observing, you rarely get a perfectly top-down view of an ellipse
- Often, you are looking at some view projected from the side
	- This does not change the measurement of the period!
	- It may make measuring the semi-major axis more difficult though
- Since many planets (and moons) have _mostly_ circular orbits, you can estimate the major axis by just taking the largest end-to-end distance that you observe


## Why Gravity?
:::::cols
::::col
- Kepler told us that there **must** be these relationships, but he couldn't say _why_
- Newton found that, if the direction that something is moving changes, then it **must** have experienced some force
	- Newton's big connection (one of them) was determining that the necessary force turns out to be from gravity (at least for most astronomic objects)
::::
::::col
\begin{tikzpicture}%%width=100%
[scale=0.7]
  \fill[inner color=yellow, outer color=orange] (0,0) circle (1cm);
  \draw[dashed, -latex, thick] (3,0) arc (0:340:3cm);
	\foreach \a in {0,45,...,350}{
	  \draw[orange, -latex] (\a:2.7) -- (\a:1.3);
	}
  \node at (3,0) {\includegraphics[width=1cm]{world.png}};
\end{tikzpicture}

::::
:::::

## Gravity
- Gravity is the universal attractor
	- Anything with mass attracts anything else with mass
	- Strength of force increases with the amount of mass involved
	- Strength of force decreases rapidly with distance between the masses

\begin{tikzpicture}%%width=70%
  \coordinate (ball1) at (0,0);
  \coordinate (ball2) at (8,0);
  \coordinate (d) at ($(ball1)-(0,1.5)$);
  \draw[Green] (ball1) --+(0,-2);
  \draw[Cyan] (ball2) --+(0,-2);
  \fill[ball color=Green!50] (ball1) circle (1cm);
  \fill[ball color=Cyan!50] (ball2) circle (1cm);
  \node at ($(ball1)+(0,1.2)$) {$M_1$};
  \node at ($(ball2)+(0,1.2)$) {$M_2$};
  \draw[<->, thick] (d) -- node[midway,below] {d} (d-|ball2);
  \draw[-latex, ultra thick,Red] (ball1) -- +(4,0);
  \draw[-latex, ultra thick,Red] (ball2) -- +(-4,0);
  \node[Red] at (4,1) {\scalebox{1.5}{$F_g = G \frac{M_1 M_2}{d^2}$}};
\end{tikzpicture}


## Newton meets Kepler's 3rd
- Kepler already had worked out
  $$ \frac{a^3}{p^2} = \text{same value for all planets orbiting Sun} $$
- Newton worked out, starting with the force, that two objects held in orbit by gravity would obey:
  $$ \frac{a^3}{p^2} = \frac{G(M_1 + M_2)}{4\pi^2} $$
  where:
  - $M_1, M_2$ are the masses of the objects in _kilograms_
  - $a$ is the average separation between the objects in _meters_
  - $p$ is the orbital period in _seconds_
  - $G$ is the gravitational constant ($6.67\times10^{-11}$)


## Some nicer units
- Put in more convenient units, Newton's formulation of Kepler's 3rd boils down to:
  $$ \frac{a^3}{p^2} = (M_1 + M_2)_\odot $$
  where
  - $M_1, M_2$ are the masses of objects in _solar masses_
  - $a$ is the average separation of the objects in _AU_
  - $p$ is the orbital period in _years_
- For the Sun and most planets, $M_1 + M_2 \approx 1 M_\odot$
- If you can measure $a$ and $p$, then you can work out the mass of the objects!
-->

## More Vocab
- _Apogee/Aphelion_: The point at which a body is at its furthest from the Earth/Sun
- _Perigee/Perihelion_: The point at which a body is at its closest to the Earth/Sun

\begin{tikzpicture}%%width=60%
[every node/.style={font=\sf, color=Black}]
  \coordinate (sun) at (1.65,0);
  \draw[very thick] (0,0) ellipse (3cm and 2.5cm);
  \fill[inner color=yellow, outer color=orange] (sun) circle (2mm);
  \node[circle, fill, inner sep=2pt, label={right: Perihelion}](p) at (3,0) {};
  \node[circle, fill, inner sep=2pt, label={left: Aphelion}](a) at (-3,0) {};
  \draw[stealth-stealth, thick, Red] (p) -- (sun);
  \draw[stealth-stealth, thick, Red] (a) -- (sun);
\end{tikzpicture}


# Ellipse Fitting
## Pesky Ellipses
- Because the elliptical orbits that space objects follow are tied directly to some physical properties, it can be very useful to determine elliptical parameters
- The classic equation of an ellipse:
  $$ \frac{(x-x_c)^2}{a^2} + \frac{(y-y_c)^2}{b^2} = 1 $$
  allows for no rotations
- The general polynomial form:
  $$ Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0 $$
  doesn't give obvious expressions for the desired parameters


## Fitting Tactics
- There are a variety of approaches, most perform best when the sampled points come from around the entire ellipse
	- Different approaches can error in different ways if all points are from one side of the ellipse
- All the pre-built options required some odd libraries and worked in very different ways between Python and R
- So I found starting code for both and then customized to create fitting options that work very similarly in both
	- In Python: [fitellipse.py](../libs/fitellipse.py)
	- In R: [fitellipse.r](../libs/fitellipse.r)


## Library Explanation
- Both libraries introduce 3 main functions. To get access to them, ensure that the `fitellipse.*` file is in the same directory as your notebook or source file, and then run or import the library
	- In Python: 
	  
	  ```python
	  run fitellipse.py
	  ```
	  or 
	  ```python
	  from fitellipse import *
	  ```
	- In R: 
	
	  ```R
	  source('fitellipse.r')
	  ```

## Library Functions
- The 3 main functions are:
	- `fit_ellipse(x, y)`: Computes the best fit ellipse and returns both the coefficients and the common fit parameters as an associative array. Takes points to fit as lists of x and y coordinates.
	- `get_ellipse(fit_params)`: Computes a series of x and y points representing the given ellipse, where `fit_params` is the output of `fit_ellipse`. Useful for visualization purposes.
	- `create_test_ellipse(Rx,Ry,Cx,Cy,Rotation,NoiseLevel)`: creates a series of x and y coordinates representing an ellipse with the given parameters


## Caveats
- Both fitting functions seem to do excellent at determining the center of the ellipse and major and minor axes
- There is some variability with the resulting angle, where it sometimes comes out $90^\circ$ off ($\pi/2$ technically since the angle is in radians)
	- I believe this stems from discrepancies with whether it ends up finding the angle from the major or minor axis, but I've spent hours trying to make it more robust and thus far failed
	- Fortunately, the angle of the ellipse has absolutely no bearing on any physical laws, and is only of interest for visualization purposes
		- If the angle of the ellipse is wrong when plotting, add or subtract $\pi/2$ from it


## Example with Comet NeoWISE
- Comet NeoWISE made an appearance several years back
	- That particular comet had a huge eccentricity and thus couldn't easily be visualized
- A different NeoWISE comet though has its positions over the next 100 years stored [here](../demos/neowise_positions.csv)
- How could I estimate the eccentricity of this comet?


## Activity!
- I have given you [here](../demos/pluto_positions.csv) the coordinates of Pluto over a 300 year span
- Your task: determining the closest approach Pluto makes with the Sun. That is to say, what is its perihelion distance?
	- Doing so will require working out the major axis and determining where the Sun would be
	- A picture might help you keep things straight as you go
- **Note**: If you compare this value to the known value, it isn't going to be perfect, owing to us not fully accounting for the tilt of Pluto's orbit when computing the ellipse