exoplanet-dev / exoplanet

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

Exposure integration fails for large exposure times #39

Closed rodluger closed 5 years ago

rodluger commented 5 years ago

When the exposure time is comparable to the duration of the transit (which can occur for ultra-short period planets), the integration algorithm in exoplanet produces spurious wiggles. The code below compares a light curve integrated with exoplanet (orange line in plot) and one using a simple moving average (green line in plot). This occurs (with varying degrees) for all three integration order values and is not related to the boundary effects: even though I'm plotting the light curve only in the vicinity of transit, I'm computing it on a much wider time grid. Note that setting use_in_transit = False does not fix this issue.

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

def moving_average(a, n):
    """
    Compute a moving average over a window
    of `n` points. Based on
    https://stackoverflow.com/a/14314054
    """
    if n % 2 != 0:
        n += 1
    ret = np.cumsum(a, dtype=float)
    ret[n:] = ret[n:] - ret[:-n]
    result = ret[n - 1:] / n
    return np.concatenate([np.ones(n // 2) * result[0],
                           result,
                           np.ones(n // 2 - 1) * result[-1]])

# Generate a starry light curve
orbit = xo.orbits.KeplerianOrbit(period=1)
t = np.linspace(-0.5, 0.5, 1000)
u = [0.3, 0.2]
starry_op = xo.StarryLightCurve(u)

# Compute the light curve
lc = starry_op.get_light_curve(orbit=orbit, r=0.1, t=t).eval()

# Compute the light curve with finite exposure time
lc_int = starry_op.get_light_curve(orbit=orbit, r=0.1, t=t, texp=0.1, order=2).eval()

# Same thing, but numerically using a moving average
lc_int_num = moving_average(lc, int(0.1 / (t[1] - t[0])))

# Plot them
plt.plot(t, lc, lw=2, label="Flux")
plt.plot(t, lc_int, lw=2, label="Integrated with exoplanet")
plt.plot(t, lc_int_num, lw=2, label="Integrated with moving average")
plt.ylabel("relative flux")
plt.xlabel("time [days]")
plt.xlim(-0.1, 0.1)
plt.margins(None, 0.3)
plt.legend()
plt.show()

image

rodluger commented 5 years ago

Actually, increasing oversample improves the behavior, so perhaps this isn't really a bug. I find that I need to go to oversample > 25 in this particular case before the wiggles disappear. I wonder if there's a more efficient integration method for these cases, but perhaps I should just suck it up :)

dfm commented 5 years ago

I bet #11 will help this, but it's a bit of a pain to get it working properly within theano. I've been working on it and this will be a great use case to demonstrate how useful it is!

dfm commented 5 years ago

I've gone back and forth on this for a while and I've found that adding in the contact points doesn't help too much. We can solve this by using an adaptive algorithm and some sort of polynomial approximation to the orbit, but I find that that breaks down pretty badly for long exposure times like this. I'm going to punt on this for now and say that you're just going to have to use a "large" number for oversample when you want to have a long exposure time compared to the ingress and egress.