Planetary phase curve routines

Contents:

Functions for fitting phase curves.

Use of ‘errfunc’ is encouraged.

REQUIREMENTS:Transit light curve routines

phasecurves.devfunc(*arg, **kw)

Generic function to give the weighted residuals on a function or functions:

INPUTS:

(fitparams, function, arg1, arg2, ... , depvar, weights)

OR:

(fitparams, function, arg1, arg2, ... , depvar, weights, kw)

OR:

(allparams, (args1, args2, ..), npars=(npar1, npar2, ...))

where allparams is an array concatenation of each functions input parameters.

If the last argument is of type dict, it is assumed to be a set of keyword arguments: this will be added to resfunc’s direct keyword arguments, and will then be passed to the fitting function **kw. This is necessary for use with various fitting and sampling routines (e.g., kapteyn.kmpfit and emcee.sampler) which do not allow keyword arguments to be explicitly passed. So, we cheat! Note that any keyword arguments passed in this way will overwrite keywords of the same names passed in the standard, Pythonic, way.

OPTIONAL INPUTS:  

jointpars – list of 2-tuples.

For use with multi-function calling (w/npars keyword). Setting jointpars=[(0,10), (0,20)] will always set params[10]=params[0] and params[20]=params[0].

gaussprior – list of 2-tuples (or None values), same length as “fitparams.”

The i^th tuple (x_i, s_i) imposes a Gaussian prior on the i^th parameter p_i by adding ((p_i - x_i)/s_i)^2 to the total chi-squared. Here in devfunc(), we _scale_ the error-weighted deviates such that the resulting chi-squared will increase by the desired amount.

uniformprior – list of 2-tuples (or ‘None’s), same length as “fitparams.”

The i^th tuple (lo_i, hi_i) imposes a uniform prior on the i^th parameter p_i by requiring that it lie within the specified “high” and “low” limits. We do this (imprecisely) by multiplying the resulting deviates by 1e9 for each parameter outside its limits.

ngaussprior – list of 3-tuples of Numpy arrays.

Each tuple (j_ind, mu, cov) imposes a multinormal Gaussian prior on the parameters indexed by ‘j_ind’, with mean values specified by ‘mu’ and covariance matrix ‘cov.’ This is the N-dimensional generalization of the ‘gaussprior’ option described above. Here in devfunc(), we _scale_ the error-weighted deviates such that the resulting chi-squared will increase by the desired amount.

For example, if parameters 0 and 3 are to be jointly constrained (w/unity means), set:

jparams = np.array([0, 3]) mu = np.array([1, 1]) cov = np.array([[1, .9], [9., 1]]) ngaussprior=[[jparams, mu, cov]] # Double brackets are key!

EXAMPLE:

from numpy import *
import phasecurves
def sinfunc(period, x): return sin(2*pi*x/period)
snr = 10
x = arange(30.)
y = sinfunc(9.5, x) + randn(len(x))/snr
guess = 8.
period = optimize.fmin(phasecurves.errfunc,guess,args=(sinfunc,x, y, ones(x.shape)*snr**2))

phasecurves.domodfit(profile, fdatc, wmat, xmat)

Helper function for fitsin

Generates takes a matrix of variables and adds (1) a phase profile to the first row and (2) a flux conservation constraint to the last column.

phasecurves.doubleexp(params, t)

Model Agol et al. 2010’s double-exponential IRAC ramp function.

params - 4- or 5-sequence, defining the function as shown below.

t - sequence. Input time (presumed: since start of observations)

Functional form:

if len(params)==4:

F’/F = 1 - p0*exp(-t/p2) - p1*exp(-t/p3)

elif len(params)>4:

F’/F = p4 * (1 - p0*exp(-t/p2) - p1*exp(-t/p3))

phasecurves.doubleexp2(params, t)

Model a different double-exponential ramp function.

INPUTS:

params : 3- or 4-sequence

Parameters that define the function as shown below.

t : sequence.

Input time (presumed: since start of observations)

Functional_form:  

if len(params)==3:

F’/F = (1 - p0*exp(-t/p1)) * exp(-t/p2)

elif len(params)>3:

F’/F = (1 - p0*exp(-t/p1)) * exp(-t/p2) * p3

phasecurves.doubleexp214(params, t)

Model a different double-exponential ramp function w/14 positional offsets.

INPUTS:

params : sequence of length [(3 or 4) + 14]

Parameters that define the function as shown below.

t : sequence.

Input time (presumed: since start of observations). If not of shape (14 x N), will be reshaped to that.

Functional_form:  

if len(params)==3:

F = (1 - p0*exp(-t/p1)) * exp(-t/p2)

elif len(params)>3:

F = (1 - p0*exp(-t/p1)) * exp(-t/p2) * p3

return (F * (1 + p[3/4::]).reshape(14,1) )

phasecurves.eclipse14_single(param, tparam, t, xyord=None, x=None, y=None)

compute 3-parameter eclipse function of a single event, and account for 14 different possible flux offsets and X/Y positional motions.

param : parameters to be fit

[Fstar, Fplanet, t_center, c0, c1, ... c13]

tparam : parameters to be held constant (from transit)

[b, v (in Rstar/day), p (Rp/Rs)]

Input data ‘t’,’x’,’y’ must therefore be of size (14xN); if not, it will be reshaped into that.

If xyord==None or 0:

f[i, j] = (1. + p[3+i]) * [eclipse light curve]

“xyord” determines the linear order of the polynomial in x and y. If xyord==1, then: f[i,j] = (...) * (1.+ p[3+i] +p[17+i]*x + p[31+i]*y)

[note that the the individual offsets will be subjected to the constraint: param[3] = 1./(1.+param[4:17]).prod() - 1. ]

phasecurves.eclipse_single(param, tparam, t)

compute 3-parameter eclipse function of a single event

param : 3 parameters to be fit

[Fstar, Fplanet, t_center]

tparam : 3 parameters to be held constant (from transit)

[b, v (in Rstar/day), p (Rp/Rs)]

Input data ‘t’ must be of size (14xN); if not, it will be reshaped

into that.

phasecurves.errfunc(*arg, **kw)

Generic function to give the chi-squared error on a generic

function or functions:

INPUTS:

(fitparams, function, arg1, arg2, ... , depvar, weights)

OR:

(fitparams, function, arg1, arg2, ... , depvar, weights, kw)

OR:

(allparams, (args1, args2, ..), npars=(npar1, npar2, ...))

where allparams is an array concatenation of each functions input parameters.

If the last argument is of type dict, it is assumed to be a set of keyword arguments: this will be added to errfunc2’s direct keyword arguments, and will then be passed to the fitting function **kw. This is necessary for use with various fitting and sampling routines (e.g., kapteyn.kmpfit and emcee.sampler) which do not allow keyword arguments to be explicitly passed. So, we cheat! Note that any keyword arguments passed in this way will overwrite keywords of the same names passed in the standard, Pythonic, way.

OPTIONAL INPUTS:  

jointpars – list of 2-tuples.

For use with multi-function calling (w/npars keyword). Setting jointpars=[(0,10), (0,20)] will always set params[10]=params[0] and params[20]=params[0].

gaussprior – list of 2-tuples (or None values), same length as “fitparams.”

The i^th tuple (x_i, s_i) imposes a Gaussian prior on the i^th parameter p_i by adding ((p_i - x_i)/s_i)^2 to the total chi-squared. Here in devfunc(), we _scale_ the error-weighted deviates such that the resulting chi-squared will increase by the desired amount.

uniformprior – list of 2-tuples (or ‘None’s), same length as “fitparams.”

The i^th tuple (lo_i, hi_i) imposes a uniform prior on the i^th parameter p_i by requiring that it lie within the specified “high” and “low” limits. We do this (imprecisely) by multiplying the resulting deviates by 1e9 for each parameter outside its limits.

ngaussprior – list of 3-tuples of Numpy arrays.

Each tuple (j_ind, mu, cov) imposes a multinormal Gaussian prior on the parameters indexed by ‘j_ind’, with mean values specified by ‘mu’ and covariance matrix ‘cov.’ This is the N-dimensional generalization of the ‘gaussprior’ option described above. Here in devfunc(), we _scale_ the error-weighted deviates such that the resulting chi-squared will increase by the desired amount.

For example, if parameters 0 and 3 are to be jointly constrained (w/unity means), set:

jparams = np.array([0, 3]) mu = np.array([1, 1]) cov = np.array([[1, .9], [9., 1]]) ngaussprior=[[jparams, mu, cov]] # Double brackets are key!

scaleErrors – bool

If True, fit for the measurement uncertainties. In this case, the first element of ‘fitparams’ (“s”) is used to rescale the measurement uncertainties. Thus weights –> weights/s^2, and chi^2 –> 2 N log(s) + chi^2/s^2 (for N data points).

EXAMPLE:

from numpy import *
import phasecurves
def sinfunc(period, x): return sin(2*pi*x/period)
snr = 10
x = arange(30.)
y = sinfunc(9.5, x) + randn(len(x))/snr
guess = 8.
period = optimize.fmin(phasecurves.errfunc,guess,args=(sinfunc,x, y, ones(x.shape)*snr**2))

phasecurves.errfunc14xymult_cfix(*arg, **kw)

Generic function to give the chi-squared error on a generic function:

INPUTS:

(fitparams, arg1, arg2, ... indepvar, depvar, weights)

phasecurves.fitcurve(pdat0, fdat0, efdat0, i2, xmat, phi=30, bsind=None, func=<ufunc 'cos'>, retall=False, args=None, nsigma=inf)

Decorrelate data with 14 XYs, 14 offsets, and a curve model of arbitrary amplitude.

INPUTS: phase data (14xN) flux data (14xN) error on flux data (14xN) boolean time index array (14) matrix generated for testing

phi – phi values to test bsind – (N). bootstrap indices. Indices, 0-N inclusive and with

repetitions allowed.

func – phase function, returning amplitude for an input phase

value in radians.

phasecurves.fithelper(param, func, pdat0, fdat0, efdat0, fdatc, wmatc, xmat, i2, i3, bsind=None, nsigma=5, retall=False)

Helper funtion for fitcurve – XXX

Param – either [phi_0] or [phi_0, inclination]

phasecurves.fitsin(pdat0, fdat0, efdat0, i2, xmat, phi=30, bsind=None)

Decorrelate data with 14 XYs, 14 offsets, and a sinusoid model.

INPUTS: phase data (14xN) flux data (14xN) error on flux data (14xN) boolean time index array (14) matrix generated for testing

phi – either phi values to test, or number of evenly space phi

values to test.

bsind – (N). bootstrap indices. Indices, 0-N inclusive and with

repetitions allowed.

phasecurves.lam4fit(param, x)

Compute labertian phase function with a fixed period=1, and x in units of orbital cycles.

param = [DC pedestal, AC amplitude, inclination (radians)]

phasecurves.lam4fit2(param, x)

Compute labertian phase function with a fixed period=1, and x in units of orbital cycles – this time with a variable phase offset.

param = [DC pedestal, AC amplitude, inclination (radians), phase offset (radians)]

phasecurves.lam4fit_noinc(param, inc, x)

Compute labertian phase function with a fixed period=1, and x in units of orbital cycles – a variable phase offset but FIXED inclination (in radians)

param = [DC pedestal, AC amplitude, phase offset (radians)]

phasecurves.lambertian(ophase, inc=1.5707963267948966)

Return a lambertian phase function with peak-to-valley amplitude unity.

INPUTS:

ophase (seq) an orbital phase (in radians). Secondary eclipse

(or ‘opposition’ for non-transiting planets) occurs at pi; Primary transit (or ‘conjuction’) occurs at 0 or 2*pi

inc (float) system inclination angle (also in radians). Edge-on

is pi/2, face-on is 0.

phasecurves.lambertian_amplitude(inc, n=5000)

Return amplitude of a nominally unity-amplitude lambertian with a given inclination angle, using function ‘lambertian’. inc is in radians.

phasecurves.lambertian_mean(inc, n=5000)

Return mean of a nominally unity-amplitude lambertian with a given inclination angle, using function ‘lambertian’. inc is in radians.

phasecurves.lnprobfunc(*arg, **kw)

Return natural logarithm of posterior probability: i.e., -chisq/2.

Inputs are the same as for errfunc().

SEE ALSO:gaussianprocess.negLogLikelihood()

phasecurves.makexmat(xpos, ypos, constraint=True)

Generate matrix for least-squares MIPS photometric detrending.

Generate the LSq dependent-variable matrix. The rows are: 0–constant pedestal, (ii*3+1)–dither position-dependent pedestal offset, (ii*3+2)–dither position-dependent x-position correlation, (ii*3+3)–dither position-dependent y-position correlation.

EXAMPLE:

xmat = makexmat() cleanedphot = dot(coef, xmat)

phasecurves.map2lightcurve(map, z, k, cumulative_ingress_masks, cumulative_egress_masks, ingress_zmasks, egress_zmasks, alt=False)

Take a 2D planet map, and convert it into an eclipse light curve.

INPUTS:

map : 2D Numpy array

[M x M] square map of planet surface brightness distribution. Note that map.sum() corresponds to the fractional eclipse depth.

z : 1D Numpy array

length-N planet crossing parameter z (i.e., distance between planet and stellar geocenters in units of stellar radii).

k : scalar

planet/star radius ratio

cumulative_ingress_masks, cumulative_egress_masks : 3D Numpy arrays

[nslice x M x M] boolean maps of planet occultation at various stages, from prepEclipseMap()

ingress_zmasks, egress_zmasks : 2D Numpy arrays

[nslice x N] boolean maps of which slices correspond to which z-indices, from prepEclipseMap()

RETURNS:

lightcurve : 1D Numpy array

length-N light curve, normalized to unity in eclipse.

phasecurves.mcmc14xy(z, t, x, y, sigma, params, stepsize, numit, nstep=1, sumtol=1e+20, xyord=1, prodtol=1e-06, fixc0=False)

Applies Markov Chain Monte Carlo model fitting using the Metropolis-Hastings algorithm.

INPUTS:

z : 1D array

Contains dependent data

t,x,y : 1D array

Contains independent data: phase, x- and y-positions

sigma : 1D array

Standard deviation of dependent (y) data

params : 1D array

Initial guesses for parameters

stepsize : 1D array

Array of 1-sigma change in parameter per iteration

numit : int

Number of iterations to perform

nstep : int

Saves every “nth” step of the chain

sumtol : float

Tolerance error, used as: abs(params[3:17].sum()/sumtol)**2

xyord : int

“xyord” determines the linear order of the polynomial in x and y.

Returns:

allparams : 2D array

Contains all parameters at each step

bestp : 1D array

Contains best paramters as determined by lowest Chi^2

numaccept: int

Number of accepted steps

chisq: 1D array

Chi-squared value at each step

REFERENCES:

Numerical Recipes, 3rd Edition (Section 15.8); Wikipedia

REVISIONS:

2008-05-02 Kevin Stevenson, UCF

kevin218@knights.ucf.edu Started converting MCMC from IDL to Python while making upgrades

2008-06-21 Kevin Stevenson

Finished updating, current stable version

2008-11-17 Kevin Stevenson

Updated docstring Simplified routine for AST3110 project and AST5765 class demo

2010-06-12 14:21 IJC: Modified for phasecurve test

2010-06-30 11:43 IJC: Added nstep option

2010-07-01 16:47 IJC: Added sum constraint for parameters 3-17.

2010-07-14 11:26 IJC: Added product constraints for parameters 3-17

phasecurves.mcmc_eclipse14_single(z, t, x, y, sigma, params, tparam, stepsize, numit, nstep=1, xyord=None)

MCMC for 17-parameter eclipse function of a single event

INPUTS:

z : 1D array

Contains dependent data

t,x,y : 1D array

Contains independent data: phase, x- and y-positions

sigma : 1D array

Contains standard deviation of dependent (z) data

params : 17 parameters to be fit

[Fstar, Fplanet, t_center, c0, ... , c13]

tparam : 3 parameters to be held constant (from transit)

[b, v (in Rstar/day), p (Rp/Rs)]

stepsize : 1D array

Array of 1-sigma change in parameter per iteration

numit : int

Number of iterations to perform

nstep : int

Saves every “nth” step of the chain

xyord : int

Highest order in x/y motions

RETURNS:

allparams : 2D array

Contains all parameters at each step

bestp : 1D array

Contains best paramters as determined by lowest Chi^2

numaccept: int

Number of accepted steps

chisq: 1D array

Chi-squared value at each step

REFERENCES:

Numerical Recipes, 3rd Edition (Section 15.8); Wikipedia

phasecurves.mcmc_eclipse_single(z, t, sigma, params, tparam, stepsize, numit, nstep=1)

MCMC for 3-parameter eclipse function of a single event

INPUTS:

z : 1D array

Contains dependent data

t : 1D array

Contains independent data: phase, x- and y-positions

sigma : 1D array

Contains standard deviation of dependent (z) data

params : 3 parameters to be fit

[Fstar, Fplanet, t_center]

tparam : 3 parameters to be held constant (from transit)

[b, v (in Rstar/day), p (Rp/Rs)]

stepsize : 1D array

Array of 1-sigma change in parameter per iteration

numit : int

Number of iterations to perform

nstep : int

Saves every “nth” step of the chain

RETURNS:

allparams : 2D array

Contains all parameters at each step

bestp : 1D array

Contains best paramters as determined by lowest Chi^2

numaccept: int

Number of accepted steps

chisq: 1D array

Chi-squared value at each step

REFERENCE:

Numerical Recipes, 3rd Edition (Section 15.8); Wikipedia

phasecurves.model_fixed_param(varparam, fixedparam, fixedindex, func, *arg, **kw)

Allow modeling with some parameters held constant.

INPUTS:

varparam : sequence

Primary parameters (which can be varied).

fixedparam : sequence

Secondary parameters (which should be held fixed.)

fixedindex : sequence of ints

Indices of parameters which should be held fixed, when passed to

func : function

Modeling function. Arguments

Thus if param = [10, 20, 50, 70] and holdfixed = [2], one would set varparam = [10, 50, 70] and fixedparam = [20].

OPTIONS:

*arg : tuple

Arguments to be passed to func

**kw : dict

Keyword (optional) arguments to be passed to func

OUTPUTS:

func(param, *arg, **kw)

phasecurves.nsin(param, x)

Compute phase function with unknown period, assuming n sinusoids:

p(x) = param[0] - param[i+1] * cos(2*pi*x/param[i+2] + param[i+3]) - param[i+3+1] * cos(2*pi*x/param[i+3+2] + param[i+3+3]) - ...

phasecurves.phase2sin(param, x, absamp=True)

compute double-sinusoid with a fixed period=1 offset (in radians):

p(x) = param[0] - param[1] * cos(2*pi*x + param[2]) - param[3] * cos(4*pi*x + param[4])

ampsamp: if True, take the absolute value of param[1] & param[3]

phasecurves.phaselamb14(param, x)

compute phase function with a fixed period=1, assuming a sinusoid, and account for 14 different possible flux offsets.

Input data ‘x’ must therefore be of size (14xN); if not, it will be reshaped into that.

p[i,j] = (1. + param[3+j]) * (param[0] + param[1]*lambertian(2*pi*x[i,j]+param[2]))

[note that the first c-parameter (individual offset) will be

constrained such that: prod(1. + param[3::]) = 1.]

phasecurves.phaselinsin(param, x)

compute phase function with a linear drift and period=1:

p(x) = param[0]+param[1]*x - param[2]*cos(2*pi*x + param[3])

phasecurves.phaselinsin2(param, phaseoffset, x)

compute phase function with a linear drift and fixed period=1 and fixed phase offset (in radians):

p(x) = param[0]+param[1]*x - param[2]*cos(2*pi*x + phaseoffset)

phasecurves.phasepoly14(param, x)

compute phase function, assuming a polynomial, and account for 14 different possible flux offsets.

Input data ‘x’ must be in units of orbital phase, and must be of size (14xN); if not, it will be reshaped into that.

For an order-N polynomial:

p[i,j] = (1. + param[N+j]) * (numpy.polyval(param[0:N], x))

[note that the first c-parameter (individual offset) will be

constrained such that: prod(1. + param[N::]) = 1.]

phasecurves.phasesin(param, x)

compute phase function with a fixed period=1, assuming a sinusoid:

p(x) = param[0] - param[1] * cos(2*pi*x + param[2])

phasecurves.phasesin14(param, x)

compute phase function with a fixed period=1, assuming a sinusoid, and account for 14 different possible flux offsets.

Input data ‘x’ must therefore be of size (14xN); if not, it will be reshaped into that.

p[i,j] = (1. + param[3+j]) * (param[0] - param[1]*cos(2*pi*x[i,j]+param[2]))

[note that the first c-parameter (individual offset) will be

constrained such that: prod(1. + param[3::]) = 1.]

phasecurves.phasesin14xymult(param, xyord, crossord, t, x, y)

compute phase function with a fixed period=1, assuming a sinusoid, and account for 14 different possible flux offsets and X/Y positional motions.

Input data ‘t’,’x’,’y’ must therefore be of size (14xN); if not, it will be reshaped into that.

“xyord” determines the linear order of the polynomial in x and y. If xyord==1, then: f[i,j] = p[0] - p[1]*cos(2*pi*x[i,j]+p[2]) + p[3+i] +p[17+i]*x + p[31+i]*y

If xyord==2, then: f[i,j] = p[0] + p[3+i] +p[17+i]*x + p[31+i]*y +p[45+i]*x**2 + p[59+i]*y**2

• p[1]*cos(2*pi*x[i,j]+p[2])

If crossord==1, then the cross-terms (x*y) will be included using

14 coefficients. crossord>1 still just uses the single cross-terms; I haven’t generalized it yet to higher orders.

[note that the the individual offsets will be subjected to the constraint: param[3::] -= param[3::].mean() ]

phasecurves.phasesin14xymult_cfix(param, xyord, crossord, t, x, y)

compute phase function with a fixed period=1, assuming a sinusoid, and account for 14 different possible flux offsets and X/Y positional motions.

Input data ‘t’,’x’,’y’ must therefore be of size (14xN); if not, it will be reshaped into that.

“xyord” determines the linear order of the polynomial in x and y. If xyord==1, then: f[i,j] = (p[0] - p[1]*cos(2*pi*x[i,j]+p[2])) * (1.+ p[3+i] +p[17+i]*x + p[31+i]*y)

If xyord==2, then: f[i,j] = (p[0] - p[1]*cos(2*pi*x[i,j]+p[2])) * (1.+ p[3+i] +p[17+i]*x + p[31+i]*y)

f[i,j] = (p[0] - p[1]*cos(2*pi*x[i,j]+p[2])) * (1.+ p[3+i] +p[17+i]*x + p[31+i]*y +p[45+i]*x**2 + p[59+i]*y**2

If crossord==1, then the cross-terms (x*y) will be included using

14 coefficients. crossord>1 still just uses the single cross-terms; I haven’t generalized it yet to higher orders.

[note that the the individual offsets will be subjected to the constraint: param[3] = 1./(1.+param[4:17]).prod() - 1. ]

phasecurves.phasesin2(param, phaseoffset, x)

compute phase function with a fixed period=1 and phase offset (in radians):

p(x) = param[0] - param[1] * cos(2*pi*x + phaseoffset)

phasecurves.phasesinsin14(param, x)

compute phase function with a fixed period=1, assuming a sinusoid and first harmonic, and account for 14 different possible flux offsets.

Input data ‘x’ must therefore be of size (14xN); if not, it will be reshaped into that.

p[i,j] = (1. + param[5+j]) * [ param[0] - param[1]*cos(2*pi*x[i,j]+param[2]) + param[3]*cos(4*pi*x[i,j]+param[4]) ]

[note that the first c-parameter (individual offset) will be

constrained such that: prod(1. + param[5::]) = 1.]

NOTE:The magnitude of the amplitudes will always be taken; they cannot be negative.

phasecurves.phasesinsin14_2(param, x)

compute phase function with a fixed period=1, assuming a sinusoid and first harmonic, and account for 14 different possible flux offsets.

Input data ‘x’ must therefore be of size (14xN); if not, it will be reshaped into that.

p[i,j] = (1. + param[5+j]) * [ param[0] - param[1]*cos(2*pi*x[i,j]) - param[2]*sin(2*pi*x[i,j]) -

param[3]*cos(4*pi*x[i,j]) - param[4]*sin(4*pi*x[i,j]) ]

[note that the first c-parameter (individual offset) will be

constrained such that: prod(1. + param[5::]) = 1.]

phasecurves.prepEclipseMap(nslice, npix0, k, b, z=None)

Set up variables necessary for (DISCRETE!) eclipse-mapping.

INPUTS:

nslice : int

Number of slices in model

npix : int

Number of pixels across the digitized planet face.

k : scalar

Planet/star radius ratio.

b : scalar

Transit impact parameter [always < (1+k) ]

z : 1D Numpy array

Transit crossing parameter (i.e., separation between geometric planet and stellar centers) at epochs of interest.

RETURNS:

slicemasks : Numpy array

[N x M x M] boolean array, one for each map cell.

timemasks : Numpy array

[N x T] boolean array. Equals true for the data indices which contribute to the corresponding slice mask.

cumulative_ingress_masks, cumulative_egress_masks : 3D Numpy arrays

[nslice x M x M] boolean maps of planet occultation at various stages, from prepEclipseMap()

ingress_zmasks, egress_zmasks : 2D Numpy arrays

[nslice x N] boolean maps of which slices correspond to which z-indices, from prepEclipseMap()

phasecurves.prolatesize(f, phase)

The projected size of a prolate ellipsoidal planet, viewed edge-on, at a given orbital phase.

INPUTS:

f : scalar

Ratio of large and small axes of the ellipsoid

phase : scalar or NumPy array

Orbital phase, in radians. 0 = transit, 1.57 ~ quadrature, etc.

OUTPUTS:

a_scale

Value by which to scale up the area of a circle (f=1).

REFERENCE:

Vickers 1996: http://dx.doi.org/10.1016/0032-5910(95)03049-2

phasecurves.ramp10(params, phase, args={'guess': [1, 0.2], 'n': 2})

Model Ramp Eq. 10 from Stevenson et al. (2011).

params: 2-sequence

parameters that define the function, as shown below.

phase: NumPy array.

Orbital phase (or more generally, ‘time’)

Functional form:

ramp = r[0] * (1. + r[1] * (phase - 0.5))

phasecurves.ramp11(params, phase, args={'guess': [1, 0.14, -1.9], 'n': 3})

Model Ramp Eq. 11 from Stevenson et al. (2011).

params: 3-sequence

parameters that define the function, as shown below.

phase: NumPy array.

Orbital phase (or more generally, ‘time’)

Functional form:

ramp = r[0] * (1. + r[1] * (phase - 0.5) + r[2] * (phase - 0.5)**2)

phasecurves.ramp2n(params, phase, args={'guess': [1, 26.6, 7.8], 'n': 3})

Model Ramp Eq. 2 (negative) from Stevenson et al. (2011).

params: 3-sequence

parameters that define the function, as shown below.

phase: NumPy array.

Orbital phase (or more generally, ‘time’)

Functional form:

ramp = r[0] * (1. - np.exp(-r[1]*phase + r[2]))

phasecurves.ramp2p(params, phase, args={'guess': [1, -0.16, 4.2], 'n': 3})

Model Ramp Eq. 2 (positive) from Stevenson et al. (2011).

params: 3-sequence

parameters that define the function, as shown below.

phase: NumPy array.

Orbital phase (or more generally, ‘time’)

Functional form:

ramp = r[0] * (1. + np.exp(-r[1]*phase + r[2]))

phasecurves.ramp3n(params, phase, args={'guess': [1, 141, 57.7, 0.123], 'n': 4})

Model Ramp Eq. 3 (negative) from Stevenson et al. (2011).

params: 4-sequence

parameters that define the function, as shown below.

phase: NumPy array.

Orbital phase (or more generally, ‘time’)

Functional form:

ramp = r[0] * (1. - np.exp(-r[1]*t + r[2]) + r[3] * (t - 0.5))

phasecurves.ramp3p(params, phase, args={'guess': [1, -0.16, 4.2, 0.1], 'n': 4})

Model Ramp Eq. 3 (positive) from Stevenson et al. (2011).

params: 4-sequence

parameters that define the function, as shown below.

phase: NumPy array.

Orbital phase (or more generally, ‘time’)

Functional form:

ramp = r[0] * (1. + np.exp(-r[1]*t + r[2]) + r[3] * (t - 0.5))

phasecurves.ramp4n(params, phase, args={'guess': [1, -0.00037, -0.94, 0.087, -1.08], 'n': 5})

Model Ramp Eq. 4 (negative) from Stevenson et al. (2011).

params: 5-sequence

parameters that define the function, as shown below.

phase: NumPy array.

Orbital phase (or more generally, ‘time’)

Functional form:

ramp = r[0] * (1. - np.exp(-r[1]*phase + r[2]) + r[3] * (phase - 0.5) + r[4] * (phase - 0.5)**2)

phasecurves.ramp4p(params, phase, args={'guess': [1, -0.068, 2.33, 0.933, -20.5], 'n': 5})

Model Ramp Eq. 4 (positive) from Stevenson et al. (2011).

params: 5-sequence

parameters that define the function, as shown below.

phase: NumPy array.

Orbital phase (or more generally, ‘time’)

Functional form:

ramp = r[0] * (1. + np.exp(-r[1]*phase + r[2]) + r[3] * (phase - 0.5) + r[4] * (phase - 0.5)**2)

phasecurves.ramp5n(params, phase, args={'guess': [1.0, 20, 83, 8.1, -0.1], 'n': 5})

Model Ramp Eq. 5 (negative) from Stevenson et al. (2011).

params: 5-sequence

parameters that define the function, as shown below.

phase: NumPy array.

Orbital phase (or more generally, ‘time’)

Functional form:

ramp = r[0] * (1. - np.exp(-r[1]*phase + r[2]) - np.exp(-r[3]*phase + r[4]))

phasecurves.ramp5p(params, phase, args={'guess': [1, -0.32, 2, -0.08, 2], 'n': 5})

Model Ramp Eq. 5 (positive) from Stevenson et al. (2011).

params: 5-sequence

parameters that define the function, as shown below.

phase: NumPy array.

Orbital phase (or more generally, ‘time’)

Functional form:

ramp = r[0] * (1. + np.exp(-r[1]*phase + r[2]) + np.exp(-r[3]*phase + r[4]))

phasecurves.ramp6(params, phase, args={'guess': [1, 0.053, 0.004, 0.4], 'n': 4})

Model Ramp Eq. 6 from Stevenson et al. (2011).

params: 4-sequence

parameters that define the function, as shown below.

phase: NumPy array.

Orbital phase (or more generally, ‘time’)

Functional form:

ramp = r[0] * (1. + r[1] * (phase - 0.5) + r[2] * np.log(phase - r[3]))

phasecurves.ramp7(params, phase, args={'guess': [1, 0.034, 0.35, 0.005, 0.35], 'n': 5})

Model Ramp Eq. 7 from Stevenson et al. (2011).

params: 5-sequence

parameters that define the function, as shown below.

phase: NumPy array.

Orbital phase (or more generally, ‘time’)

Functional form:

ramp = r[0] * (1. + r[1] * (phase - 0.5) + r[2] * (phase - 0.5)**2 + r[3] * np.log(phase - r[4]))

phasecurves.ramp8(params, phase, args={'guess': [1, 0.0096, 0.35, 0.00053], 'n': 4})

Model Ramp Eq. 8 from Stevenson et al. (2011).

params: 4-sequence

parameters that define the function, as shown below.

phase: NumPy array.

Orbital phase (or more generally, ‘time’)

Functional form:

ramp = r[0] * (1. + r[1] * np.log(phase - r[2]) + r[3] * np.log(phase - r[2])**2)

phasecurves.ramp9(params, phase, args={'guess': [1, 0.003, 0.6, 0.009, 0.35, 0.0004], 'n': 6})

Model Ramp Eq. 9 from Stevenson et al. (2011).

params: 6-sequence

parameters that define the function, as shown below.

phase: NumPy array.

Orbital phase (or more generally, ‘time’)

Functional form:

ramp = r[0] * (1. + r[1] * (phase - 0.5) + r[2] * (phase - 0.5)**2 + r[3] * np.log(phase - r[4]) + r[5] * np.log(phase - r[4])**2)

phasecurves.resfunc(*arg, **kw)

Generic function to give the error-weighted deviates on a function or functions:

INPUTS:

(fitparams, function, arg1, arg2, ... depvar, errs)

OR:

(allparams, (args1, args2, ..), npars=(npar1, npar2, ...))

where allparams is an array concatenation of each functions input parameters.

OPTIONAL INPUTS:  

jointpars – list of 2-tuples.

For use with multi-function calling (w/npars keyword). Setting jointpars=[(0,10), (0,20)] will always set params[10]=params[0] and params[20]=params[0].

EXAMPLE:

SEE ALSO:

resfunc()

phasecurves.rotmod(param, airmass, rotang, phase=None)

Model the Bean & Seifert rotation angle effect.

INPUTS:

param : 3- or 4-sequence

airmass : NumPy array

rotang : NumPy array

Instrument rotator angle, in radians.

phase : NumPy array, or None

Orbital phase (or time, etc.)

OUTPUTS:

param[0] + param[1] * airmass * np.cos(param[2] + rotang) + param[3] * (phase - phase.mean())

phasecurves.sin2_errs(params, eparams, nphi=10000.0, ntrials=10000.0)

Estimate the uncertainties from a double-sinusoid fit.

FUNCTION:

p[0] - p[1]*cos(phi + p[2]) + p[3]*cos(2*phi + p[4])

INPUTS:

params : 5-sequence

parameters for the function, as defined immediately above.

eparams : 5-sequence

1-sigma uncertainties on parameters

OPTIONS:

ntrials : float

number of Monte Carlo trials to run

nphi : float

number of points in phase curve (0-1, inclusive)

RETURNS:

(visibilities, peak_offset (rad), trough_offset (rad), true_vals)

SEE_ALSO:

phasesinsin14()

phasecurves.singleexp(params, t)

Model a simple, single-exponential function.

params: 2- or 3-sequence, defining the function as shown below.

t: sequence. Input time (presumed: since start of observations)

Functional form:

if len(params)==2:

F’/F = 1 - p0*exp(-t/p1)

elif len(params)>2:

F’/F = p2 * (1 - p0*exp(-t/p1))

phasecurves.singleexp14(params, t)

Model a simple, single-exponential function.

params: 16- or 17-sequence, defining the function as shown below.

t: sequence. Input time (presumed: since start of observations)

Functional form:

if len(params)==2:

F’/F = 1 - p0*exp(-t/p1)

elif len(params)>2:

F’/F = p2 * (1 - p0*exp(-t/p1))

... with fourteen additional sensitivity parameters.

phasecurves.slicemodel(param, xi)

Compute a slice model via Cowan & Agol (2008).

xi is from 0 to 2*pi

phasecurves.subfit_kw(params, input, i0, i1)

Parse keyword inputs (for errfunc(), etc.) for multiple concatenated inputs (e.g., with input key ‘npars’ set).

INPUTS:

params : 1D NumPy array of parameters

input : dict of keywords

i0 : int, first index of current parameters

i1 : int, last index of current parameters (e.g.: params[i0:i1])

phasecurves.transit_single(param, t)

compute 6+L-parameter eclipse function of a single event

velocity, impact parameter, stellar flux, planet/star radius ratio, time of center transit, period, limb darkening

param : parameters to be fit:

[Fstar, t_center, b, v (in Rstar/day), p (Rp/Rs), per (days)]

Up to two additional parameter can be concatenated onto the end, respectively representing linear and quadratic limb-darkening.

Input data ‘t’ must be of size (14xN); if not, it will be reshaped

into that.

phasecurves.tsin(param, x)

Compute phase function with unknown period, assuming a sinusoid:

p(x) = param[0] - param[1] * cos(2*pi*x/param[2] + param[3])

phasecurves.unwrap_joint_params(params, jfw_indices)

Unwrap parameters that are jointly constrained.

INPUTS:

params – 1D NumPy array

The P non-redundant parameters to the input function –

i.e., any parameters which are to be jointly fit (i.e., both held to the same, floating, value) are included only once.

jfw_indices : sequence of scalars and sequences

A length-P sequence of scalars and sequences. Each element jfw_indices[i] indicates the indices in the unwrapped set of parameters that will be assigned the value of params[i].

The final value of jfw_indices should be an integer equal to the length of the final set of unwrapped parameters.

EXAMPLE:

import tools
import numpy as np

npts = 100
snr = 20.
params = [1, 0.5, 1, 1]
x = np.linspace(0, 1, npts)
y = np.polyval(params, x) + np.random.randn(npts)/snr
jointpars = [(0, 2), (0, 3)]

joint_guess = np.array([1, 0.5])
jfw_indices = [[0, 2, 3], [1], 4]
full_params = tools.unwrap_joint_params(joint_guess, jfw_indices)

SEE_ALSO:

wrap_joint_params()

phasecurves.visit_offsets(visitcoef, masks)

INPUTS:

visitcoef : 1D NumPy array

offsets from unity for each of N HST visits

masks : 2D NumPy array, N x M

Set of boolean masks for each of N visits (not orbit!), assuming M total observations.

NOTES:Note that visitcoef[0] will be set so that the quantity (1. + visitcoef).prod() always equals unity.

phasecurves.wrap_joint_params(params, jointpars)

Wrap parameters that are jointly constrained.

INPUTS:

params – 1D NumPy array

All parameters, some of which may be jointly constrained.

jointpars – list of 2-tuples.

For use with multi-function calling (w/npars keyword). Setting jointpars=[(0,10), (0,20)] will always set params[10]=params[0] and params[20]=params[0].

EXAMPLE:

import tools
import numpy as np

npts = 100
snr = 20.
params = [1, 0.5, 1, 1]
x = np.linspace(0, 1, npts)
y = np.polyval(params, x) + np.random.randn(npts)/snr
jointpars = [(0, 2), (0, 3)]

all_params, joint_indices = tools.wrap_joint_params(full_params, jointpars)
wrapped_params = tools.unwrap_joint_params(all_params, joint_indices)

SEE_ALSO:

unwrap_joint_params()