Contents:
Functions for fitting phase curves.
Use of ‘errfunc’ is encouraged.
REQUIREMENTS: | Transit light curve routines |
---|
Generic function to give the weighted residuals on a function or functions:
INPUTS: |
OR:
OR:
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: | |
|
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))
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.
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)
Model a different double-exponential ramp function.
INPUTS: |
|
---|---|
Functional_form: | |
|
Model a different double-exponential ramp function w/14 positional offsets.
INPUTS: |
|
---|---|
Functional_form: | |
return (F * (1 + p[3/4::]).reshape(14,1) ) |
compute 3-parameter eclipse function of a single event, 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] = (...) * (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. ]
compute 3-parameter eclipse function of a single event
INPUTS: |
OR:
OR:
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: | |
|
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))
Generic function to give the chi-squared error on a generic function:
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.
Helper funtion for fitcurve – XXX
Param – either [phi_0] or [phi_0, inclination]
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
Compute labertian phase function with a fixed period=1, and x in units of orbital cycles.
param = [DC pedestal, AC amplitude, inclination (radians)]
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)]
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)]
Return a lambertian phase function with peak-to-valley amplitude unity.
Return amplitude of a nominally unity-amplitude lambertian with a given inclination angle, using function ‘lambertian’. inc is in radians.
Return mean of a nominally unity-amplitude lambertian with a given inclination angle, using function ‘lambertian’. inc is in radians.
Return natural logarithm of posterior probability: i.e., -chisq/2.
Inputs are the same as for errfunc().
SEE ALSO: | gaussianprocess.negLogLikelihood() |
---|
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.
Take a 2D planet map, and convert it into an eclipse light curve.
INPUTS: |
|
---|
Applies Markov Chain Monte Carlo model fitting using the Metropolis-Hastings algorithm.
INPUTS: |
|
---|---|
Returns: |
|
REFERENCES: | Numerical Recipes, 3rd Edition (Section 15.8); Wikipedia |
REVISIONS: |
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 |
MCMC for 17-parameter eclipse function of a single event
INPUTS: |
|
---|---|
RETURNS: |
|
REFERENCES: | Numerical Recipes, 3rd Edition (Section 15.8); Wikipedia |
MCMC for 3-parameter eclipse function of a single event
INPUTS: |
|
---|---|
RETURNS: |
|
REFERENCE: | Numerical Recipes, 3rd Edition (Section 15.8); Wikipedia |
Allow modeling with some parameters held constant.
INPUTS: |
fixedindex : sequence of ints
Thus if param = [10, 20, 50, 70] and holdfixed = [2], one would set varparam = [10, 50, 70] and fixedparam = [20]. |
---|---|
OPTIONS: | |
OUTPUTS: |
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]) - ...
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]
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]))
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])
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)
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.
compute phase function with a fixed period=1, assuming a sinusoid:
p(x) = param[0] - param[1] * cos(2*pi*x + param[2])
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]))
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])
[note that the the individual offsets will be subjected to the constraint: param[3::] -= param[3::].mean() ]
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
[note that the the individual offsets will be subjected to the constraint: param[3] = 1./(1.+param[4:17]).prod() - 1. ]
compute phase function with a fixed period=1 and phase offset (in radians):
p(x) = param[0] - param[1] * cos(2*pi*x + phaseoffset)
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: | The magnitude of the amplitudes will always be taken; they cannot be negative. |
---|
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.
Set up variables necessary for (DISCRETE!) eclipse-mapping.
INPUTS: |
|
---|---|
RETURNS: |
|
The projected size of a prolate ellipsoidal planet, viewed edge-on, at a given orbital phase.
INPUTS: |
|
---|
Model Ramp Eq. 10 from Stevenson et al. (2011).
Model Ramp Eq. 11 from Stevenson et al. (2011).
Model Ramp Eq. 2 (negative) from Stevenson et al. (2011).
Model Ramp Eq. 2 (positive) from Stevenson et al. (2011).
Model Ramp Eq. 3 (negative) from Stevenson et al. (2011).
Model Ramp Eq. 3 (positive) from Stevenson et al. (2011).
Model Ramp Eq. 4 (negative) from Stevenson et al. (2011).
Model Ramp Eq. 4 (positive) from Stevenson et al. (2011).
Model Ramp Eq. 5 (negative) from Stevenson et al. (2011).
Model Ramp Eq. 5 (positive) from Stevenson et al. (2011).
Model Ramp Eq. 6 from Stevenson et al. (2011).
Model Ramp Eq. 7 from Stevenson et al. (2011).
Model Ramp Eq. 8 from Stevenson et al. (2011).
Model Ramp Eq. 9 from Stevenson et al. (2011).
Generic function to give the error-weighted deviates on a function or functions:
INPUTS: |
OR:
|
---|---|
OPTIONAL INPUTS: | |
|
EXAMPLE:
Model the Bean & Seifert rotation angle effect.
INPUTS: | param : 3- or 4-sequence airmass : NumPy array
|
---|---|
OUTPUTS: | param[0] + param[1] * airmass * np.cos(param[2] + rotang) + param[3] * (phase - phase.mean()) |
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: |
|
OPTIONS: |
|
RETURNS: | (visibilities, peak_offset (rad), trough_offset (rad), true_vals) |
SEE_ALSO: |
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)
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)
F’/F = p2 * (1 - p0*exp(-t/p1))
... with fourteen additional sensitivity parameters.
Compute a slice model via Cowan & Agol (2008).
xi is from 0 to 2*pi
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]) |
---|
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.
Compute phase function with unknown period, assuming a sinusoid:
p(x) = param[0] - param[1] * cos(2*pi*x/param[2] + param[3])
Unwrap parameters that are jointly constrained.
INPUTS: |
|
---|---|
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: |
INPUTS: |
|
---|
Wrap parameters that are jointly constrained.
INPUTS: |
|
---|---|
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: |