exoplanet-dev / exoplanet

Fast & scalable MCMC for all your exoplanet needs!
https://docs.exoplanet.codes
MIT License
206 stars 52 forks source link

get_true_anomaly when mean anomaly is pi #31

Closed danhey closed 5 years ago

danhey commented 5 years ago

I've been using exoplanet quite a bit the past few days and noticed some weird behaviour in exoplanet.orbits.get_true_anomaly. When the mean anomaly is exactly pi, get_true_anomaly returns 0, whereas my implementation (and some others I've tried) return pi. You can see this in action below.

from exoplanet.utils import eval_in_model
from exoplanet.orbits import get_true_anomaly
import pymc3 as pm
import theano.tensor as tt

period = 10.
eccen = 0.5
asini = 100.
varpi = 1.

# TESS time sampling
time = np.arange(0, 27, 1.0 / (24 * 30))

with pm.Model() as model:
    M = tt.zeros_like(tt.constant(time)) + 2.0 * np.pi * (tt.constant(time)) / period
    f = get_true_anomaly(M, eccen + tt.zeros_like(M))
    tau_tens = (- (1 - tt.square(eccen)) * tt.sin(f+varpi) / (1 + eccen*tt.cos(f))) * (asini / 86400.)

    tau_val = eval_in_model(tau_tens)
    true_anom_val = eval_in_model(f)

mean_anom = np.zeros_like(time) + 2.0 * np.pi * (time) / period
# kepler() is a basic Newton's method solver I implemented
ecc_anom = kepler(mean_anom, eccen + np.zeros_like(mean_anom))
true_anom_python = 2.0 * np.arctan2(np.sqrt(1.0+eccen)*np.tan(0.5*ecc_anom),np.sqrt(1.0-eccen) + np.zeros_like(time))
tau_python = -(asini/86400) * (1.0 - np.square(eccen)) * np.sin(true_anom_python + varpi) / (1.0 + eccen*np.cos(true_anom_python))

residual

This flips the sign on the cos term in the tau variable, and makes the plot look a bit funky:

tau_compare

dfm commented 5 years ago

Is this with the GitHub master? I fixed this bug and added a test a few days ago.

danhey commented 5 years ago

I didn't think of that, thanks Dan!