exoplanet-dev / exoplanet

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

Error in radial velocity for `t~p`, `omega=0.5*pi` #145

Closed rodluger closed 3 years ago

rodluger commented 3 years ago

Describe the bug The radial velocity returned by exoplanet.orbits.KeplerianOrbit().get_radial_velocity(t) jumps discontinuously when t is very close to t0 + period.

A few things I've observed:

To Reproduce

import numpy as np
import matplotlib.pyplot as plt
import exoplanet as xo

period = 365.25

orbit = xo.orbits.KeplerianOrbit(
    period=period,
    t0=0.0,
    incl=0.5 * np.pi,
    ecc=0.0,
    omega=0.5* np.pi,
    m_planet=3e-06,
    m_star=1.0,
    r_star=1.0,
)

t = np.linspace(0, period, 1000)
rv = orbit.get_radial_velocity(t=t).eval()
plt.plot(t, rv)
plt.xlabel("t [days]")
plt.ylabel("rv [m/s]");

The above code produces the following plot, where the discontinuity is evident at 365.25 days. The RV at that time should be zero, but it's equal to -0.06383679.

image

The sensitivity of 1e-13 can be seen in this figure:

delta = np.logspace(-15, -12, 1000)
rv = orbit.get_radial_velocity(t=period + delta).eval()
plt.plot(np.log10(delta), rv)
plt.xlabel("log |t - p|")
plt.ylabel("rv [m/s]");

image

Expected behavior For omega = 0.5 * np.pi, the RV at t = t0 + p should be exactly zero.

Your setup (please complete the following information):

Additional context This issue is unlikely to be encountered in practice because it happens at a pathological point. But that point (omega=0.5 * np.pi) is the default in starry, so I'm getting this discontinuity in several of my tutorial notebooks. If you're swamped I can try to dig into the Kepler solver code and see if I can find anything.

dfm commented 3 years ago

Thanks - all this info is super useful to track down the issue. I bet it has something to do with numerical issues the mod implementation in the kepler solver. I'll get it sorted!

ericagol commented 3 years ago

What is the behavior if you do not set the eccentricity? The Kepler solver shouldn't be needed for e=0, and it looks like the code should work fine for ecc=None: https://github.com/exoplanet-dev/exoplanet/blob/55d2252c71191044613fabb9c8bd3062aca3bc1b/src/exoplanet/orbits/keplerian.py#L188-L191

rodluger commented 3 years ago

Yep, it works fine with ecc=None.