exoplanet-dev / exoplanet

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

Discontinuous behavior at large b when specifying both b and duration #130

Closed elnjensen closed 3 years ago

elnjensen commented 3 years ago

When using xo.orbits.KeplerianOrbit and specifying both duration and impact parameter b, and then generating a limb-darkened light curve with xo.LimbDarkLightCurve, some combinations of parameters give eclipses as expected, but at a slightly different impact parameter will switch over to return a flat light curve without any warnings or exceptions.

To Reproduce This sample code illustrates the issue:

import exoplanet as xo
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import warnings
warnings.filterwarnings("ignore", category=FutureWarning)

rp_over_rstar = 1
u = (0.5, 0.2)
t = np.linspace(0.25, 0.75, 100)
for b in (0.93, 0.94, 0.95, 0.96, 0.97, 0.98):
    orbit1 = xo.orbits.KeplerianOrbit(period=1, 
                                      t0=0.5, 
                                      b=b, 
                                      duration=0.05,
                                     )
    light_curve = xo.LimbDarkLightCurve(u).get_light_curve(
                            orbit=orbit1, 
                            r=rp_over_rstar, 
                            t=t
        ).eval()
    plt.plot(t, light_curve, label=f"$b={b}$")

plt.xlabel("Time (JD)")
plt.ylabel("Relative flux")
plt.legend()
plt.show()

That generates this plot:

exoplanet_b_duration

As shown, the curves with b = 0.96 and lower show eclipses, but at b = 0.97 and higher the behavior switches to a flat light curve.

Expected behavior Most likely this comes from over-determining the light curve so that it is no longer possible to meet both the b and duration constraints. If that is the case, some kind of warning or exception should be generated.

Further comments As I look at the plot, it makes me wonder what "duration" signifies in this case, even for the light curves that aren't flat. The eclipses have very different start to end durations. There is a crossing point of all of the eclipses, and the width there indeed matches the specified duration, but it's not clear to me how that is defined. Naively it seems that giving both b and duration should always overdetermine the light curve; the implicit M_star combined with the period will set the semimajor axis, and the specified R_planet / R_star, combined with the implicit R_star, will set R_planet. Thus, it would seem that a given impact parameter would uniquely determine the transit duration.

Versions: Python version : 3.8.3 IPython version : 7.16.1

exoplanet : 0.4.3 matplotlib: 3.2.2 theano : 1.0.4 numpy : 1.18.5

Mac OS X 10.15.7

Additional questions

When not specifying the duration while fitting a transit, is it possible to deterministically store the transit duration at each step in the chain?

Is duration defined as first to fourth contact, or in some other way?

Thanks for your attention to this!

ericagol commented 3 years ago

Maybe an issue with starry? @rodluger ?

dfm commented 3 years ago

It's defined using Eq. 14 from Winn (2010). I don't think it'll have anything to do with starry! I expect that there will be inconsistencies for large radius ratios but I'll take a closer look tomorrow.

elnjensen commented 3 years ago

Thanks for the quick attention to this! FYI, it can fail for small radius ratios too if you change the duration: exoplanet_b_duration2

(I realize that that's a pretty un-physical set of parameters when combined with Mstar = 1 and Rstar = 1, but part of my question is whether impossible orbits should be caught or flagged somehow.)

dfm commented 3 years ago

here's how it's defined btw:

img

elnjensen commented 3 years ago

Right, I can see that those points define a time interval that's the same across the light curves, but Equation 14 in Winn 2010 is for the total transit duration, first to fourth contact, while those points are all well after ingress / well before egress. So I'm not sure how to square those two things.

dfm commented 3 years ago

This is how it's calculated: https://github.com/exoplanet-dev/exoplanet/blob/d28c2a2949edd87dfaa355f0664260e591891b6a/src/exoplanet/orbits/keplerian.py#L829

There might be a bug, but I'm offline tonight so I can't check!

dfm commented 3 years ago

It's not a bug - it's because you need to pass ror to the orbit as well:

    orbit1 = xo.orbits.KeplerianOrbit(period=1, 
                                      t0=0.5, 
                                      b=b, 
                                      duration=0.05,
                                      ror=rp_over_rstar
                                     )
dfm commented 3 years ago

It would be good to add a warning to the code if duration is used when ror is not provided - good catch!

elnjensen commented 3 years ago

Ah, great, thanks. I can confirm that the example code gives sensible-looking lightcurves with that change.

elnjensen commented 3 years ago

I'd be happy to put together a simple PR for this if that's helpful. If so, should it be a warning or an exception? Seems like probably ValueError, since the results are incorrect otherwise and so it probably shouldn't go through.

And is master the right branch, or do you have a development branch that's better for a PR?

It may be simplest for you to just implement yourself, but I'd be happy to help if it's useful.

dfm commented 3 years ago

Thanks for the offer, Eric! Sure - why don't you open a PR to the main branch to raise an exception when the radius ratio is not provided? I'm also working on a major refactor of how all this parameterization stuff is handled for the next major release, but it still might be some time before that's ready so I think that it would be worth having this in the meantime!

dfm commented 3 years ago

Fixed in #161