---
title: Exoplanets 558 - Kepler Solutions
author: Andy Tzanidakis
tags: Exoplanets
---

# Exoplanets: Keplerian Orbit Solutions

Table of Content:
[TOC]

## Problem 1

### Part (a)
Here I present a simple Python function that solves Kepler's equation using Newtons method[^1]:

def kepler_equation(e_dat, E_dat):
    """ Analytic solution to Kepler's equation. Returns the mean anomaly.
        
        Parameters:
        -----------
        e_dat (float or array): Eccentricity of system 
        E_Dat (float or array): Eccentric anomaly of system
        
        Returns:
        -------
        mean_anomaly (float or array): Mean anomaly of orbital system (i.e the fraction of an elliptical orbit's period that has elapsed since the orbiting                                        body passed periapsis)
    """
    return E_dat - e_dat*np.sin(E_dat)

def kepler_solver(e, M, thresh=1e-10, timeout=5):
    """
        Solve the kepler equation assuming Newton's method. Returns the Eccentric anomaly.
        
        Parameters:
        -----------
        e: eccentricity
        M: Mean anomaly 
        thresh: solution threshold
        timeout: timeout if solution does not exits (returns False;) time is in seconds.
    """
    
    E = M + 0.85*e*np.sign(np.sin(M))  
    
    t0 = time() # initial time
    solv = 0 
    while np.abs(kepler_equation(e, E)-M)>thresh:
        E -= (kepler_equation(e, E)-M)/(1-e*np.cos(E)) # difference in each step
        tt = time() # calculate time 
        if tt-t0>timeout:
            solv = 1
            break
    
    if not solv:  
        return (E)
    else:
        return False

[^1]: I add a `timeout` argument to my Kepler function because some of my walkers in the Monte Carlo analysis seemed to give indefinite solutions that make the while loop never satisfied.

### Part (b)
Now I can check to see that my Kepler solver function works on the boundaries of eccentricity and mean anomaly

$$
\tag{1}
0 < e < 1
$$

$$
\tag{2}
0< M < 2\pi
$$

To check if the solver works, I build a uniform grid of 100x100 pixels in ecentricity and mean anomaly space. At each cell I run the `Kepler_solver` function to find the solution to the eccentric anomaly. In [Figure 1](https://i.imgur.com/6HKFK9L.png) we show the the values of the eccentric anomaly for each cell and see that our `Kepler_Solver` works well within the specified boundaries.

|<img src="https://i.imgur.com/6HKFK9L.png" alt="Image" width="410" height="350" style="display: block; margin: 0 auto" />|
|:--:|
| Figure 1. --Numerical solutiont grid to Kepler's equation via Newtons methood tested for limits in ecentricity (0<e<1) and mean anomaly (0<M<2$\pi$)|

## Problem 2
To calculate the radial velocity of the star, I will be using the funcitonal form:

$$
\tag{3}
V_{*,Z} = K_* cos(\omega) cos(f) - K_* sin(\omega) sin(f) + K_* e cos(\omega) + \gamma
$$

where the the true anomaly (f) is defined as:
$$
\tag{4}
f = 2 tan^{-1}\bigg{[} \bigg{(} \frac{1+e}{1-e}\bigg{)}^{1/2} tan\bigg{(}\frac{E}{2}\bigg{)}\bigg{]}
$$

Below is my implementation in Python:

def true_anomaly(e, M):
    """Calculate the true anomally.
        
        Parameters:
        -----------
        e: cccentricity
        M: Mean anomaly

    """
    E_model = kepler_solver(e, M) # calculate kelpers eqn
    return 2*np.arctan(np.sqrt((1+e)/(1-e)) * np.tan(E_model/2))

def Vr_star(t, Kstar, omega, P, tp, e, gamma):
    """Calculate the radial velocity for each step of time"""
    
    # Assuming that the time is an array
    calc_model = []
    
    # Solve Kepler's equation for each time step
    for time in (t):
        mean_anomaly = ((2*np.pi)/(P)) * (time-tp)

        ft = true_anomaly(e, mean_anomaly)
        vr = Kstar * np.cos(omega) * np.cos(ft) - Kstar*np.sin(omega)*np.sin(ft) + Kstar*e*np.cos(omega) + gamma
        calc_model.append(vr)
    
    return np.array(calc_model)

Using the equations we defined from above, we can test our solutions to a hypoethetical orbit to esnure that the equations are working. Here I will test my radial velcoty measurements as a function of time by varying the eccentricity.

|<img src="https://i.imgur.com/eF0OkvF.png" alt="Image" width="560" height="350" style="display: block; margin: 0 auto" />|
|:--:|
| Figure 1. -- Phased folded radial velocity curve using Kepler's equation. We vary the eccentricity for each orbit in bins of $\delta e$~0.1. |

## Problem 3
To find the period of the mysterious planet, I will be using the Lomb-Scargle Peridoogram (LSP) implementation from the open-source Python package [`gatspy`]([google](https://github.com/astroML/gatspy)) that supports the floating-mean implementation of LSP.

After experimenting with various frequency grids, for this analysis I fixed my frequency grid from 5 days to 150 days. It also appears that the underlying signal required higher order Fourier terms in the periodogram. For the final periodogram seen Figure 3, I chose N=10 fourier terms.

|<img src="https://i.imgur.com/NWZjLhh.png" alt="Image" style="display: block; margin: 0 auto" />|
|:--:|
| Figure 3. -- Periodogram of mystery planet on a period grid of 10 to 140 days. The yellow highlighted lines are the period of the highest three peaks in the periodogram.

In my LSP analysis, I found several high peaks all of which were approximately harmonics of a period of 55.7 and 111.45 days. Highlighted in yellow in Figure 3, I show the top three highest peaks in the periodogram. After testing and phase-folding the radial velocity curve, I found that the period that minimized the dispersion is $P=111.45$ days. This period is consisent with known exoplanet [HD 80606 b](https://en.wikipedia.org/wiki/HD_80606_b).

## Problem 4
For this problem, I now aim to constraint the orbital parameters of HD 80606b with the given radial velocity measurements. To constraint the parameters, I will be using a Markov chain Monte Carlo (MCMC) analysis (via `emcee`). Before I fit the data, I explored several combination of parameters to find the best start on initial values.

For this problem, I will assume a standard likelihood function without any additional scattered noise. I simultaneously fit all six orbital parameters parameters to the data.

In Figure 4, I show the corner plot of the posterior distribution of all six parameters. Overall, I find a good agreement in the recovery of parameters compared to [Naef et al. 2001](https://arxiv.org/pdf/astro-ph/0106256.pdf) -- particularly the very high value of eccentricity!

|<img src="https://i.imgur.com/vIGLKci.png" alt="Image" style="display: block; margin: 0 auto" />|
|:--:|
| Figure 4. -- Corner plot of the posterior distribution for all six orbital parameters constrained on the simulated radial velocity measurements from our Keplerian system. The dashed lines in each histogram represent the 0.18 and 0.84 quantiles of each distribution respectively.

In Figure 5, I show the best fit from the posterior distribution on a baseline of 1000 days including the residuals. I found that the residuals are roughly Gaussian and hence conclude that my fitting routine did a decent job fitting the RV data.

|<img src="https://i.imgur.com/Bq5cvLJ.png" alt="Image" style="display: block; margin: 0 auto" />|
|:--:|
| Figure 5. -- RV data plotted with the best fit model from the posterior distribution. Bottom figure is the residual between the observed and modeled RV values.

## Problem 6

### Part (a)

We want to show the conservation of angular momentum simplifies to:

$$
f\dot{g} - \dot{f}g=1
$$

$$
\vec{x} = f \vec{x_0} + g \vec{v_0}
$$

$$
\vec{v} = \dot{f} \vec{x_0} + \dot{g} \vec{v_0}
$$

We will begin with the conservation of angular momentum in the plane of the sky:

$$
\begin{equation}
\vec{x} \times \vec{v} = \vec{x_0} \times \vec{v_0}
\end{equation}
$$

First, let us begin by expanding Gauss's equations to the conservation of angular momentum:
$$
\vec{x} \times \vec{v} = (f \vec{x_0} + g \vec{v_0}) \times (\dot{f} \vec{x_0}+\dot{g}\vec{v_0})
$$

We can now begin to expand the cross product:
$$
\vec{x} \times \vec{v} = (f \vec{x_0} \times \dot{f} \vec{x_0}) + (f \vec{x_0} \times \dot{g} \vec{v_0}) + (g \vec{v_0} \times \dot{f} \vec{x_0}) + (g \vec{v_0} \times \dot{g} \vec{v_0})
$$

Note that a general property of cross product holds true:
$$
f\vec{x_0} \times \dot{f} \vec{x_0} =0, \space \space g \vec{v_0}\times\dot{g}\vec{v_0}=0
$$

Hence by applying the cross product rules we can further expand our terms:
$$
\vec{x} \times \vec{v} = (f \vec{x_0} \times \dot{g} \vec{v_0}) + (g \vec{v_0} \times \dot{f} \vec{x_0}) = (f \vec{x_0} \times \dot{g} \vec{v_0}) - (\dot{f}\vec{v_0} \times g \vec{v_0}) = (f\dot{g} - \dot{f} g) (\vec{x_0} \times \vec{v_0})
$$

Since we know from the conservation of angular momentum that we can simplify the term:

$$
\vec{x} \times \vec{v} = (f\dot{g} - \dot{f} g) (\vec{x} \times \vec{v})
$$

$$
\boxed{(f\dot{g} - \dot{f} g) =1}
$$

### Part (b)
Now we want to show that part (a) holds. We will assume that $E_0$=$M_0$=0,
and $M=n(t-$t_0$)=E - esin(E)$ (via Kepler's equation)

First, I will begin by simplifying the Gauss's terms:
$$
f = \frac{a}{r_0} (cos(E)-1)+1 = \frac{a}{r_0} (cos(E)-1+\frac{r_0}{a})
$$

$$
g=(t-t_0) + \frac{1}{n} \bigg{(} sin(E) - E\bigg{)} = \frac{1}{n} \bigg{(} M + sin(E) - E\bigg{)}
$$

$$
\dot{f} = - \frac{a^2 n sin(E)}{r r_0}
$$

$$
\dot{g} = \frac{a}{r} (cos(E)-1)+1 = \frac{a}{r} (cos(E)-1+\frac{r}{a})
$$

I will use the above relationships to further simplify. First, I assume $\delta=cos(E)-1$. Then we can continue expanding:

$$
f \dot{g} - \dot{f} g
$$

$$
\bigg{[} (\frac{a^2}{r_0}(\delta - \frac{r_0}{a})) * (\frac{a}{r}(\delta - \frac{r}{a}))\bigg{]} - \bigg{[} (\frac{-a^2 n sin(E)}{r r_0}) * (\frac{1}{n}(M+sin(E)-E))\bigg{]}
$$

$$
\bigg{[} (\frac{a^2}{r_0}(\delta - \frac{r_0}{a})) * (\frac{a}{r}(\delta - \frac{r}{a}))\bigg{]} - \bigg{[} \frac{-a^2 sin(E)(M+sin(E)-E)}{r r_0} \bigg{]}
$$

$$
\frac{a^2}{r r_0} \bigg{[} (\delta - \frac{r_0}{a} ) (\delta - \frac{r}{a}) + (sin(E)(M + sin(E)-E))\bigg{]}
$$

Now let's substitute back:

$$
\frac{a^2}{r r_0} \bigg{[} (cos(E)-1+\frac{r_0}{a})(cos(E)-1+\frac{r}{a})+sin(E)(E-esin(E)+sin(E)-E) \bigg{]}
$$

$$
\frac{a^2}{r r_0} \bigg{[} (cos(E)-1+\frac{a(1-e)}{a})(cos(E)-1+\frac{a(1-ecos(E))}{a})+sin(E)(E-esin(E)+sin(E)-E) \bigg{]}
$$

$$
\frac{a^2}{r r_0} \bigg{[} cos^2(E) - ecos^2(E) - ecos(E) + e^2cos(E) + sin^2(E)(1-e) \bigg{]}
$$

$$
\frac{a^2}{r r_0} \bigg{[} cos^2(E) + sin^2(E) - e(cos^2(E) + sin^2(E)) - ecos(E) + e^2 cos(E) \bigg{]}
$$

$$
\frac{a^2}{r r_0} \bigg{[} 1-e-e cos(E) + e^2 cos(E)\bigg{]}
$$

$$
\frac{a^2}{r r_0} \bigg{[} (1-e)(1-ecos(E)) \bigg{]}
$$

Which by substitution shows us that
$$
\frac{a^2}{r r_0} \frac{r r_0}{a^2} = 1 = f\dot{g} - \dot{f} g
$$

Hence we have sucessfully demonstrated that angular momentum is conserved!

## Acknowledgements
For this homework I worked in collaboration with Tom Wagg, Bonni Choi, Sam Garza, David Wang, and Megan Gialluca.