adding rv tutorial

docs/_static/notebooks/citation.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [],
    "source": [

docs/index.rst

@@ -20,6 +20,7 @@ Tutorials
 
    tutorials/intro-to-pymc3
    tutorials/pymc3-extras
+   tutorials/rv
    tutorials/gp
 
 

exoplanet/estimators.py

@@ -2,8 +2,8 @@
 
 from __future__ import division, print_function
 
-__all__ = ["estimate_minimum_mass", "lomb_scargle_estimator",
-           "autocorr_estimator"]
+__all__ = ["estimate_semi_amplitude", "estimate_minimum_mass",
+           "lomb_scargle_estimator", "autocorr_estimator"]
 
 import numpy as np
 
@@ -29,8 +29,8 @@ def _get_design_matrix(periods, t0s, x):
     ] + [np.ones((1, len(x)))], axis=0).T
 
 
-def estimate_minimum_mass(periods, x, y, yerr=None, t0s=None, m_star=1):
-    """Estimate the minimum mass(es) for planets in an RV series
+def estimate_semi_amplitude(periods, x, y, yerr=None, t0s=None):
+    """Estimate the RV semi-amplitudes for planets in an RV series
 
     Args:
         periods: The periods of the planets. Assumed to be in ``days`` if not
@@ -42,20 +42,16 @@ def estimate_minimum_mass(periods, x, y, yerr=None, t0s=None, m_star=1):
         yerr (Optional): The uncertainty on ``y``.
         t0s (Optional): The time of a reference transit for each planet, if
             known.
-        m_star (Optional): The mass of the star. Assumed to be in ``M_sun``
-            if not an AstroPy Quantity.
 
     Returns:
-        An estimate of the minimum mass of each planet as an AstroPy Quantity
-        with units of ``M_jupiter``.
+        An estimate of the semi-amplitude of each planet in units of ``m/s``.
 
     """
     if yerr is None:
         ivar = np.ones_like(y)
     else:
         ivar = 1.0 / yerr**2
 
-    m_star = u.Quantity(m_star, unit=u.M_sun)
     periods = u.Quantity(np.atleast_1d(periods), unit=u.day)
     if t0s is not None:
         t0s = u.Quantity(np.atleast_1d(t0s), unit=u.day).value
@@ -71,7 +67,32 @@ def estimate_minimum_mass(periods, x, y, yerr=None, t0s=None, m_star=1):
     else:
         w = w[:-1]
         K = np.sqrt(w[::2]**2 + w[1::2]**2)
+    return K
+
+
+def estimate_minimum_mass(periods, x, y, yerr=None, t0s=None, m_star=1):
+    """Estimate the minimum mass(es) for planets in an RV series
+
+    Args:
+        periods: The periods of the planets. Assumed to be in ``days`` if not
+            an AstroPy Quantity.
+        x: The observation times. Assumed to be in ``days`` if not an AstroPy
+            Quantity.
+        y: The radial velocities. Assumed to be in ``m/s`` if not an AstroPy
+            Quantity.
+        yerr (Optional): The uncertainty on ``y``.
+        t0s (Optional): The time of a reference transit for each planet, if
+            known.
+        m_star (Optional): The mass of the star. Assumed to be in ``M_sun``
+            if not an AstroPy Quantity.
 
+    Returns:
+        An estimate of the minimum mass of each planet as an AstroPy Quantity
+        with units of ``M_jupiter``.
+
+    """
+    m_star = u.Quantity(m_star, unit=u.M_sun)
+    K = estimate_semi_amplitude(periods, x, y, yerr=yerr, t0s=t0s)
     m_J = K / 28.4329 * m_star.value**(2./3)
     m_J *= (periods.to(u.year)).value**(1./3)
 

exoplanet/light_curve.py

@@ -39,7 +39,7 @@ def __init__(self, u, r_star=1.0, model=None):
         self.c = get_cl(u_ext)
         self.c_norm = self.c / (np.pi * (self.c[0] + 2 * self.c[1] / 3))
 
-    def get_light_curve(self, r, orbit, t, texp=None, oversample=7, order=2,
+    def get_light_curve(self, r, orbit, t, texp=None, oversample=7, order=0,
                         use_approx_in_transit=True, duration_factor=3):
         """Get the light curve for an orbit at a set of times
 

exoplanet/orbits/keplerian.py

@@ -282,6 +282,10 @@ def get_planet_velocity(self, t):
     def get_star_velocity(self, t):
         """Get the star's velocity vector
 
+        .. note:: For a system with multiple planets, this will return one
+            column per planet with the contributions from each planet. The
+            total velocity can be found by summing along the last axis.
+
         Args:
             t: The times where the velocity should be evaluated.
 
@@ -293,19 +297,36 @@ def get_star_velocity(self, t):
         return tuple(tt.squeeze(x)
                      for x in self._get_velocity(self.m_planet, t))
 
-    def get_radial_velocity(self, t, output_units=None):
+    def get_radial_velocity(self, t, K=None, output_units=None):
         """Get the radial velocity of the star
 
         Args:
             t: The times where the radial velocity should be evaluated.
+            K (Optional): The semi-amplitudes of the orbits. If provided, the
+                ``m_planet`` and ``incl`` parameters will be ignored and this
+                amplitude will be used instead.
             output_units (Optional): An AstroPy velocity unit. If not given,
-                the output will be evaluated in ``m/s``.
+                the output will be evaluated in ``m/s``. This is ignored if a
+                value is given for ``K``.
 
         Returns:
-            The radial velocity evaluated at ``t`` in units of
-            ``output_units``.
+            The reflex radial velocity evaluated at ``t`` in units of
+            ``output_units``. For multiple planets, this will have one row for
+            each planet.
 
         """
+
+        # Special case for K given: m_planet, incl, etc. is ignored
+        if K is not None:
+            f = self._get_true_anomaly(t)
+            if self.ecc is None:
+                return tt.squeeze(K * tt.cos(f))
+            # cos(w + f) + e * cos(w) from Lovis & Fischer
+            return tt.squeeze(
+                K * (self.cos_omega*tt.cos(f) - self.sin_omega*tt.sin(f) +
+                     self.ecc * self.cos_omega))
+
+        # Compute the velocity using the full orbit solution
         if output_units is None:
             output_units = u.m / u.s
         conv = (1 * u.R_sun / u.day).to(output_units).value

exoplanet/theano_ops/kepler/solver.cc

@@ -84,7 +84,7 @@ int APPLY_SPECIFIC(solver)(
     M = M_in[n];
     e = e_in[n];
 
-    if (e >= 1) {
+    if (e > 1) {
       PyErr_Format(PyExc_ValueError, "eccentricity must be 0 <= e < 1");
       return 1;
     }