Here we'll briefly discuss how to account for finite exposure time integration, which tends to smooth out occultation light curves. Unfortunately, there's no analytic way (that I know of) to tackle this, so the thing to do is to oversample the light curve on a fine time grid and numerically integrate over the exposure window.
%matplotlib inline
%run notebook_setup.py
import numpy as np
import starry
starry.config.lazy = False
starry.config.quiet = True
Let's instantiate a simple Primary object:
star = starry.Primary(starry.Map(ydeg=0, udeg=2, amp=1.0), m=1.0, r=1.0)
And let's give it some limb darkening:
star.map[1] = 0.40
star.map[2] = 0.26
Here's what that looks like:
star.map.show()
Let's now create a featureless hot Jupiter:
planet = starry.kepler.Secondary(
starry.Map(2, amp=0), m=0, r=0.1, porb=1.0, ecc=0.3, w=30, t0=0,
)
Now that we have a star and a planet, we can instantiate the planetary system. By default, exposure time integration is disabled.
sys = starry.System(star, planet)
We're ready to compute a transit light curve. Let's do this over 1000 cadences between $t=-0.1$ and $t=0.1$.
# HACK: run this to pre-compile the flux method
sys.flux(0.0);
%%time
time = np.linspace(-0.1, 0.1, 1000)
flux = sys.flux(time)
Cool -- starry computed 1,000 cadences in just a few ms. Let's check it out:
plt.plot(time, flux)
plt.xlabel("time [days]")
plt.ylabel("system flux");
Ok, now let's enable exposure time integration. We have to instantiate a new System object for this:
sys_exp = starry.System(star, planet, texp=0.01, oversample=9, order=0)
We passed in the three keywords controlling exposure time integration. The first is texp, the length of the exposure window in sys_exp.time_unit (usually days). The second is oversample, the factor by which the light curve is oversampled. The larger this number, the more accurate the model will be, at the cost of extra computation time. Finally, order controls the order of the numerical integration: 0 for a centered Riemann sum
(equivalent to the "resampling" procedure suggested by Kipping 2010), 1 for the trapezoid rule, or 2 for Simpson’s rule.
# HACK: run this to pre-compile the flux method
sys_exp.flux(0.0);
Let's compute the flux (and time the computation):
%%time
flux_exp = sys_exp.flux(time)
That was a factor of a few slower than the original evaluation, but it's not bad. Here's the comparison of the two light curves:
plt.plot(time, flux, label=r"$t_{exp} = 0$")
plt.plot(time, flux_exp, label=r"$t_{exp} = 0.01$")
plt.legend()
plt.xlabel("time [days]")
plt.ylabel("system flux");
Exposure time integration also affects phase curves. Let's give the planet a random map and compute its phase curve with and without exposure time integration. We'll make the planet's rotational period really short so we can clearly tell the difference between the two:
planet.map.amp = 0.1
planet.prot = 0.05
planet.map[1:, :] = 0.1 * np.random.randn(planet.map.Ny - 1)
planet.map.show()
Let's grab just the planet flux (by passing total=False to the flux method and keeping only the second vector):
%%time
flux = sys.flux(time, total=False)[1]
%%time
flux_exp = sys_exp.flux(time, total=False)[1]
Here are the two light curves; it's clear that the finite exposure time has smoothed out the phase curve.
plt.plot(time, flux, label=r"$t_{exp} = 0$")
plt.plot(time, flux_exp, label=r"$t_{exp} = 0.01$")
plt.legend()
plt.xlabel("time [days]")
plt.ylabel("planet flux");