orbital motion - Calculating the eccentricity of an exoplanet

I'm wondering how to calculate the eccentricity of an exoplanet by its radial velocity vs. phase graph. To clarify my question I will take an exoplanet called WASP-14b 2 as an example (http://exoplanets.org/detail/WASP-14_b).

A plot of the radial velocity of the star vs the phase is displayed in the upper left corner. I am wondering how I could possibly calculate the eccentricity of the exoplanet using this graph (or some other values given in the original measurements). I found a few ways to calculate the eccentricity:

$$ e = \left | \mathrm{e} \right| $$

This uses the eccentricity vector which is calculated using this formula:

$$ \mathrm{e} = \frac{v \times h}{\mu}-\frac{r}{\left | r \right |} $$

The problem here is that this formula needs the specific angular momentum vector and the position vector, which I do not know given only the measurements. However, there is another way to calculate the eccentricity:

$$ e = 1 - \frac{2}{(r_a/r_p) + 1} $$

where $r_a$ is the radius of the apoapsis and $r_p$ the radius of the periaosis. These values are not known using only the measurements, but I believe it should be possible to calculate them by taking the integral of the sine function (radial velocity vs. phase). This would give me the position of the star at any given moment. The problem is that I cannot find the exact points displayed in the graph anywhere, let alone a sine function that would fit them.

When I do get an integral of the function I still have to create one for the planet itself, since this describes the movement of the star. I am able to calculate the mass of the planet using the following formula:

$$ r^3 = \frac{GM_{star}}{4\pi^2}P_{star}^2 $$

which gives me the distance between the star an the planet. Next I can calculate the velocity of the planet using:

$$ V_{PL} = \sqrt{GM_{star}/r} $$

And after that I can calculate the mass of the planet using this formula:

$$ M_{PL} = \frac{M_{star}V_{star}}{V_{PL}} $$

But this is where another problem comes up as a Wikipedia article on Doppler Spectroscopy states: "Observations of a real star would produce a similar graph, although eccentricity in the orbit will distort the curve and complicate the calculations below."

Where do I find the corrected calculations and how can I possibly calculate the eccentricity of this planet using only these values ($M_{star}$ and the plot, of which I cannot find the exact points)?

Additional sources: http://adsabs.harvard.edu/abs/2009MNRAS.392.1532J

Answer

There are a number of options if you want an off-the-shelf solution to fitting RV curves. Perhaps the best free one is Systemic Console.

However, it is not too hard to do something basic yourself.

First define some terms:

$\nu(t)$ is the true anomaly - the angle between the pericentre and the position of the body around its orbit, measured from the centre of mass focus of the ellipse.

$E(t)$ is the eccentric anomaly and is defined through the equation $$\tan \frac{E(t)}{2} = \left(\frac{1+e}{1-e}\right)^{-1/2} \tan \frac{\nu(t)}{2}$$

The mean anomaly $M(t)$ is given by $$M(t) = \frac{2\pi}{p}(t - \tau),$$ where $p$ is the orbital period and $\tau$ is the time of pericentre passage.

"Kepler's equation" tells us that $$M(t) = E(t) - e \sin E(t)$$

Finally, the radial velocity is given by $$V_r(t) = K\left[\cos(\omega + \nu(t)) +e \cos \omega \right] + \gamma,$$ where $K$ is the semi-amplitude, $\gamma$ is the centre of mass radial velocity and $\omega$ is the usual angle defining the argument of the pericentre measured from the ascending node.

OK, so the problem is that the radial velocity does not depend explicitly on $t$, but rather on $\nu$. So what you do is the following:

1. 
   

   Choose values for $K, \gamma$, $\omega$, $\tau$, $p$ and $e$; these are your "free parameters that describe the orbit. The closer you can get your initial guess, the better.

   
   



2. 
   

   You use these parameters to predict what the radial velocities would be at the times of observation of your RV datapoints. You do this by calculating $\nu(t)$ using the equations above. Start with the second equation and calculate $M(t)$. Then you have to solve the third equation to get $E(t)$. This is transcendental, so you have to use a Newton-Raphson method or something similar to find the solution. Once you have $E(t)$ then you use the first equation to find $\nu(t)$. Then use 4th equation to calculate $V_r(t)$ at each of your datapoint times.

   
   


3. 
   

   Calculate a chi-squared (or similar figure of merit) from comparing the predicted and measured values of $V_r(t)$.

   
   


4. 
   

   Iterate the values of the free parameters and go back to step 2. Continue till your fit converges.