TTVOrbit for StarryLightCurve about exoplanet

[dfm]

This functionality has not been carefully tested -- that's also why it's not documented :) -- so use at your own risk!

That being said, the bug here is that the transit_times or ttvs arguments must be lists of arrays (or list of tensors, or list of lists), one for each body. You could change it to transit_times=[[1, 2, 3]] and it should work.

[jiayindong]

Cool, thanks!

[gjgilbert]

@jiayindong @dfm Just a heads up that even using the list-of-arrays workaround, I've been running into issues with TTVOrbit. For some reason, not all the transits are being generated, or are perhaps offset by an unusually large TTV amplitude (see top figure below). The reference planet here is Kepler-18c (P=7.6 days), and should have more or less evenly spaced transits, with a TTV amplitude of ~10 min. As you can see, many transits are missing.

Of course, it's possible I've made a mistake in my implementation somewhere, although everything seems to work fine when I use KeplerianOrbit (bottom figure).

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

def calculate_model_flux(pshape, Rstar, tts, time_stamps, cadences):
    '''
    pshape: array of transit parameters -- [T0, P, rp, zeta, b^2, esinw, ecosw, u1, u2]
    Rstar: stellar radius [solar radii]
    tts: list of transit times
    time_stamps: list of time stamps (one per transit)
    cadences: cadence of each transit
    '''
    # check that vector lengths are consistent
    if len(time_stamps) != len(tts):
        raise ValueError('inconsistent input lengths')
    if len(cadences) != len(tts):
        raise ValueError('inconsistent input lengths')

    # convert pshape parameters to a form more easily used by starry
    T0, P, rp, a, inc, ecc, w, u1, u2 = pshape_to_pstarry(pshape, Rstar)
    b = np.sqrt(pshape[4])

    # setup exo.StarryLightCurve() object
    exoSLC = exo.StarryLightCurve([u1,u2])

    orbit = exo.orbits.TTVOrbit(transit_times=[tts], , a=a, r_star=Rstar)    

    model = exoSLC.get_light_curve(orbit=orbit, r=rp, t=np.hstack(time_stamps)).eval()

    return model

# 'planets' is a list of objects which hold relevant planet parameters, ttv info, etc    
p = planets[0]

tts = p.tts[p.quality]                      # a list of transit times
time_stamps = p.grab_stamps('time')         # a list of 'stamps' cut for 3 transit durations around each tt
flux_stamps = p.grab_stamps('flux')         # a list of correspnding flux values
cadences = p.stamp_cadence[p.quality]       # Kepler cadence of each stamp ('short' or 'long')

model = wings.calculate_model_flux(p.pshape, RSTAR, tts, time_stamps, cadences)

plt.figure(figsize=(18,6))
plt.plot(np.hstack(time_stamps), np.hstack(flux_stamps), 'k.')
plt.plot(np.hstack(time_stamps), model+1, 'orange')
plt.show()

[dfm]

Can you send a simple executable example that doesn't have the behavior that you expect and I'll take a look? Thanks!

On Thu, Jul 18, 2019 at 5:11 PM Jiayin Dong ***@***.***> wrote: Reopened #49 <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_dfm_exoplanet_issues_49&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=AvKCBU-zyph1qdf6NaatMQ&m=Ww5gFI6EIGxiLCRLoEp28hwySFMZ9G9HM-xHyvfGn-g&s=HbWs5F3btRM92qAdikA5M9mr53zeBeOz8bWwmhUrj1Q&e=> . — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_dfm_exoplanet_issues_49-3Femail-5Fsource-3Dnotifications-26email-5Ftoken-3DAACVQSSFXHDMQT6XVNFVRPTQADL6ZA5CNFSM4IEWPNLKYY3PNVWWK3TUL52HS4DFWZEXG43VMVCXMZLOORHG65DJMZUWGYLUNFXW5KTDN5WW2ZLOORPWSZGOSSUAFOA-23event-2D2494038712&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=AvKCBU-zyph1qdf6NaatMQ&m=Ww5gFI6EIGxiLCRLoEp28hwySFMZ9G9HM-xHyvfGn-g&s=bRtuSOZ42_zgcb6GDudHQXMx_pSfve8INlLTn6ab-uE&e=>, or mute the thread <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_AACVQSSIAI7QU4TGCES5USDQADL6ZANCNFSM4IEWPNLA&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=AvKCBU-zyph1qdf6NaatMQ&m=Ww5gFI6EIGxiLCRLoEp28hwySFMZ9G9HM-xHyvfGn-g&s=tfnDSnQhet5kJPx1D9mP1AX2AWjgkxNfqjLL34Np3lY&e=> .

-- Dan Foreman-Mackey Associate Research Scientist Flatiron Institute http://dfm.io

[gjgilbert]

I've uploaded my code here: https://github.com/gjgilbert/pterodactyl. Everything should be executable from the jupyter notebook. You'll need to specify a few paths in the first two cells but everything else should run without issue. The behavior in question is in the last few cells under the section "do a model fit" and the relevant function in pterodactyl_wings is "calculate_model_flux"

The lightcurve detrending takes a few minutes, so if you'd prefer me to put together a slimmed-down version of this code let me know and I'll be happy to do it.

[dfm]

Thanks! Yeah - it would be best if you could put together a simple standalone version that isolates the problem. I have a few ideas, but like I said this definitely hasn't been well tested so it'll take some poking!

[gjgilbert]

Okay, a notebook 'dfm_exoplanet_testing' is uploaded here: https://github.com/gjgilbert/pterodactyl, plus three .fits files containing relevant data for the three planets in the Kepler-18 system. All you need to do is specify which file you want to read in (keyword "FITSFILE" in the second cell). FYI, the naming convention I've used sorts the planets by orbital period, which is not necessarily their actual KOI or Kepler ID.

[dfm]

Hi Greg, sorry for the slow response over here! Can you make a little toy example that has the same behavior without all the surrounding infrastructure. Like one single standalone function with no dependencies or files. Then we'll be able to see more clearly where the issue is coming from. Thanks!

[gjgilbert]

No worries about the response time, Dan!

Here's a completely stripped-down notebook that reproduces the issue: https://github.com/gjgilbert/pterodactyl/blob/master/dfm_exoplanet_testing.ipynb

[dfm]

@gjgilbert: I got a chance to look at this and I worked out what was going on. I'd developed the model in the limit where the ttvs were less than the period, but here you have some transits that come more than a period late so that was messing up all the bookkeeping. I just pushed a change that looks a bit better, but there's still an issue when the transit comes >2 "periods" late. In the example you sent, one of your transits is 10 days after the previous one - are you sure that there shouldn't be another transit in there that you're missing? The way the code is set up, you need to provide a transit time for every transit (even the ones that aren't observed). I'll work on making this clearer, but let me know if it's reasonable to expect that some systems will have TTVs larger than the reference period....

[dfm]

Edit: change all occurrences of period to period/2 in the above comment - sorry!