Analytical transit model (Mandel & Agol 2002)

The forTrans module provides easy access to the FORTRAN routines by Mandel & Agol 2002 1, which can be used to calculate transit light curves. While the FORTRAN code ensures fast computation of the light curve, PyAstronomy supplies funcFit’s fitting framework.

Note

The code contained within the files occultnl.f and occultquad.f is used with kind permission of the authors. It is described in 1; further explanations are given in Eastman et al. 2013 (PASP 125, 83).

1(1,2)

Mandel & Agol 2002, “Analytic Light Curves for Planetary Transit Searches”, ApJ, 580, 171.

The MandelAgolLC class

The MandelAgolLC model can be used to calculate transit light curves with quadratic or non-linear limb darkening based on a circular or general Keplerian orbit.

class PyAstronomy.modelSuite.XTran.forTrans.mandelAgol.MandelAgolLC(orbit='circular', ld='quad', collCheck=True, pyfo=False)

Analytical transit light-curves using the formulae provided by Mandel & Agol 2002.

Note

The computation of transit light curves is done using the external occultquad FORTRAN library.

This library can be installed, e.g., via

pip install PyAstronomy_ext

It can also be compiled manually using SciPy’s f2py wrapper (http://www.scipy.org/F2py). Simply go to the forTrans directory of the source distribution of PyAstronomy, then invoke

f2py -c occultquad.pyf occultquad.f

f2py -c occultnl.pyf occultnl.f

If no FORTRAN implementation is available, a python re-implementation is used. Performance may be impacted.

Model parameters

The set of parameters specifying this model depends on: the type of orbit chosen (circular or keplerian) and the type of limb darkening chosen (quadratic or non-linear).

More information on the Keplerian orbit can be found here: Calculate a Keplerian (two body) orbit

Orbital model parameters (circular orbit):
  • p - Radius ratio between planet and star.

  • a - Semi-major axis of planetary orbit [stellar radii].

  • i - Inclination of orbit in degrees (90 deg is edge on view).

  • T0 - Time offset of transit center.

  • per - Period of planetary orbit.

  • b - Describes the flux ratio between a stellar companion and the main star (default is 0).

Orbital model parameters (Keplerian orbit):
  • p - Radius ratio between planet and star.

  • a - Semi-major axis of planetary orbit [stellar radii].

  • i - Inclination of orbit in degrees (90 deg is edge on view).

  • per - Period of planetary orbit.

  • b - Describes the flux ratio between a stellar companion and the main star (default is 0).

  • tau - Time of periapsis passage.

  • Omega - Longitude of the ascending node [deg].

  • w - Argument of periapsis [deg]. Note that the longitude if periapsis is given by Omega+w.

  • e - Orbital eccentricity (0-1).

  • sed - Secondary eclipse depth

Limb darkening parameters (quadratic):
  • linLib - Linear limb-darkening coefficient.

  • quadLimb - Quadratic limb-darkening coefficient.

Limb darkening parameters (non-linear):
  • a1 - Non-Linear limb-darkening coefficient.

  • a2 - Non-Linear limb-darkening coefficient.

  • a3 - Non-Linear limb-darkening coefficient.

  • a4 - Non-Linear limb-darkening coefficient.

Limb-darkening laws

The quadratic limb-darkening law is given by:

\[\frac{I(\mu)}{I(1)}= 1 - linLimb \times (1-\mu) - quadLimb \times (1-\mu)^2\]

The non-linear limb-darkening law is given by:

\[\frac{I(\mu)}{I(1)}= 1 - \sum_{n=1}^{4}{a_n(1-\mu^{n/2})}\]
Modeling of secondary eclipse

If the parameter sed is not equal zero, a secondary eclipse (occultation) is added to the model. The secondary eclipse is modeled via a purely geometric eclipse of the stellar and planetary disk, which would then be behind the star (no limb darkening). The parameter sed specifies the depth of the secondary eclipse when the planetary disk is completely covered by the stellar disk, i.e., the depth may not be reached in a grazing configuration.

Warning

Time units have to be consistent.

Parameters
orbitstring, {“circular”, “keplerian”}, optional

Determines whether a circular or full keplerian planetary orbit is used in the calculations. The default is a circular orbit.

ldstring, {“quad”, “nl”}

Determines whether quadratic or non-linear limb darkening is used. The default is quadratic limb darkening.

collCheckboolean, optional

If True (default), the model will check whether there is a physical collision between the star and the planet on evaluating the model and raises an exception when there is one.

pyfoboolean, optional

Use python implementation of FORTRAN routines. Default is False; if FORTRAN routines are unavailable, this is the fallback option. Performance may be impacted.

Methods

MCMCautoParameters(ranges[, picky, ...])

Convenience function to generate parameters for MCMC fit.

addConditionalRestriction(*args)

Define a conditional restriction.

assignValue(specval)

Assign new values to variables.

assignValues(specval)

Assign new values to variables.

autoFitMCMC(x, y, ranges[, picky, stepsize, ...])

Convenience function to using auto-generated sampling parameters in MCMC.

availableParameters()

Provides a list of existing parameters.

backendStatus()

Print information on the code being used for evaluation

delRestriction(parName)

Delete restriction

description([parenthesis])

Returns a description of the model based on the names of the individual components.

errorConfInterval(par[, dstat, statTol, ...])

Calculate confidence interval for a parameter.

evaluate(time)

Calculate a light curve according to the analytical models given by Mandel and Agol.

fit(x, y[, yerr, X0, minAlgo, mAA, ...])

Carries out a fit.

fitEMCEE([x, y, yerr, nwalker, priors, ...])

MCMC sampling using emcee package.

fitMCMC(x, y, X0, Lims, Steps[, yerr, ...])

Carry out MCMC fit/error estimation.

freeParamNames()

Get the names of the free parameters.

freeParameters()

Get names and values of free parameters.

freeze(specifiers)

Consider variables free to float.

frozenParameters()

Get names and values of frozen parameters.

getRelationsOf(specifier)

Return relations of a variable.

getRestrictions()

Get all restrictions.

hasVariable(specifier)

Determine whether the variable exists.

numberOfFreeParams()

Get number of free parameters.

parameterSummary([toScreen, prefix, sorting])

Writes a summary of the parameters in text form.

parameters()

Obtain parameter names and values.

relate(dependentVar, independentVars[, func])

Define a relation.

removeConditionalRestriction(*args)

Remove an existing conditional constraint.

renameVariable(oldName, newName)

Change name of variable.

restoreState(resource)

Restores parameter values from file or dictionary.

saveState(*args, **kwargs)

Save the state of the fitting object.

setObjectiveFunction([miniFunc])

Define the objective function.

setPenaltyFactor(penalFac)

Change the penalty factor.

setRestriction(restricts)

Define restrictions.

setRootName(root[, rename])

Define the root name of the model.

showConditionalRestrictions(**kwargs)

Show conditional restrictions.

steppar(pars, ranges[, extractFctVal, quiet])

Allows to step a parameter through a specified range.

thaw(specifiers)

Consider variables fixed.

untie(parName[, forceFree])

Remove all relations of parameter parName, i.e., the parameter is not dependend on other parameters.

updateModel()

Recalculate the model using current settings.

Example: Calculate model light curve (circular orbit, quadratic limb darkening)

# Import some unrelated modules
import numpy as np
import matplotlib.pylab as plt
# ... and now the forTrans module
from PyAstronomy.modelSuite import forTrans as ft


# Create MandelAgolLC object with
# circular orbit and quadratic limb darkening
ma = ft.MandelAgolLC(orbit="circular", ld="quad")

# See the available parameters and their current values
ma.parameterSummary()

# Set parameters
ma["per"] = 0.2
ma["i"] = 90.0
ma["a"] = 6.5
ma["T0"] = 0.0
ma["p"] = 0.16
ma["linLimb"] = 0.47
ma["quadLimb"] = 0.24
ma["b"] = 0.

# Choose some time axis
time = np.linspace(0, 0.5, 1000)

# ... and calculate model
y = ma.evaluate(time)

# Let's see what happened ...
plt.plot(time, y, 'b.')
plt.show()

Example: Calculate model light curve (keplerian orbit, quadratic limb darkening)

# Import some unrelated modules
import numpy as np
import matplotlib.pylab as plt
# ... and now the forTrans module
from PyAstronomy.modelSuite import forTrans as ft


# Create MandelAgolLC object with
# keplerian orbit and quadratic limb darkening
ma = ft.MandelAgolLC(orbit="keplerian", ld="quad")

# See the available parameters and their current values
ma.parameterSummary()

# Set parameters
ma["per"] = 0.2
ma["i"] = 88.76
ma["a"] = 6.5
ma["p"] = 0.16
ma["linLimb"] = 0.47
ma["quadLimb"] = 0.24
ma["b"] = 0.
ma["e"] = 0.75
ma["w"] = 127.
ma["Omega"] = 3.9

# Choose some time axis
time = np.linspace(0, 0.5, 1000)

# ... and calculate model
y = ma.evaluate(time)

# Let's see what happened ...
plt.plot(time, y, 'b.')
plt.show()

Example: Comparing limb-darkening laws

# Import some modules
import numpy as np
import matplotlib.pylab as plt
# ... and now the forTrans module
from PyAstronomy.modelSuite import forTrans as ft

# First, let's compute a transit model using
# quadratic limb-darkening prescription.
ma = ft.MandelAgolLC(ld="quad")

# Set parameters. The LD coefficients are taken
# from Claret 2011 for a solar-metallicity star
# with Teff=6000 K and logg=4.5.
ma["per"] = 0.2
ma["i"] = 90.0
ma["a"] = 6.5
ma["T0"] = 0.5
ma["p"] = 0.16
ma["linLimb"] = 0.0479
ma["quadLimb"] = 0.2716
ma["b"] = 0.

# Choose some time axis
time = np.linspace(0, 0.2, 1000)

# ... and calculate model
yQLD = ma.evaluate(time)

# Now, let's compute a transit model with
# non-linear limb-darkening prescription
# for the same stellar parameters.
maNL = ft.MandelAgolLC(ld="nl")
maNL["per"] = 0.2
maNL["i"] = 90.0
maNL["a"] = 6.5
maNL["T0"] = 0.5
maNL["p"] = 0.16
maNL["a1"] = 0.5335
maNL["a2"] = 0.0793
maNL["a3"] = -0.3466
maNL["a4"] = 0.1609
maNL["b"] = 0.

yNLLD = maNL.evaluate(time)

# Let's compare both models...
plt.plot(time, yQLD, '-', label="Quadratic LD")
plt.plot(time, yNLLD, 'd', label="Non-linear LD")
plt.legend()
plt.show()

Non-linear limb-darkening - The MandelAgolNLLC class

Warning

The mandelAgolNL is considered deprecated. Its functionality has been absorbed in the MandelAgolLC model.

class PyAstronomy.modelSuite.XTran.forTrans.mandelAgolNL.MandelAgolNLLC

Calculate and fit analytical transit light-curves using the formulae provided by Mandel & Agol 2002. This model uses the non-linear limb-darkening prescription by Claret 2011:

\[\frac{I(\mu)}{I(1)}= 1 - \sum_{n=1}^{4}{a_n(1-\mu^{n/2})}\]

Note

The computation of transit light curves is done using the external occultnl FORTRAN library. This library must be compiled manually using SciPy’s f2py wrapper (http://www.scipy.org/F2py). Simply go to the forTrans directory of the source distribution of PyAstronomy, then invoke

f2py -c occultnl.pyf occultnl.f

Model parameters:
  • p - Radius ratio between planet and star.

  • a - Semi-major axis of planetary orbit [stellar radii].

  • i - Inclination of orbit in degrees (90 deg is edge on view).

  • a1 - Non-Linear limb-darkening coefficient.

  • a2 - Non-Linear limb-darkening coefficient.

  • a3 - Non-Linear limb-darkening coefficient.

  • a4 - Non-Linear limb-darkening coefficient.

  • T0 - Time offset of transit center.

  • per - Period of planetary orbit.

  • b - Describes the flux ratio between a stellar companion and the main star (default is 0).

Methods

MCMCautoParameters(ranges[, picky, ...])

Convenience function to generate parameters for MCMC fit.

addConditionalRestriction(*args)

Define a conditional restriction.

assignValue(specval)

Assign new values to variables.

assignValues(specval)

Assign new values to variables.

autoFitMCMC(x, y, ranges[, picky, stepsize, ...])

Convenience function to using auto-generated sampling parameters in MCMC.

availableParameters()

Provides a list of existing parameters.

delRestriction(parName)

Delete restriction

description([parenthesis])

Returns a description of the model based on the names of the individual components.

errorConfInterval(par[, dstat, statTol, ...])

Calculate confidence interval for a parameter.

evaluate(time)

Calculate a light curve according to the analytical models given by Mandel & Agol 2002.

fit(x, y[, yerr, X0, minAlgo, mAA, ...])

Carries out a fit.

fitEMCEE([x, y, yerr, nwalker, priors, ...])

MCMC sampling using emcee package.

fitMCMC(x, y, X0, Lims, Steps[, yerr, ...])

Carry out MCMC fit/error estimation.

freeParamNames()

Get the names of the free parameters.

freeParameters()

Get names and values of free parameters.

freeze(specifiers)

Consider variables free to float.

frozenParameters()

Get names and values of frozen parameters.

getRelationsOf(specifier)

Return relations of a variable.

getRestrictions()

Get all restrictions.

hasVariable(specifier)

Determine whether the variable exists.

numberOfFreeParams()

Get number of free parameters.

parameterSummary([toScreen, prefix, sorting])

Writes a summary of the parameters in text form.

parameters()

Obtain parameter names and values.

relate(dependentVar, independentVars[, func])

Define a relation.

removeConditionalRestriction(*args)

Remove an existing conditional constraint.

renameVariable(oldName, newName)

Change name of variable.

restoreState(resource)

Restores parameter values from file or dictionary.

saveState(*args, **kwargs)

Save the state of the fitting object.

setObjectiveFunction([miniFunc])

Define the objective function.

setPenaltyFactor(penalFac)

Change the penalty factor.

setRestriction(restricts)

Define restrictions.

setRootName(root[, rename])

Define the root name of the model.

showConditionalRestrictions(**kwargs)

Show conditional restrictions.

steppar(pars, ranges[, extractFctVal, quiet])

Allows to step a parameter through a specified range.

thaw(specifiers)

Consider variables fixed.

untie(parName[, forceFree])

Remove all relations of parameter parName, i.e., the parameter is not dependend on other parameters.

updateModel()

Recalculate the model using current settings.