modelLightCurve.py

Tutorial Contents

• Generating a transit light curve for an exoplanet, HAT-P-7 b
• Extra: Generating a transit light curve for Earth

Generating a transit light curve with OSCAAR

Get the source for modelLightCurve.py here.

OSCAAR-extras/examples/modelLightCurve.py demonstrates how to generate a model light curve using the analytical formulation of Mandel & Agol (2002) with OSCAAR.

First, we will import the required packages to run this example.

import numpy as np
from matplotlib import pyplot as plt
import oscaar

Then we need to define the planetary system parameters according to the Mandel & Agol formalism. We enter the ratio of the radius of the planet to the radius of the star (often denoted R_p/R_s), the ratio of the semi-major axis to the radius of the star (a/R_s), the period, inclination (typically i), limb darkening coefficients (gamma_1, gamma_2), eccentricity, longitude of pericenter and the time of midtransit (t_0, or the "epoch"). The light curve that we'll replicate is that of transiting hot Jupiter HAT-P-7b.

## Define system parameters for planet HAT-P-7 b. 
## Citation: Morris et al. 2013 (http://arxiv.org/abs/1301.4503)
RpOverRs = 0.07759				## R_p/R_s = ratio of planetary radius to stellar radius
aOverRs = 4.0					## a/R_s = ratio of semimajor axis to stellar radius
period = 2.204737				## Period [days]
inclination = 83.111			        ## Inclination [degrees]
gamma1 = 0.3525					## Linear limb-darkening coefficient
gamma2 = 0.168					## Quadratic limb-darkening coefficient
eccentricity = 0.0				## Eccentricity
longPericenter = 0.0			        ## Longitude of pericenter
epoch = 2454954.357463			        ## Mid transit time [Julian date]

Now we'll input the times that we're simulating in Julian Dates for the light curve in units of days. We'll pretend to observe for 6 hours centered on mid-transit

## Generate times at which to calculate the flux
durationOfObservation = 6./24	## Elapsed time [days]
times = np.arange(epoch-durationOfObservation/2,epoch+durationOfObservation/2,1e-5)

We collect the system parameters in a list in the following order, and feed them as an argument into the oscaar.occultquad() function, along with the times. We'll save the model light curve to the array flux:

modelParams = [RpOverRs,aOverRs,period,inclination,gamma1,gamma2,eccentricity,longPericenter,epoch]

## Generate light curve
flux = oscaar.occultquad(times,modelParams)

oscaar.occultquad() is a Python wrapper around a C script that computes the equations for the quadratic limb-darkening case in Section 4 of Mandel & Agol (2002). This function can be a little cumbersome and slow when written in Python, so we took the liberty of wrapping a C version for you.

Let's plot the output!

## Plot light curve
plt.plot(times,flux,linewidth=2)
plt.title('Transit light curve for HAT-P-7 b')
plt.ylabel('Relative Flux')
plt.xlabel('Time (JD)')
plt.ylim([0.9925,1.0008])
plt.show()

What we see:

===

Generating A Transit Light Curve For Earth

Get the source code for this example,modelLightCurveEarth.py, here

If you'd like to see what the light curves would look like for other exoplanets, you can look up system parameters for most published exoplanets on exoplanets.org, and change the values in the modelParams accordingly. For example, what would Earth's transit light curve look like to an observer on a distant exoplanet? Surely the Earth is one planet for which we know the system parameters well, let's plot it's light curve! Let's tweak the parameters:

## Define system parameters for planet Earth, pulled from Wolfram Alpha
## Limb-darkening coeffs loosely based on Sanchez-Bajo et al. (2002) http://iopscience.iop.org/0143-0807/23/3/311
RpOverRs = 6367.5/695500.0      ## R_p/R_s = ratio of planetary radius to stellar radius
aOverRs = 1.49597887e8/695500.0 ## a/R_s = ratio of semimajor axis to stellar radius
period = 365.25                 ## Period [days]
inclination = 90.0              ## Inclination [degrees]
gamma1 = 0.7			## Linear limb-darkening coefficient
gamma2 = 0.0			## Quadratic limb-darkening coefficient
eccentricity = 0.016710220	## Eccentricity
longPericenter = 0.0		## Longitude of pericenter
epoch = 2454954.357463		## Mid transit time [Julian date]

In the previous example we used an observation of 6 hours, which is typical of most observations of transiting hot Jupiters, which orbit very closely to their stars. Since Earth orbits the Sun at a much greater distance than those hot Jupiters, and it has a correspondingly longer period, it takes Earth much longer to transit the Sun. You could see this as a direct effect of Johannes Kepler's "third law". So let's change the observation duration.

durationOfObservation = 16./24	## Elapsed time [days]

Then, as you might expect since RpOverRs is going to be so small, we'll want to zoom in close on the flux axis with the plt.ylim() plotting command:

plt.title('Transit light curve for Earth')
plt.ylabel('Relative Flux')
plt.xlabel('Time (JD)')
plt.ylim([0.999,1.001])
plt.show()

And we'll see the Earth's transit light curve:

Notice how much shallower this transit is. This is why NASA's Kepler mission has a very sensitive photometer in order to have any hopes of finding Earth-like planets.