Contents:
| FUNCTIONS: | occultuniform() – uniform-disk transit light curve occultquad() – quadratic limb-darkening occultnonlin() – full (4-parameter) nonlinear limb-darkening
|
|---|---|
| REQUIREMENTS: | |
| NOTES: | Certain values of p (<0.09, >0.5) cause some routines to hang; your mileage may vary. If you find out why, please let me know!
For occultquad() I relied heavily on the IDL code of E. Agol and J. Eastman. Function appellf1() comes from the mpmath compilation, and is adopted (with modification) for use herein in compliance with its BSD license (see function documentation for more details). |
| REFERENCE: | The main reference is that seminal work by Mandel and Agol (2002). |
| LICENSE: | Created by Ian Crossfield at UCLA. The code contained herein may be reused, adapted, or modified so long as proper attribution is made to the original authors. |
| REVISIONS: |
2011-04-25 17:32 IJMC: Fixed bug in ellpic_bulirsch().
|
Generate a JKTEBOP light curve using the F2Py-compiled library.
| INPUTS: | The inputs to the Fortran function ‘jktebop.getmodel’ are many, and their formatting is complicated. You may want to examine the Fortran source code for more insight.
|
|---|---|
| OPTIONS: | None (so far!) |
| NOTES: | See the JKTEBOP documentation for more details. JKTEBOP is currently available online at http://www.astro.keele.ac.uk/~jkt/codes/jktebop.html To compile JKTEBOP (v34) with F2Py, I had to rename the source file to be “jktebop_orig.f90”, and then I ran the following command:
|
| EXAMPLE: | import pylab as py
import jktebop_f2py
import transit
v = np.zeros(138, dtype=float)
v[[1,2,3,4,5]] = [0.211, 0.154, 0.3, 0.3, 88.59]
v[12] = 0.0013 # Mass ratio
v[17] = 1 # Integration ring size; cannot be zero.
v[18] = 1.3382320363 # Orbital period
v[19] = 54740.62 # Transit ephemeris
vary = np.zeros(138, dtype=int)
ldtype = [4, 1] # Quadratic for star, linear for planet.
nsine, npoly = 0, 0
psine, ppoly = np.zeros(9), np.zeros(9)
times = py.linspace(v[19]-0.1, v[19]+0.1, 100)
dtype1 = 1 # To compute a light curve
la, lb = 0, 0 # (????)
numint = 1 # Cannot be zero
ninterval = 0 # Irrelevant, because numint=1
magout_direct = py.array([jktebop_f2py.getmodel(v, vary, ldtype, nsine, psine, npoly, ppoly, time, dtype1, la, lb, numint, ninterval) for time in times])
magout_alt = transit.JKTEBOP_lightcurve(v, vary, ldtype, nsine, psine, npoly, ppoly, times, dtype1, la, lb, numint, ninterval)
|
| SEE_ALSO: |
vv = ((v, vary, ldtype, nsine, psine, npoly, ppoly, bjd, 1, 0, 0, numint, ninterval))
out0 = np.array([jktebop_f2py.getmodel(v, vary, ldtype, nsine, psine, npoly, ppoly, time, 1, 0, 0, numint, ninterval) for time in bjd]) out00 = jktebop_mod.getmodelarray(v, vary, ldtype, nsine, psine, npoly, ppoly, bjd, 1, 0, 0, numint, ninterval, bjd.size) out1 = transit.JKTEBOP_lightcurve(*vv) out2 = transit.JKTEBOP_lightcurve_helper(vv)
out3 = np.array(map(transit.JKTEBOP_lightcurve_helper, [(v, vary, ldtype, nsine, psine, npoly, ppoly, time0, 1, 0, 0, numint, ninterval) for time0 in bjd])).squeeze() out4 = np.array(pool.map(transit.JKTEBOP_lightcurve_helper, [(v, vary, ldtype, nsine, psine, npoly, ppoly, time0, 1, 0, 0, numint, ninterval) for time0 in bjd])).squeeze()
1e4 1e5 1e6
0 0.104 0.984 9.762 1 0.103 0.979 9.878 2 0.101 0.999 9.801 3 0.252 3.028 31.8 4 0.228 3.006 30.7
This method doesn’t seem any faster!
Determine best fit and uncertainties on transits, eclipses, phasecurves.
| SEE_ALSO: | blender.modeltransit_jktebop |
|---|
returns: bestfit, sampler, weights, bestmod
Determine best fit and uncertainties on secondary eclipse data.
| SEE_ALSO: | modeleclipse_simple |
|---|
returns: bestfit, sampler, weights, bestmod
Fit transit to data, and estimate uncertainties on the fit.
| INPUTS: |
|
|---|---|
| OUTPUTS: | (eventually, some object with useful fields) |
| SEE_ALSO: |
Give the Appell hypergeometric function of two variables.
| INPUTS: | six parameters, all scalars. |
|---|---|
| OPTIONS: | eps – scalar, machine tolerance precision. Defaults to 1e-10. |
| NOTES: | Adapted from the mpmath module, but using the scipy (instead of mpmath) Gauss hypergeometric function speeds things up. |
| LICENSE: | MPMATH Copyright (c) 2005-2010 Fredrik Johansson and mpmath contributors. All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
Implement a simple/stupid version of the Box-Least-Squares algorithm.
| INPUTS: |
|
|---|---|
| OUTPUTS: | (test_periods, reduction_in_dispersion) |
| NOTES: | Transits shorter than (prange[0]/nbins) and longer than (maxwid*prange[1]/nbins) may not be correctly modelled. This routine is my own crude attempt at a box-fitting least-squares algorithm. Note however that this is NOT the well-publicized (and more rigorous) version of Kovacs et al. (2002). Here for each trial period, I bin the data and construct simple box-shaped transit models for combinations of ingress and egress times. For each model the reduction in the standard deviation of (binned data - model) is calculated; for each period, the greatest reduction value is reported. |
REturn a boolean mask, True wherever the planet is in transit.
Create an input file suitable for running JKTEBOP.
| INPUT: | A long sequence of numerical & string inputs, defined as follows. If simulating a transiting planet, object “A” indicates the star and object “B” indicates the planet.
|
|---|---|
| OPTIONS: |
|
| OUTPUT: | If no ‘filename’ option was passed in, returns a list of strings, suitable to writing to disk via the usual:
|
| EXAMPLE: | vals = [2, 1, 0.21, 0.15, 88.5, 0.0013, 0, 0, 1, 1, 0, 0, 'quad', 'lin', 0.3, 0, 0.3, 0, 0, 0, 0, 0.6]
output = transit.createJKTEBOPinput(vals, filename='test.in', clobber=False)
|
| NOTES: | Put a negative number for the mass ratio to force the stars to be spherical. The mass ratio will then be irrelevant (it is only used to get deformations). To input R_A/a and R_B/a (instead of [R_A+R_B]/a and R_B/R_A), give a negative value for [R_A+R_B]/a. Then it will be interpreted to mean R_A/a, and R_B/R_A will be interpreted as R_B/a. If eccentricity < 10 then e and omega will be assumed to be e*cos(omega) and e*sin(omega). If e >= 10 then e and omega will be assumed to be (e+10) and omega (degrees). The first option is often better unless eccentricity is larger or fixed. See the JKTEBOP documentation for more details on all these parameters. JKTEBOP is currently available online at http://www.astro.keele.ac.uk/~jkt/codes/jktebop.html |
| TO_DO: | Add other optional parameters: TMIN, LRAT, THDL, ECSW, ENSW, SINE, POLY, NUMI, RV1 & RV2, orbital period, reference epoch... Allow fitting (i.e., enable Tasks 3-9). |
Compute Hasting’s polynomial approximation for the complete elliptic integral of the first (ek) and second (kk) kind.
| INPUTS: | k – scalar or Numpy array |
|---|---|
| OUTPUTS: | ek, kk |
| NOTES: | Adapted from the IDL function of the same name by J. Eastman (OSU). |
Compute complete elliptic integrals of the first kind (K) and second kind (E) using the series expansions.
Compute the complete elliptical integral of the third kind using the algorithm of Bulirsch (1965).
| INPUTS: | n – scalar or Numpy array k– scalar or Numpy array |
|---|---|
| NOTES: | Adapted from the IDL function of the same name by J. Eastman (OSU). |
Evaluate the integral at a specified limit (upper or lower)
data: time series to fit using least-squares.
sv: state vectors (e.g., various instrumental parameters)
ords: orders to raise each sv vector to: e.g., [1, [1,2], 3]
tlc: eclipse light curve
edata: error on the data (for chisq ONLY! No weighted fits.)
index: array index to apply to data, sv, and tlc
dopb: do prayer-bead uncertainty analysis
dotransit: include tlc in the fitting; otherwise, leave it out.
Helper function for bls_simple.
args = thisperiod, times, flux, phasebins, maxwid):
Return the integral in I*(z) in Eqn. 8 of Mandel & Agol (2002). – Int[I(r) 2r dr]_{z-p}^{1}, where:
| INPUTS: |
p = scalar. Planet/star radius ratio.
lower, upper – floats. Limits of integration in units of mu |
|---|---|
| RETURNS: | value of the integral at specified z. |
MCMC for 3-parameter eclipse function with KNOWN orbit
| INPUTS: |
|
|---|---|
| RETURNS: |
|
| REFERENCES: | Numerical Recipes, 3rd Edition (Section 15.8); Wikipedia |
MCMC for 5-parameter eclipse function of transit with KNOWN period
| INPUTS: |
|
|---|---|
| RETURNS: |
|
| REFERENCES: | Numerical Recipes, 3rd Edition (Section 15.8); Wikipedia |
MCMC for 5-parameter eclipse function of transit with KNOWN period
| INPUTS: |
|
|---|---|
| RETURNS: |
|
| REFERENCES: | Numerical Recipes, 3rd Edition (Section 15.8); Wikipedia |
Model an eclipse light curve of arbitrary type to a flux time series, assuming zero eccentricity and a fixed, KNOWN period.
| INPUTS: |
func – function to fit to data; presumably transit.occultuniform()
|
|---|
Model an eclipse light curve of arbitrary type to a flux time series, assuming zero eccentricity and a fixed, KNOWN orbit.
| INPUTS: |
func – function to fit to data; presumably transit.occultuniform() t – numpy array. Time of observations. |
|---|
Model an eclipse light curve of arbitrary type to a flux time series, assuming zero eccentricity and a fixed, KNOWN orbit.
| INPUTS: |
func – function to fit to data; presumably transit.occultuniform()
|
|---|
Model a full planetary light curve: transit, eclipse, and (sinusoidal) phase variation. Accept independent eclipse and transit times-of-center, but otherwise assume a circular orbit (and thus symmetric transits and eclipses).
| INPUTS: | params – (M+10+N)-sequence with the following:
|
|---|---|
| OPTIONS: |
|
| EXAMPLE: | TBW |
| NOTES: | This should be updated to use the new ‘transitonly’ options in t2z() |
Model a transit light curve of arbitrary type to a flux time series, assuming zero eccentricity and a fixed, KNOWN period.
| INPUTS: |
func – function to fit to data, e.g. transit.occultquad per – float. Orbital period, in days. t – numpy array. Time of observations. |
|---|
Model a transit light curve of arbitrary type to a flux time series, assuming zero eccentricity and a fixed, KNOWN period, and assuming MIPS-type data with 14 separate sensitivity dependencies.
| INPUTS: | params – (14+5+N)-sequence with the following:
func – function to fit to data, e.g. transit.occultquad per – float. Orbital period, in days.
|
|---|---|
| SEE ALSO: |
Model a transit light curve of arbitrary type to a flux time series, assuming zero eccentricity.
| INPUTS: |
t – numpy array. Time of observations.
|
|---|---|
| NOTES: | If quadratic or linear limb-darkening (L.D.) is used, the sum of the L.D. coefficients cannot exceed 1. If they do, this routine normalizes the coefficients [g1,g2] to: g_i = g_i / (g1 + g2). If “Rp/R*”, or “R*/a” are < 0, they will be set to zero. If “P” < 0.01, it will be set to 0.01. If “inc” > 90, it will be set to 90. |
Add a simple ad-hoc haze model to a planet’s radius spectrum.
| INPUTS: |
|
|---|---|
| NOTES: | Because this routine makes use of the numpy.interp function, using a smaller input model grid can significantly speed things up. |
Nonlinear limb-darkening light curve; cf. Section 3 of Mandel & Agol (2002).
| INPUTS: | z – sequence of positional offset values p0 – planet/star radius ratio cn – four-sequence. nonlinear limb darkening coefficients |
|---|---|
| EXAMPLE: | # Reproduce Figure 2 of Mandel & Agol (2002):
from pylab import *
import transit
z = linspace(0, 1.2, 50)
cns = vstack((zeros(4), eye(4)))
figure()
for coef in cns:
f = transit.occultnonlin(z, 0.1, coef)
plot(z, f)
|
| SEE ALSO: | |
| NOTES: | Scipy is much faster than mpmath for computing the Beta and Gauss hypergeometric functions. However, Scipy does not have the Appell hypergeometric function – the current version is not vectorized. |
Nonlinear limb-darkening light curve in the small-planet approximation (section 5 of Mandel & Agol 2002).
| INPUTS: | z – sequence of positional offset values p – planet/star radius ratio
|
|---|---|
| NOTE: | I had to divide the effect at the near-edge of the light curve by pi for consistency; this factor was not in Mandel & Agol, so I may have coded something incorrectly (or there was a typo). |
| EXAMPLE: | # Reproduce Figure 2 of Mandel & Agol (2002):
from pylab import *
import transit
z = linspace(0, 1.2, 100)
cns = vstack((zeros(4), eye(4)))
figure()
for coef in cns:
f = transit.occultnonlin_small(z, 0.1, coef)
plot(z, f, '--')
|
| SEE ALSO: |
Quadratic limb-darkening light curve; cf. Section 4 of Mandel & Agol (2002).
| INPUTS: | z – sequence of positional offset values p0 – planet/star radius ratio
|
|---|---|
| OPTIONS: |
|
| EXAMPLE: | # Reproduce Figure 2 of Mandel & Agol (2002):
from pylab import *
import transit
z = linspace(0, 1.2, 100)
gammavals = [[0., 0.], [1., 0.], [2., -1.]]
figure()
for gammas in gammavals:
f = transit.occultquad(z, 0.1, gammas)
plot(z, f)
# Calculate the same geometric transit with two different
# sets of limb darkening coefficients:
from pylab import *
import transit
p, b = 0.1, 0.5
x = (arange(300.)/299. - 0.5)*2.
z = sqrt(x**2 + b**2)
gammas = [.25, .75]
F1, Funi, lambdad, etad = transit.occultquad(z, p, gammas, retall=True)
gammas = [.35, .55]
F2 = 1. - ((1. - gammas[0] - 2.*gammas[1])*(1. - F1) +
(gammas[0] + 2.*gammas[1])*(lambdad + 2./3.*(p > z)) + gammas[1]*etad) /
(1. - gammas[0]/3. - gammas[1]/6.)
figure()
plot(x, F1, x, F2)
legend(['F1', 'F2'])
|
| SEE ALSO: | |
| NOTES: | In writing this I relied heavily on the occultquad IDL routine by E. Agol and J. Eastman, especially for efficient computation of elliptical integrals and for identification of several apparent typographic errors in the 2002 paper (see comments in the source code). From some cursory testing, this routine appears about 9 times slower than the IDL version. The difference drops only slightly when using precomputed quantities (i.e., retall=True). A large portion of time is taken up in ellpic_bulirsch() and ellke(), but at least as much is taken up by this function itself. More optimization (or a C wrapper) is desired! |
Uniform-disk transit light curve (i.e., no limb darkening).
| INPUTS: |
p – scalar; planet/star radius ratio.
|
|---|---|
| SEE ALSO: |
Simple toy model for Pixel-Level-Decorrelation testing.
tparams and time are for modeleclipse_simple()
vecs are: np.vstack((phat.T, othervecs.T)); of shape (Nvec x Nobs)
params are [eclipseCenter, ... then all the PLD coefs]
Run JKTEBOP simulation from the command line, and return results.
| INPUTS: | See createJKTEBOPinput() for a description of the inputs. For now, only “Task 2” can be run through this Python interface. |
|---|---|
| OPTIONS: |
|
| NOTES: | See the JKTEBOP documentation for more details. JKTEBOP is currently available online at http://www.astro.keele.ac.uk/~jkt/codes/jktebop.html |
| EXAMPLE: | vals = [2, 1, 0.21, 0.15, 88.5, 0.0013, 0, 0, 1, 1, 0, 0, 'quad', 'lin', 0.3, 0, 0.3, 0, 0, 0, 0, 0.6]
inputFile = 'JKTEBOP_test.in'
outputFile = 'JKTEBOP_test.out'
period = 9.8
exe = os.path.expanduser('~')+'/proj/transit/jktebop/jktebop_orig'
time, mag = transit.runJKTEBOP(vals, infile=inputFile, outfile=outputFile, period=period, clobber=True, exe=exe)
## THIS IS A SPECIAL EXAMPLE FOR MY CODE ONLY --
## I HAVE MODIFIED 'jktebop' SO DON'T EXPECT THE FOLLOWING TO
## WORK PROPERLY ON YOUR COMPUTER!
import transit
import pylab as py
vals = [2, 1, 0.211, 0.154, 88.59, 0.0013, 0, 0, 1, 1, 0, 0, 'quad', 'lin', 0.3, 0., 0.3, 0, 0, 0, 0, 0.0]
inputFile = 'JKTEBOP_test_mod.in'
outputFile = 'JKTEBOP_test_mod.out'
dataFile = 'wasp4.dat' # data file; only its timestamps are used!
period = 1.3382320363 # Orbital period, in days
t0 = 54740.62 # Time of central transit, in days.
exe = os.path.expanduser('~')+'/proj/transit/jktebop/jktebop_mod'
time, mag = transit.runJKTEBOP(vals, infile=inputFile, outfile=outputFile, datfile=dataFile, period=period, t0=t0, clobber=True, exe=exe)
exampleData = py.loadtxt(dataFile)
examplePhase = (exampleData[:,0] - t0) / period
py.figure()
py.plot(examplePhase, 10**(-0.4*exampleData[:,1]), 'ob')
py.plot(time/period, 10**(-0.4*mag), '--r', linewidth=2)
py.xlabel('Orbital Phase')
py.ylabel('Normalized Flux')
py.minorticks_on()
py.legend(['Sample Observations', 'JKTEBOP Model'], 4)
|
| SEE_ALSO: | |
| TO_DO: | Enable light-curve simulations (Tasks 3-9). |
Placeholder for backwards compatibility with my old code. The function is now called occultnonlin_small().
Convert HJD (time) to transit crossing parameter z.
| INPUTS: | tt – scalar. transit ephemeris per – scalar. planetary orbital period (in days) inc – scalar. orbital inclination (in degrees)
ars – scalar. ratio a/Rs, orbital semimajor axis over stellar radius ecc – scalar. orbital eccentricity. longperi=0 scalar. longitude of periapse (in radians)
|
|---|---|
| ALGORITHM: | At zero eccentricity, z relates to physical quantities by: z = (a/Rs) * sqrt(sin[w*(t-t0)]**2+[cos(i)*cos(w*[t-t0])]**2) |
Placeholder for my old code; the new function is called occultuniform().
Convert transit crossing parameter z to a time offset for circular orbits.
| INPUTS: | per – scalar. planetary orbital period inc – scalar. orbital inclination (in degrees) ars – scalar. ratio a/Rs, orbital semimajor axis over stellar radius z – scalar or array; transit crossing parameter z. |
|---|---|
| RETURNS: |
|