Examples

The following examples should give a reasonably impression of how the transit model class can be used.

Calculating a model

In this first example, we demonstrate how to calculate a model transit light-curve.

from PyAstronomy.modelSuite import palTrans
import matplotlib.pylab as plt
import numpy as np

# Create a PalLC instance
plc = palTrans.PalLC()

# Set parameter values
plc["p"] = 0.1    # Planet radius / Stellar radius
plc["per"] = 1.0  # Orbital period
plc["a"] = 2.0    # Large semi major axis [R_S]
plc["i"] = 90.0   # Orbital inclination [deg]
# Specify limb darkening
# (quadratic limb-darkening law)
plc["linLimb"] = 0.4
plc["quadLimb"] = 0.2
# Specify T0 (central time of transit)
plc["T0"] = -0.1
# Specify binary contribution
plc["b"] = 0.0

# Check the parameters
plc.parameterSummary()

# Create a time axis
time = np.arange(1000)/1000.0 - 0.5

# Calculate the light curve using above set
# model parameters
lightcurve = plc.evaluate(time)

# Plot the result
plt.plot(time, lightcurve, 'bp')
plt.show()

Fitting a transit

The leading half of this example is very similar to the first one; we only use a smaller number of point for the time axis. The second half demonstrates how the fitting interface is invoked (more examples using simpler models can are shown in the funcFit tutorial).

Warning

If you are using mpmath for evaluation of the elliptical integrals, the calculations will be quite slow.

from __future__ import print_function, division
from PyAstronomy.modelSuite import palTrans
import matplotlib.pylab as plt
import numpy as np

# Create a PalLC instance
plc = palTrans.PalLC()

# Set parameter values
plc["p"] = 0.1    # Planet radius / Stellar radius
plc["per"] = 1.0  # Orbital period
plc["a"] = 7.5    # Large semi major axis [R_S]
plc["i"] = 90.0   # Orbital inclination [deg]
# Specify limb darkening
# (quadratic limb-darkening law)
plc["linLimb"] = 0.4
plc["quadLimb"] = 0.2
# Specify T0 (central time of transit)
plc["T0"] = -0.1
# Specify binary contribution
plc["b"] = 0.0

# Check the parameters
print("Input parameters: ")
plc.parameterSummary()

# Create a time axis
time = np.arange(100)/100.0 * 0.2 - 0.2

# Calculate the light curve using above set
# model parameters
lc = plc.evaluate(time)

# Save the result and add some noise
flux = lc + np.random.normal(0.0, 0.002, time.size)

# Now lets try to recover what we put in
# Choose some "guess" parameters
plc["p"] = 0.1     # Planet radius / Stellar radius
plc["per"] = 1.0   # Orbital period
plc["a"] = 7.5     # Large semi major axis [R_S]
plc["i"] = 90.0    # Orbital inclination [deg]
# Specify limb darkening
# (quadratic limb-darkening law)
plc["linLimb"] = 0.4
plc["quadLimb"] = 0.2
# Specify T0 (central time of transit)
plc["T0"] = -0.08
# Specify binary contribution
plc["b"] = 0.0

# Assume we want to fit "p", "a", "i", and "T0"
plc.thaw(["T0", "i"])

# Before we start fitting, check how the elliptical integrals
# are evaluated (mpmath or Boost)
print("Which elliptical integrals are used?: ", plc.whichEllInts())

# Carry out the fit
plc.fit(time, flux, yerr=np.ones(time.size)*0.002)

print("Fit parameters: ")
plc.parameterSummary()

# Plot the result
plt.plot(time, flux, 'bp')
plt.plot(time, plc.model, 'r-')
plt.show()

Obtain a model taking finite integration time into account

This example shows how to use the PalLC_Rebin class to take finite integration times and the resulting light-curve distortion into account. This example is very similar to the first one.

from PyAstronomy.modelSuite import palTrans
import matplotlib.pylab as plt
import numpy as np

# Create a PalLC_Rebin instance
plc = palTrans.PalLC_Rebin()

# Set parameter values
plc["p"] = 0.1    # Planet radius / Stellar radius
plc["per"] = 1.0  # Orbital period
plc["a"] = 2.0    # Large semi major axis [R_S]
plc["i"] = 90.0   # Orbital inclination [deg]
# Specify limb darkening
# (quadratic limb-darkening law)
plc["linLimb"] = 0.8
plc["quadLimb"] = 0.2
# Specify T0 (central time of transit)
plc["T0"] = -0.1
# Specify binary contribution
plc["b"] = 0.0

# Check the parameters
plc.parameterSummary()

# Create a time axis
time = np.arange(50)/50.0 - 0.51

# Specify oversampling parameters.
# Here use 10 points per observed bin.
plc.setRebinArray_Ndt(time, 10, time[1]-time[0])

# Calculate the light curve using above set
# model parameters
lc = plc.evaluate(time)

# Plot the result (both the overbinned and final
# model light-curves)
plt.plot(plc.rebinTimes, plc.unbinnedModel, 'b.-')
plt.plot(time, plc.binnedModel, 'rd--')
plt.legend(["Overbinned LC", "Averaged LC"])
plt.show()