Discontinuities in light curve at high f values

[ben-cassese]

When generating transit light curves for planets with extreme flattening values (f roughly >0.6), unexpected and isolated spikes occasionally appear in the output:

import jax
jax.config.update("jax_enable_x64", True)
import jax.numpy as jnp
from squishyplanet import OblateSystem
import matplotlib.pyplot as plt

state = {
    "t_peri" : 0.0,
    "times" : jnp.linspace(0.0, 5, 2000), # units of days
    "a" : 2.0, # units of R_star
    "period" : 10, # units of days
    "projected_effective_r" : 0.1, # units of Rp/Rs
    "ld_u_coeffs":jnp.array([0.3, 0.1]), # u coefficients for a polynomial (here,
    # quadratic u1 and u2) limb darkening law
    "tidally_locked":False,
    "parameterize_with_projected_ellipse":True
}
planet = OblateSystem(**state)

spherical_lc = planet.lightcurve()
oblate_lc = planet.lightcurve({"projected_f" : 0.6})

fig, ax = plt.subplots(nrows=2, sharex=True, figsize=(8, 6))
ax[0].plot(planet.state["times"], spherical_lc, lw=2, label="spherical")
ax[0].plot(planet.state["times"], oblate_lc, lw=1, linestyle="--", label="oblate")
ax[0].set(ylabel="relative flux")
ax[0].legend()
ax[1].plot(planet.state["times"], (spherical_lc - oblate_lc)*1e6)
ax[1].set(xlabel="time [days]", ylabel="difference [ppm]")
fig.subplots_adjust(hspace=0);

These values are likely nonphysical (see Berardo and de Wit 2022), but could still cause problems for users who allow exploration of a wide parameter space.

Many thanks to @CalebLammers for pointing this out! Will try to hunt down the cause, but contributions are welcome.

[ben-cassese]

When working on #17, found some similar spikes and so implemented a fix there. The issue was cartesian_intersection_to_parametric_angle in squishyplanet.engine.parametric_ellipse. Previously, we were solving for this angle numerically, and at high $f$ values it sometimes did not converge. When revisiting it, realized there was actually a pretty simple analytic solution, so implemented that instead.

So, should be all set now! I've modified test_projected_parameterization to test $f$ values up to 0.99, and to look for single isolated spikes. The same code above now produces:

And changing projected_f to 0.99 now produces: