This page was created by the IDL library routine
mk_html_help. For more information on
this routine, refer to the IDL Online Help Navigator
or type:
? mk_html_help
at the IDL command line prompt.
Last modified: Tue Sep 16 17:50:51 2014.
NAME:
BJD2TARGET
PURPOSE:
Iteratively calls TARGET2BJD to convert a BJD in Barycentric
Dynamical Time (BJD_TDB) to a BJD in the target barycenter
time (BJD_TARGET) within TOL days.
DESCRIPTION:
The opposite of TARGET2BJD; see description there.
INPUTS:
BJD_TDB - A scalar or array of BJDs in TDB. Must be double
precision.
INCLINATION - The inclination of the orbit
A - The semi-major axis of the orbit (AU)
TP - The time of periastron of the orbit (BJD_TARGET)
PERIOD - The period of the orbit (days)
E - Eccentricity of the orbit
OMEGA - Argument of Periastron of the orbit (radians)
OPTIONAL INPUTS:
Q - The mass ratio of the targets (M1/M2). If not
specified, an infinite mass ratio is assumed (M1 is
stationary at the barycenter) (8 ms effect for Hot
Jupiters).
TOL - The tolerance, in days. The iterative procedure will
stop after the worst-case agreement is better than this.
Default = 1d-8 (1 ms).
OPTIONAL KEYWORDS:
PRIMARY - If set, the information comes from the position of
the primary (as in RV), and therefore the correction
will be the light travel time from the center of the
primary to the Barycenter -- analagous to the
difference between HJD and BJD in our solar system
(only ~8 ms for Hot Jupiters, but increasing with a).
Otherwise, the correction will be the light
travel time from the smaller body to the barycenter
(as in transits) -- analagous to the difference
between JD and BJD in the Solar System.
NOTE: if Q is not specified and PRIMARY is, no
correction is applied.
OUTPUTS:
BJD_TARGET - The time as it would flow in the Barycenter of the target.
LIMITATIONS:
We ignore the distance to the object (plane parallel waves), which
should have a similar effect as the distance plays in the BJD
correction (< 1 ms). We also ignore the systemic velocity, which
will compress/expand the period by a factor gamma/c.
REVISION HISTORY:
2011/06: Written by Jason Eastman (OSU)
(See bjd2target.pro)
NAME:
ELLKE
PURPOSE:
Computes Hasting's polynomial approximation for the complete
elliptic integral of the first (ek) and second (kk) kind. Combines
the calculation of both so as not to duplicate the expensive
calculation of alog10(1-k^2).
CALLING SEQUENCE:
ellke, k, ek, kk
INPUTS:
k - The elliptic modulus.
OUTPUTS:
ek - The elliptic integral of the first kind
kk - The elliptic integral of the second kind
MODIFICATION HISTORY
2009/04/06 -- Written by Jason Eastman (Ohio State University)
(See ellke.pro)
NAME:
ELLPIC_BULIRSCH
PURPOSE:
Computes the complete elliptical integral of the third kind using
the algorithm of Bulirsch (1965):
Bulirsch 1965, Numerische Mathematik, 7, 78
Bulirsch 1965, Numerische Mathematik, 7, 353
CALLING SEQUENCE:
result = ellpic_bulirsch(n, k)
INPUTS:
n,k - int(dtheta/((1-n*sin(theta)^2)*sqrt(1-k^2*sin(theta)^2)),0, pi/2)
RESULT:
The complete elliptical integral of the third kind
MODIFICATION HISTORY
2009/03 -- Written by Eric Agol
(See ellpic_bulirsch.pro)
NAME:
EXOFAST
PURPOSE:
Simultaneously fits RV and/or transit data for a single
planet. Please cite Eastman et al., 2013
(http://adsabs.harvard.edu/abs/2013PASP..125...83E) if you make
use of this routine in your research. Please report errors or bugs
to jeastman@lcogt.net
CALLING SEQUENCE:
exofast, [RVPATH=, TRANPATH=, BAND=, PRIORS=, PREFIX=, /CIRCULAR,
/NOSLOPE, /SECONDARY, /UPDATE, PNAME=, SIGCLIP=, NTHIN=,
MAXSTEPS=, MINPERIOD=, MAXPERIOD=, NMIN=, /DISPLAY,
/DEBUG, RANDOMFUNC=, SEED=, /SPECPRIORS, /BESTONLY,
NINTERP=, EXPTIME=, /LONGCADENCE]
INPUTS:
RVPATH - The path to the RV data file. The file must have 3 columns:
1) Time (BJD_TDB -- See Eastman et al., 2010)
2) RV (m/s)
3) err (m/s)
NOTE: The units must be as specified, or the fit
will be wrong or fail.
NOTE 2: If omitted, just the transit data will be fit
TRANPATH - The path to the transit data file. The file must
have at least 3 columns:
1) Time (BJD_TDB -- See Eastman et al., 2010)
2) Normalized flux
3) err
4) Detrend parameter 1
....
N+3) Detrend parameter N
NOTE: The units must be as specified, or the fit
will be wrong or fail.
NOTE 2: If omitted, just the RV data will bit fit
BAND - The bandpass of the observed transit (see quadld.pro
for allowed values).
NOTE: only required if TRANPATH is specified.
PRIORS - Priors on each of the parameters. See EXOFAST_CHI2
for parameter definitions. Must be an N x 2
elements array, where N is 15 + the number of
detrending variables. If set, all non-zero values of
priors[0,*] are used to start the fit, and a penalty
equal to:
total(((pars-priors[0,*])/priors[1,*])^2)
is added to the total chi^2.
NOTE 1: For no prior, set the width to infinity. i.e.,
priors[1,*] = !values.d_infinity.
NOTE 2: Typically, the default starting guesses are
good enough and need not be changed.
TRANSIT+RV or TRANSIT-ONLY FITS: Priors on Teff and
[Fe/H] (priors[*,11:12]) must be specified, either by
declaring the priors array or specifying a planet
name and setting the /SPECPRIORS keyword.
RV-ONLY FITS: Priors on logg, Teff and [Fe/H]
(priors[*,10:12]) must be specified, either by
declaring the priors array or specifying a planet
name and setting the /SPECPRIORS keyword.
OPTIONAL INPUTS:
PREFIX - Each of the output files will have this string as a
prefix. Default is RVFILE without the
extension. Cannot contain an underscore, "_".
MINPERIOD - The minimum period to consider. The default is 1 day.
MAXPERIOD - The maximum period to consider. The default is the
range of the RV input times.
NMIN - The number of minima in the Lomb-Scargle Periodogram
to fit a full Keplerian orbit. Default is 5. If the
eccentricity is large, this may need to be
increased. The execution time of the initial global
fit is directly proportional to this value (though
is a small fraction of the total time of the MCMC fit).
PNAME - If set, starting values from exoplanets.org for the
planet PNAME will be used to seed the fit. It is
insensitive to extra spaces and capitalization.
e.g., "WASP-12 b" is the same as "wasp-12b".
SIGCLIP - If set, an iterative fit will be performed excluding
data points more than SIGCLIP*error from the best
fit. Default is no clipping.
NTHIN - If set, only every NTHINth element will be
kept. This typically doesn't affect the
resultant fit because there is a high correlation
between adjacent steps and has the advantage of
improved memory management and faster generation of
the final plots.
MAXSTEPS - The maximum number of steps to take in the MCMC
chain. Note that a 32-bit installation of IDL
cannot allocate more than 260 million
double-precision numbers, and redundant copies of
each parameter are required. A very large number
will cause memory management problems. Default is
100,000.
RANDOMFUNC - A string specifying the name of the random number
generator to use. This generator must be able to
return 1,2 or 3 dimensional uniform or normal random
deviates. Default is 'EXOFAST_RANDOM' (which is slow
but robust).
SEED - The seed to the random number generator used for the
MCMC fits. Be sure not to mix seeds from different
generators. The default is -systime(/seconds).
NINTERP - If set, the each model data point will be an average
of NINTERP samples over a duration of EXPTIME
centered on the input time.
EXPTIME - If set, the each model data point will be an average
of NINTERP samples over a duration of EXPTIME
centered on the input time.
For a dicussion on binning, see:
http://adsabs.harvard.edu/abs/2010MNRAS.408.1758K
OPTIONAL KEYWORDS:
CIRCULAR - If set, the fit will be forced to be circular (e=0,
omega_star=pi/2)
NOSLOPE - If set, it will not fit a slope to the RV data.
SECONDARY - If set, fit a secondary eclipse. This feature is not
well-tested -- use at your own risk.
UPDATE - Update the local copy of the exoplanets.org file (only
applied if PNAME is specified).
DISPLAY - If set, the plots (below) will be displayed via
ghostview (gv).
SPECPRIORS- If set (and PNAME is specified), spectroscopic priors
from exoplanets.org will be used for logg, Teff, and
[Fe/H].
BESTONLY - If set, only the best fit (using AMOEBA) will be
performed.
DEBUG - If set, various debugging outputs will be printed and
plotted.
LONGCADENCE - If set, EXPTIME=29.425 and NINTERP=10 are set to handle
long cadence data from Kepler.
OUTPUTS:
Each of the output files will be preceeded by PREFIX (defined
above). The first "?" will be either "c" for circular if /circular
is set or "e" for eccentric. The second "?" will be either "f" for
flat (if /noslope is set) or "m "for slope.
mcmc.?.?.idl - An IDL save file that contains the full chains
for each parameter, including derived paramters,
the corresponding names, the chi2 at each link,
and the index of the burn-in period.
pdfs.?.?.ps - A postscript plot of each posterior distribution
function, 8 to a page.
covars.?.?.ps - A postscript plot of the covariances between each
parameter, 16 to a page. The title of each plot is the
correlation coefficient.
median.?.?.tex - The LaTeX source code for a deluxe table of the median
values and 68% confidence interval, rounded appropriately.
best.?.?.tex - The LaTeX source code for a deluxe table of the best fit
values and 68% confidence interval, rounded appropriately.
model.?.?.ps - An postscript plot of the best-fit models and residuals.
model.?.?.rv - A text file containing the best-fit model RV for each
input time
model.?.?.flux - A text file containing the best-fit model flux for each
input time
NOTE: To extract a single page out of a multi-page PS file, use:
psselect -p# input.ps output.ps
where # is the page number to extract.
COMMON BLOCKS:
CHI2_BLOCK:
RV - A structure with three tags, corresponding to each
column of RVFILE.
TRANSIT - A structure with 3 + N tags, corresponding to each
column of the TRANFILE.
PRIORS - The priors on each of the parameters
BAND - The band of the observed transit
DEBUG - A flag to specify if the data/fit should be
plotted. This is intended for debugging only; it
dramatically increases the execution time.
NORV - A boolean to indicate the RV should not be fit
NOTRAN - A boolean to indicate the transit should not be fit
RV_BLOCK:
RVDATA - A structure with three tags, corresponding to each
column of RVFILE. This is to prefit the data with a
lomb-scargle periodogram
EXAMPLES:
;; fit HAT-P-3b example data (using priors from exoplanets.org):
exofast, rvpath='hat3.rv',tranpath='hat3.flux',pname='HAT-P-3b',$
band='Sloani',/circular,/noslope,/specpriors,minp=2.85,maxp=2.95
;; fit HAT-P-3b without using exoplanets.org
;; include a prior on the period too
;; this reproduces results from Eastman et al., 2013
per = 2.899703d0
uper = 5.4d-5
priors = dblarr(2,19)
priors[1,*] = !values.d_infinity
priors[*, 3] = [alog10(per),uper/(alog(10d0)*per)] ;; logp prior
priors[*,10] = [4.61d0,0.05d0] ;; logg prior
priors[*,11] = [5185d0, 80d0] ;; Teff prior
priors[*,12] = [0.27d0,0.08d0] ;; logg prior
exofast,rvpath='hat3.rv',tranpath='hat3.flux',band='Sloani',/circular,$
/noslope,priors=priors,minp=2.85,maxp=2.95,prefix='example.'
MODIFICATION HISTORY
2012/06 -- Public Release, Jason Eastman (LCOGT)
2012/12 -- Upgrades for Kepler data, JDE, LCOGT:
Change treatment of limb darkening
Now fits limb darkening using Claret tables as a prior
Add LongCadence flag, exptime, niterp keywords
2012/12 -- Update documentation with corrected priors array
2012/12 -- Now cosi scale (for AMOEBA) is rstar/a instead of hard
coded at 0.1. More robust for long-period planets.
Added additional output files
2012/12 -- Added Carter error estimates for BESTONLY transit fits
Fit eccentricity for transit-only fits
requires logg prior
ecosw fixed at zero for "best" fit
ecosw is only constrained by |ecosw| <= sqrt(1-esinw^2)
takes a long time for MCMC due to ecosw degeneracy
2013/01 -- Fix bugs with prior array misuse (works with example
in README again)
2013/02 -- Fixed bug in primary transit probabilities (extra factor
of Rstar/Rsun)
2013/06 -- Fixed bug when using /specpriors keyword (which broke example)
(See exofast.pro)
NAME:
EXOFAST_AMOEBA
PURPOSE:
Multidimensional minimization of a function FUNC(X), where
X is an N-dimensional vector, using the downhill simplex
method of Nelder and Mead, 1965, Computer Journal, Vol 7, pp 308-313.
This routine is based on the AMOEBA routine, Numerical
Recipes in C: The Art of Scientific Computing (Second Edition), Page
411, and is used by permission.
This is a debugged version of IDL's built-in AMOEBA, changing:
fac1 = (1.0 - fac) / n_elements(psum)
to
fac1 = (1.d0 - fac) / n_elements(psum)
CATEGORY:
Function minimization/maximization. Simplex method.
CALLING SEQUENCE:
Result = EXOFAST_AMOEBA(Ftol, ....)
INPUTS:
FTOL: the fractional tolerance to be achieved in the function
value. e.g. the fractional decrease in the function value in the
terminating step. This should never be less than the
machine's single or double precision.
KEYWORD PARAMETERS:
FUNCTION_NAME: a string containing the name of the function to
be minimized. If omitted, the function FUNC is minimized.
This function must accept an Ndim vector as its only parameter and
return a scalar single or double precision floating point value as its
result.
FUNCTION_VALUE: (output) on exit, an Ndim+1 element vector
containing the function values at the simplex points. The first
element contains the function minimum.
NCALLS: (output) the of times the function was evaluated.
NMAX: the maximum number of function evaluations allowed
before terminating. Default = 5000.
P0: Initial starting point, an Ndim element vector. The starting
point must be specified using either the keyword SIMPLEX, or P0 and
SCALE. P0 may be either single or double precision floating.
For example, in a 3-dimensional problem, if the initial guess
is the point [0,0,0], and it is known that the function's
minimum value occurs in the interval: -10 <
X(0) < 10, -100 < X(1) < 100, -200 < X(2) < 200, specify: P0=[0,0,0],
SCALE=[10, 100, 200].
SCALE: a scalar or Ndim element vector contaiing the problem's
characteristic length scale for each dimension.
SCALE is used with P0 to form an initial (Ndim+1) point simplex.
If all dimensions have the same scale, specify a scalar.
SIMPLEX: (output and/or optional input) On input, if P0 and SCALE
are not set, SIMPLEX contains the Ndim+1 vertices, each of
Ndim elements, of starting simplex, in either single or double
precision floating point, in an (Ndim, Ndim+1) array. On output,
SIMPLEX contains the simplex, of dimensions (Ndim, Ndim+1), enclosing
the function minimum. The first point, Simplex(*,0), corresponds to
the function's minimum.
OUTPUTS:
Result: If the minimum is found, an Ndim vector, corresponding to
the Function's minimum value is returned. If a function minimum
within the given tolerance, is NOT found in the given number of
evaluations, a scalar value of -1 is returned.
COMMON BLOCKS:
None.
SIDE EFFECTS:
None.
PROCEDURE:
This procedure implements the Simplex method, described in
Numerical Recipes, Section 10.4. See also the POWELL procedure.
Advantages: requires only function evaluations, not
derivatives, may be more reliable than the POWELL method.
Disadvantages: not as efficient as Powell's method, and usually
requires more function evaluations.
Results are performed in the mode (single or double precision)
returned by the user-supplied function. The mode of the inputs P0,
SCALE, or SIMPLEX, should match that returned by the function. The
mode of the input vector supplied to the user-written function, is
determined by P0, SCALE, or SIMPLEX.
EXAMPLE:
Use Amoeba to find the slope and intercept of a straight line fitting
a given set of points minimizing the maximum error:
The function to be minimized returns the maximum error,
given p(0) = intercept, and p(1) = slope:
FUNCTION FUNC, p
COMMON FUNC_XY, x, y
RETURN, MAX(ABS(y - (p(0) + p(1) * x)))
END
Put the data points into a common block so they are accessible to the
function:
COMMON FUNC_XY, x, y
Define the data points:
x = findgen(17)*5
y = [ 12.0, 24.3, 39.6, 51.0, 66.5, 78.4, 92.7, 107.8, 120.0, $
135.5, 147.5, 161.0, 175.4, 187.4, 202.5, 215.4, 229.9]
Call the function. Fractional tolerance = 1 part in 10^5,
Initial guess = [0,0], and the minimum should be found within
a distance of 100 of that point:
r = EXOFAST_AMOEBA(1.0e-5, SCALE=1.0e2, P0 = [0, 0], FUNCTION_VALUE=fval)
Check for convergence:
if n_elements(r) eq 1 then message,'AMOEBA failed to converge'
Print results.
print, 'Intercept, Slope:', r, 'Function value (max error): ', fval(0)
Intercept, Slope: 11.4100 2.72800
Function value: 1.33000
MODIFICATION HISTORY:
DMS, May, 1996. Written.
May 2010 - Bug fix for double precision steps. Jason Eastman, OSU
(See exofast_amoeba.pro)
NAME:
EXOFAST_CHI2
PURPOSE:
Computes the chi2 for a transit and/or RV for a single planet
CALLING SEQUENCE:
chi2 = exofast_chi2(pars)
INPUTS:
PARS - a parameter array containing all of the parameters in the
model.
gamma = pars[0] ;; systemic velocity
slope = pars[1] ;; slope in RV
tc = pars[2] ;; transit center time
logP = pars[3] ;; alog10(Period/days)
sqrtecosw = pars[4] ;; eccentricity/arg of periastron
sqrtesinw = pars[5] ;; eccentricity/arg of periastron
logK = pars[6] ;; alog10(velocity semi-amplitude/(m/s))
cosi = pars[7] ;; cosine of inclination of the orbit
p = pars[8] ;; rp/rstar
log(ar) = pars[9] ;; alog10(a/rstar)
logg = pars[10] ;; stellar surface gravity
teff = pars[11] ;; stellar effective temperature
feh = pars[12] ;; stellar metallicity
depth2 = pars[13] ;; secondary eclipse depth
u1 = pars[14] ;; linear limb darkening coeff
u2 = pars[15] ;; quadratic limb darkening coeff
u3 = pars[16] ;; 1st non-linear limb darkening coeff (not supported)
u4 = pars[17] ;; 2nd non-linear limb darkening coeff (not supported)
F0 = pars[18] ;; baseline flux
coeffs = pars[19:npars-1] ;; detrending variables
OPTIONAL INPUTS:
PSNAME - The name of a PS file. If set, a plot the
data/model will be written to this file.
OPTIONAL OUTPUTS:
DETERMINANT - The determinant of the parameterization above and
the uniform priors we wish to impose. In this case,
it is always 1d0 (but is required by EXOFAST_DEMC).
MODELRV - The RV model at each time (rv.bjd).
MODELFLUX - The model light curve at each time (transit.bjd).
RESULT:
The chi^2 of the model given the data and parameters.
COMMON BLOCKS:
CHI2_BLOCK - See exofast.pro for definition
MODIFICATION HISTORY
2012/06 -- Public release -- Jason Eastman (LCOGT)
2012/07 -- Fixed major bug in mstar/rstar prior width derivation
2012/12 -- Add Long cadence, quadratic limb darkening fit.
2012/12 -- Changed eccentricity constraint to e < (1-Rstar/a)
2013/02 -- Fixed bug that broke detrending, introduced in 2012/12
(See exofast_chi2.pro)
NAME:
exofast_demc
PURPOSE:
Make NCHAINS independent Markov Chain Monte Carlo chains to fit data
and return the parameter distributions.
DESCRIPTION:
Begins by determining the correct stepping scale by finding
differences in each parameters that yield a delta chi^2 = 1, then
uses 5*scale*randomn offsets in each parameter for the starting
points of each chain.
It then begins a Differential Evolution Markov Chain Monte Carlo
fit (ter Braak, 2006)
http://www.stat.columbia.edu/~gelman/stuff_for_blog/cajo.pdf. The
only slight modification to this basic alogorithm is that we
determine the magnitude of the uniform random deviate for each
parameter. We have found the dynamic range of our parameters can
be large, and thus not well suited to the one-size fits all
approach recommended there.
After taking 5% of MAXSTEPS, the program will estimate how many steps
will be necessary to be well-mixed. If it is not expected to be
well-mixed before taking the maximum number of steps, it will
output a warning and a recommended setting for NTHIN. This should
be accurate at the factor of 2-3 level.
The program stops when the chains are well-mixed, as defined by
Ford 2006 (http://adsabs.harvard.edu/abs/2006ApJ...642..505F)
using the Gelman-Rubin Statistic and number of independent draws,
or when each chain has taken MAXSTEPS, whichever is first.
Every step in the chain before all chains have crossed below the
median chi^2 will be considered the (the "burn-in"), and not used
for convergence tests. BURNNDX will specify the index of the first
useable chain so as not to be biased by the starting criteria.
CALLING SEQUENCE:
exofast_demc, bestpars, 'mychi2', pars [,CHI2=, TOFIT=,$
SCALE=, SEED=,RANDOMFUNC=,NTHIN=, MAXSTEPS=,/DONTSTOP,$
NCHAINS=, ANGULAR=,BURNNDX=,/REMOVEBURN]
INPUTS:
BESTPARS - Array of the best-fit parameters. An accurate initial
fit is required to find the correct step size. MCMC
is not ideal for fitting data, just characterizing
the errors.
CHI2FUNC - A string that specifies the name of the user-defined
function that calculates the chi^2 of a given
parameter set. Required data should be put in a
COMMON block. Its calling sequence must be:
chi2 = chi2func(pars, determinant=determinant)
OPTIONAL INPUTS:
TOFIT - An array that indexes the parameters to fit. If not
specified, all parameters are fit.
ANGULAR - An array that indexes the parameters that are
angular (must be in radians). This will enable
special ways to calculate the median and convergence
statistics, which may otherwise fail. Default is none.
SCALE - An NPARS array containing the stepping scales for
each parameter. This is not recommended for normal
use. If not specified, it will be automatically
determined using EXOFAST_GETMCMCSCALE.
NTHIN - Saves only every NTHINth link in each chain. This
is only recommended when the autocorrelation between
steps is large and memory management issues may
arise. Results will be more accurate the smaller this
is -- as long as the chain is well-mixed (pay
attention to the warnings).
MAXSTEPS - The maximum number of steps (after thinning) for each
chain to take. The default is 100,000. Take care when
increasing this number, as memory allocation issues
may arise (especially with 32 bit machines). If the
chain is not well-mixed, it is usually better to
increase NTHIN instead.
NCHAINS - The number of independent chains to run. The
execution time scales ~linearly with this number, and
must be at least 3 to calculated convergence
statistics. The default is 10.
SEED - A seed for the random number generator. Do not mix
seeds from different random number
generators, and do not use several different seeds
for the same generator. If you use a random number
generator elsewhere, you should call EXOFAST_DEMC with
SEED. Default is -systime(/seconds).
RANDOMFUNC - A string for the name of the random number generator
function. You can use your own as long as it returns
a uniform double precision deviate between 0 and 1,
accepts the /NORMAL keyword to generate random
gaussian deviates, and can return at least a two
dimensional array of deviates. You may also choose
between "EXOFAST_RANDOM" (default) or "RANDOMU".
- EXOFAST_RANDOM is an IDL wrapper for pg_ran. This is
a state of the art, bullet-proof generator with no
known statistical flaws, based on the third edition
of numerical recipes, pg 342. It is ~120x slower than
RANDOMU, though for typical DEMC fits, this accounts
for ~1% of the total runtime.
- RANDOMU is IDL's built-in, ran1-based generator with a
periodicity of ~10^9. ran1 is formally admonished by
its author, Numerical Recipes, but shown to give
identical results as the state of the art generator
See Eastman et. al., 2013 for more details.
http://adsabs.harvard.edu/abs/2013PASP..125...83E
OPTIONAL KEYWORDS:
DONTSTOP - If set, the program will not stop when it is
converged, but will instead take MAXSTEPS steps in
each chain.
REMOVEBURN - If set, the burn-in will be removed in the returned
parameter array.
OUTPUTS:
PARS - NFIT x NSTEPS x NCHAINS array of parameters.
NOTE: NSTEPS will be the lesser of when it is
converged or MAXSTEPS.
CHI2 - NSTEPS x NCHAINS array with the chi^2 for each set of
parameters
REVISION HISTORY:
2012/06 - Public Release - Jason Eastman (LCOGT)
2012/12 - When parameters aren't mixed, display the index
to the (more) static PARS array, not the index to the
dynamic TOFIT array.
- Removed memory-hungry vestigial code
(See exofast_demc.pro)
NAME:
exofast_errell
PURPOSE:
Gives the path of the error ellipses at constant probability given
two parameter distributions.
PROCEDURE
1. Creates a 2D histogram of the data with NBINS bins per
parameter.
2. Calculates the fraction of data points above each value in the
histogram.
3. Compares each calculated fraction with the probabilities,
and defines the closest as the levels to pass to CONTOUR.
4. Creates a contour plot of the 2D histogram
5. Scales the path to be in the same units as the input and
recreates the path, adding NaNs between each level for easy
plotting
CALLING SEQUENCE:
exofast_errell, xpar, ypar [,NXBIN=nxbin,NYBIN=nybin,PROBS=PROBS]$
[,XPATH=xpath][,YPATH=ypath][,/PLOT][,/OPLOT]
INPUTS:
XPAR - A vector of parameters to be plotted on the x axis.
YPAR - A vector of parameters to be plotted on the y axis.
OPTIONAL INPUTS:
NXBIN - The number of bins to take in each direction. Default
is 100. If you have only a small number of values for
xpar and ypar (10000), this should be smaller.
NYBIN - The number of bins to take in each direction. Default
is NXBIN.
XMIN/YMIN - The minimum x/y value to make the contours
XMAX/YMAX - The maximum x/y value to make the contours
PROBS - An array containing the probability contours to
plot. The default is
probs = erf([1.d0,2.d0]/sqrt(2.d0)), or
probs = [0.6826894921370859d0,$
0.9544997361036416d0]
Note: to have a meaningful contour enclosing 99.7% of
the points, you must have ~1 million points or make
NBINS smaller.
OPTIONAL KEYWORDS:
PLOT - Plot the contours on a new plot
OPLOT - Overplot the contours on an existing plot
OPTIONAL OUTPUTS:
XPATH/YPATH - The x/y values of the contours. For more
customizable plotting, return these and plot inside
your own program. All contours can be plotted simply
with IDL> plot, xpath, ypath
NOTE: XPATH/YPATH contain NaNs between levels to make the
plotting easier.
OUTSIDENDX - The indices of the parameters outside of the
outer-most contour (requires David Fanning's Inside
procedure)
EXAMPLE:
plot the 68% and 95% contours, and every point outside them
x = randomn(seed,100000)
y = randomn(seed,100000)
exofast_errell, x,y,out=out,probs=[0.68,0.95],xpath=xpath,ypath=ypath
plot, xpath, ypath, /iso
if out[0] ne -1 then oplot, x[out],y[out],psym=3
DEPENDENCIES:
INSIDE (Coyote library):
http://www.idlcoyote.com/programs/inside.pro
NOTE: there is non-vectorized version of his code online, which will
not work. Use the one at the link above.
REVISION HISTORY:
2010/03/01 - Written: Jason Eastman - The Ohio State University
2010/11/05 - Added OUTSIDENDX output (inspired by Ben Shappee, OSU)
Added, NXBIN, XMIN, XMAX, NYBIN, YMIN, YMAX inputs.
2011/03/25 - Better memory management for INSIDE -- split up
array to check (takes a little longer, but actually
works with large arrays!)
(See exofast_errell.pro)
NAME:
exofast_gelmanrubin
PURPOSE:
Calculates the Gelman-Rubin Statistic for several parallel Markov
chains, and determines whether or not the chains are well mixed
according to the Ford 2006 prescription:
http://adsabs.harvard.edu/abs/2006ApJ...642..505F
This is meant to be a dependency of EXOFAST_DEMC. An example of
its use can be seen there.
PROCEDURE
Calculates the Gelman-Rubin statistic and the number of
independent draws for each parameter, as defined by Ford, 2006. The
chain is considered well-mixed if all parameters have a
Gelman-Rubin statistic of <= 1.01 and >= 1000 independent draws.
CALLING SEQUENCE:
ismixed = exofast_gelmanrubin(pars)
INPUTS:
PARS0 - A 3 dimensional array (NPARS,NSTEPS,NCHAINS) of
parameter values
OPTIONAL OUTPUTS:
GELMANRUBIN - An NPARS element array containing the Gelman-Rubin
statistic for each parameter (equation 25).
TZ - An NPARS element array containing the number of
independent draws for each parameter (equation 26).
OPTIONAL KEYWORDS:
ANGULAR - An array of indices (with up to NPARS elements)
indicating which elements are angular values. If
specified, these parameters will be re-centered
about the mode value before calculating the
statisitics (the original array will not be
modified). If left blank, no scaling will be performed.
NOTE: angular values must be in radians.
REVISION HISTORY:
2010/03/01 - Written: Jason Eastman - The Ohio State University
(See exofast_gelmanrubin.pro)
NAME:
EXOFAST_GETB
PURPOSE:
This function returns the impact parameter as a function of BJD
given orbital elements of the planet, assuming a keplerian orbit.
Optionally, the 3 space coordinates can be returned.
CALLING SEQUENCE:
result = getb(jds, i=i,a=a,tperiastron=tperiastron,period=p,$
e=e,omega=omega, x=x, y=y, z=z)
INPUTS:
bjd - Barycentric Julians dates for desired impact
parameters (ideally in the target's
barycentric frame; see BJD2TARGET)
i - inclination of orbit (radians)
a - semi-major axis (in units of R_*)
tperiastron - periastron passage time (BJD)
Period - Period of orbit (days)
OPTIONAL INPUTS:
e - eccentricity of orbit (0 if not specified)
omega - argument of periastron of the star's orbit, in radians
-- omega_* is typically quoted from RV
-- required if e is specified
-- assumed to be pi/2 if e not specified
-- omega_* = omega_planet + !dpi
lonascnode - The Longitude of the ascending node
(radians). Assumed to be !dpi if not specified.
OUTPUTS:
result - the impact parameter as a function of BJD, in units
of R_*.
OPTIONAL OUTPUTS:
x,y,z - Arrays of the cartesian coordinates of the planet at
each BJD, in the units of R_*.
+X is right
+Y is up
+Z is out of the page (primary transit)
MODIFICATION HISTORY
2009/05/01 -- Jason Eastman (Ohio State University)
(See exofast_getb.pro)
NAME:
EXOFAST_GETCHI2_RV_FIT
PURPOSE:
Computes the chi^2 of a single planet decribed by PARS, while
analytically fitting for K, gamma, slope
CALLING SEQUENCE:
chi2 = exofast_getchi2_rv_fit(pars)
INPUTS:
The best-fit parameters for the RV fit of a single planet.
pars[0] = time of transit center
pars[1] = period
pars[2] = e*cos(omega)
pars[3] = e*sin(omega)
pars[4] = K (will be overwritten)
pars[5] = gamma (will be overwritten)
pars[6] = slope (will be overwritten)
RESULT:
The chi^2 of the parameters.
COMMON BLOCKS:
RV_BLOCK - See exofast.pro for definition
MODIFICATION HISTORY
2012/06 -- Public release -- Jason Eastman (LCOGT)
(See exofast_getchi2_rv_fit.pro)
NAME:
EXOFAST_GETCHI2_RV_FITCIR
PURPOSE:
Computes the chi^2 of a single, circular planet decribed by PARS,
while analytically fitting for Tc, K, gamma, slope
CALLING SEQUENCE:
chi2 = exofast_getchi2_rv_fitcir(pars)
INPUTS:
The best-fit parameters for the RV fit of a single planet.
pars[0] = time of transit center (will be overwritten)
pars[1] = period
pars[2] = e*cos(omega) (assumed 0)
pars[3] = e*sin(omega) (assumed 0)
pars[4] = K (will be overwritten)
pars[5] = gamma (will be overwritten)
pars[6] = slope (will be overwritten)
RESULT:
The chi^2 of the parameters.
COMMON BLOCKS:
RV_BLOCK - See exofast.pro for definition
MODIFICATION HISTORY
2012/06 -- Public release -- Jason Eastman (LCOGT)
(See exofast_getchi2_rv_fitcir.pro)
NAME:
EXOFAST_GETCHI2_RV_FITCIRNOM
PURPOSE:
Computes the chi^2 of a single, circular planet decribed by PARS, while
analytically fitting for Tc, K, gamma
CALLING SEQUENCE:
chi2 = exofast_getchi2_rv_fitcirnom(pars)
INPUTS:
The best-fit parameters for the RV fit of a single planet.
pars[0] = time of transit center (will be overwritten)
pars[1] = period
pars[2] = e*cos(omega) (assumed 0)
pars[3] = e*sin(omega) (assumed 0)
pars[4] = K (will be overwritten)
pars[5] = gamma (will be overwritten)
pars[6] = slope (assumed 0)
RESULT:
The chi^2 of the parameters.
COMMON BLOCKS:
RV_BLOCK - See exofast.pro for definition
MODIFICATION HISTORY
2012/06 -- Public release -- Jason Eastman (LCOGT)
(See exofast_getchi2_rv_fitcirnom.pro)
NAME:
EXOFAST_GETCHI2_RV_FITNOM
PURPOSE:
Computes the chi^2 of a single planet decribed by PARS, while
analytically fitting for K, gamma
CALLING SEQUENCE:
chi2 = exofast_getchi2_rv_fitnom(pars)
INPUTS:
The best-fit parameters for the RV fit of a single planet.
pars[0] = time of transit center
pars[1] = period
pars[2] = e*cos(omega)
pars[3] = e*sin(omega)
pars[4] = K (will be overwritten)
pars[5] = gamma (will be overwritten)
pars[6] = slope (assumed 0)
RESULT:
The chi^2 of the parameters.
COMMON BLOCKS:
RV_BLOCK - See exofast.pro for definition
MODIFICATION HISTORY
2012/06 -- Public release -- Jason Eastman (LCOGT)
(See exofast_getchi2_rv_fitnom.pro)
NAME:
exofast_getmcmcscale
PURPOSE:
Returns the optimal scale for MCMC chains by varying it until the
delta chi^2 = 1.
PROCEDURE:
Calculates the chi^2 of the best fit model. Then, for each
parameter, takes a small positive step (defined by seedscale) and
re-calculates the chi^2. It brackets the value that yields delta
chi^2 = 1 by doubling the seedscale until delta chi^2 > 1, then
does a binary search to find the value that is as close as
possible to delta chi^2 = 1. It repeats with a small negative
step. The average of the positive and negative excursion is
returned as the optimal scale.
If the best fit is accurate and the errors are
gaussian, uncorrelated, and accurately describe the data
(chi^2/dof = 1), this will result in the optimal acceptance rate
for MCMC fitting, 44% for 1 parameter, approaching 23% as the
number of parameters is large.
INPUTS:
BESTPARS - An NPARS element array of the best fit parameters
CHI2FUNC - A string of the named function that calculates the
chi^2 of your model given the parameters specified by
BESTPARS.
OPTIONAL INPUTS:
SEEDSCALE - An NPARS arrray that contains the small step that is
initially taken. If the chi^2 surface is a smooth,
monotonically increasing surface, then it only
affects execution time of the program. Since the
chi^2 surface is rarely so forgiving, it's better to
err on the side of too small. The default is 1d-3 for
each parameter.
BESTCHI2 - The chi^2 of the best fit parameters. If not given, it
is calculated from BESTPARS using CHI2FUNC.
ANGULAR - An array of indices (with up to NPARS elements)
indicating which elements are angular values. If
specified, and the step sizes determined are larger
than 2*pi, the step size will be set to pi.
NOTE: angular values must be in radians.
OPTIONAL KEYWORDS:
DEBUG - If set, it will print progress to the screen that
will help identify problems.
REVISION HISTORY:
2009/11/23 - Written: Jason Eastman - The Ohio State University
2013/02/26 - Change print statement to message statement
(See exofast_getmcmcscale.pro)
NAME:
EXOFAST_GETPHASE
PURPOSE:
Calculates the phase (mean anomaly/(2pi)) corresponding to a
particular true anomaly.
CALLING SEQUENCE:
phase = exofast_getphase(eccen,omega,trueanom [,/PRIMARY]$
[,/SECONDARY][,/L4][,/L5][,/PERIASTRON])
EXAMPLE:
;; find the phase of a primary transit described by e and omega
phase = exofast_getphase(e,omega,/primary)
INPUTS:
ECCEN: The eccentricity of the orbit, must be between 0 and 1
OMEGA: The argument of periastron of the star, in radians
TRUEANOM: The true anomaly (in radians) for which the mean anomaly
is desired.
OPTIONAL KEYWORDS:
PERIASTRON: Find the phase of Periastron (0 by definition)
PRIMARY: Find the phase of the primary transit
SECONDARY: Find the phase of the secondary eclipse
L4: Find the phase of the L4 Langrangian point (leading)
L5: Find the phase of the L5 Langrangian point (trailing)
ASCENDINGNODE: Find the phase of Ascending node (maximum RV signal)
DESCENDINGNODE: Find the phase of the Descending node (minimum RV signal)
DEPENDENCIES:
EXOFAST_KEPLEREQ.PRO
Modification History:
2010/06 - Rewritten by Jason Eastman (OSU)
(See exofast_getphase.pro)
NAME:
exofast_keplereq
PURPOSE:
Solve Kepler's Equation
DESCRIPTION:
Solve Kepler's Equation. Method by S. Mikkola (1987) Celestial
Mechanics, 40 , 329-334.
result from Mikkola then used as starting value for
Newton-Raphson iteration to extend the applicability of this
function to higher eccentricities
CATEGORY:
Celestial Mechanics
CALLING SEQUENCE:
eccanom=exofast_keplereq(m,ecc)
INPUTS:
m - Mean anomaly (radians; can be an array)
ecc - Eccentricity
OPTIONAL INPUT PARAMETERS:
KEYWORD INPUT PARAMETERS:
thresh: stopping criterion for the Newton Raphson iteration; the
iteration stops once abs(E-Eold)
(See exofast_keplereq.pro)
EXOFAST_LATEXTAB
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
EXOFAST_LATEXTAB
PURPOSE:
Prints the latex source to create a delux table of parameter values and
their errors. The errors are automatically rounded to two
significant digits and the values are rounded to the last
significant digit of the error.
This is intended to reduce typos and alleviate the tedium of making
tables.
CALLING SEQUENCE:
exofast_latextab, array, ['TEXFILE', PARNAMES=,
UNITS=,TITLE=,LABEL=, CAPTION=,SIDELABLES=,ORDER]
INPUTS:
ARRAY - A 2 or 3 column by NPARS rows array with the value,
+error, -error (e.g., the output from EXOFAST_PLOTDIST). If
only two columns are specified, all errors are assumed
to be symmetric.
OPTIONAL INPUTS:
TEXFILE - If specified, the program will write to this named file
instead of the screen
PARNAMES - If specified, the table will include the parameter
names as specified by this keyword. Otherwise they will
be left blank. Parameter names will be bracketed by
"$" so they are written in math mode.
UNITS - If specified, the table will include this
discriptor/units next to the parameters. Otherwise they will
be left blank.
TITLE - The title of the table
LABEL - The label for the table
CAPTION - The caption of the table
ORDER - An array of indices that specify the order of
outputs. Default is the input order. If SIDELABLES
are used, -1 should be inserted in each place the
next label should be printed.
SIDELABELS - An array of side labels to print. The number of
elements must match the number of -1s in the ORDER
array (this is where each lable will be inserted).
OUTPUTS:
Prints the latex source code to generate a table of values.
Modification History:
2011/03/10 - Written by Jason Eastman (OSU)
(See exofast_latextab.pro)
EXOFAST_LOMBSCARGLE
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
EXOFAST_LOMBSCARGLE
PURPOSE:
This function returns the NMIN best periods determined by a
lomb-scargle periodogram.
CALLING SEQUENCE:
periods = exofast_lombscargle(time, rv [BESTPARS=,
SCALE=, NMIN=, MINPERIOD=, MAXPERIOD=,NOSLOPE=, /NYQUIST, /PLOT])
INPUTS:
TIME - A time (eg, BJD_TDB) for each of the RV data points
RV - The RV data points
OPTIONAL OUTPUTS:
BESTPARS - The best (analytically fit) parameters determined from
the lomb-scargle periodogram (assumes circular orbit).
SCALE - The scale (~error) of the period, determined by the
spacing of the sampling of the periodogram.
NMIN - The number of periods corresponding to chi^2 minima
to return (default=5)
MINPERIOD - The minimum period to explore. Default = 1 unit of
time (typically days). This value is not used if
/NYQUIST is set.
MAXPERIOD - The maximum period to explore. Default is the full
range of dates.
OPTIONAL KEYWORDS
NOSLOPE - If set, no slope is assumed. Otherwise, a linear
term is analytically fit to the RV data.
NYQUIST - If set, the minimum period explored is the nyquist
frequency (half the minimum spacing between two
times). This takes precedence over MINPERIOD. Unless
the times are uniformly spaced (which is not
typical), use of this keyword is not recommended.
PLOT - A diagnostic tool. If set, the lomb-scargle
periodogram will be plotted, with the best period(s)
overplotted. This is useful for determining good
bounds, or if the observed behavior is unexpected.
DEPENDENCIES:
BUIELIB (http://www.boulder.swri.edu/~buie/idl/)
OUTPUTS:
result - An NMIN elements array containing the best periods
determined from the lomb-scargle periodogram.
MODIFICATION HISTORY
2009/05/01 -- Jason Eastman (Ohio State University)
(See exofast_lombscargle.pro)
EXOFAST_OCCULTQUAD
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
EXOFAST_OCCULTQUAD
PURPOSE:
This routine computes the lightcurve for occultation of a
quadratically limb-darkened source without microlensing. Please
cite Mandel & Agol (2002) and Eastman et al., (2013) if you make use
of this routine in your research. Please report errors or bugs to
agol@astro.washington.edu and jeastman@lcogt.net
Limb darkening has the form:
I(r)=[1-u1*(1-sqrt(1-(r/rs)^2))-u2*(1-sqrt(1-(r/rs)^2))^2]/(1-u1/3-u2/6)/pi
CALLING SEQUENCE:
occultquad, z0, u1, u2, p0, muo1, mu0, d=d
INPUTS:
z0 - impact parameter in units of rs
u1 - linear limb-darkening coefficient (gamma_1 in paper)
u2 - quadratic limb-darkening coefficient (gamma_2 in paper)
p0 - occulting star size in units of rs
OUTPUTS:
muo1 - fraction of flux at each z0 for a limb-darkened source
OPTIONAL OUTPUTS:
mu0 - fraction of flux at each z0 for a uniform source
d - The coefficients required to analytically calculate the
limb darkening parameters (see Eastman et al, 2013). For
backward compatibility, u1 and u2 are required, but not
necessary if d is used.
EXAMPLES:
; Calculate the same geometric transit with two different sets of
; limb darkening coefficients
p = 0.1d0
b = 0.5d0
x = (dindgen(300)/299d0 - 0.5d0)*2d0
z = sqrt(x^2 + b^2)
u1 = 0.25d0
u1 = 0.75d0
exofast_occultquad, z, u1, u2, p, muo1, muo, d=d
MODIFICATION HISTORY
2002 -- Eric Agol 2002
2009/04/06 -- Eric Agol (University of Washington) and
Jason Eastman (Ohio State University)
fixed bugs for p > 1 (e.g., secondary eclipses)
used double precision constants throughout
replaced rj,rc,rf with ellpic_bulirsch
faster, more robust integration for specialized case of
elliptic integral of the third kind
vectorized
more efficient case handling
combined ellk and ellec into ellke for speed
200x speed improvement over previous IDL code in typical case
allow negative values of p (anti-transit) to avoid Lucy-Sweeney like bias
(See exofast_occultquad.pro)
EXOFAST_PLOTDIST
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
exofast_plotdist
PURPOSE:
Plots the distributions from MCMC chains (Probability Distribution
Functions (PDFs) and covariances; returns median values and 68%
confidence intervals.
DESCRIPTION:
CALLING SEQUENCE:
exofast_plotdist, pars [,MEDIANPARS, ANGULAR=, PARNAMES=, UNITS=,
/NOCOVAR, PDFNAME=, COVARNAME=, BESTPARS=, PROBS=, /DEGREES]
INPUTS:
PARS - Array of NPARS x NSTEPS for each step in a Markov
Chain (with all chains concatenated together).
OPTIONAL INPUTS:
ANGULAR - An array that indexes the parameters that are
angular (radians unless /DEGREES is set). This will enable
special ways to calculate the median and plots, which
may otherwise fail. Default is none.
PARNAMES - An NPARS string array containing the parameter names
to be plotted on the axes.
HINT: TeXtoIDL is awesome.
UNITS - An NPARS string array containing the units to be
displayed on the X-axis of the PDFs.
PDFNAME - The filename of the output postscript plot of
PDFs. Default is pdf.ps
COVARNAME - The filename of the output postscript plot of
covariances. Default is covar.ps.
BESTPARS - An NPARS array that specifies the "best" parameters
of each value. If specified, a line will be drawn on
each PDF and a dot will be drawn on each covariance
plot at these values (this doesn't have to be
the best values).
PROBS - Array of probability contours to draw on the covariance
plots. Default is erf([1d,2d]/sqrt(2d0)) ~ [0.68,0.95].
OPTIONAL KEYWORDS:
DEGREES - If set, ANGULAR parameters are in degrees, not radians.
NOCOVAR - If set, no covariances plots will be made. Since
there will be NPARS*(NPARS-1)/2 plots, this can be
time consuming.
OUTPUTS:
MEDIANPARS - A 3 x NPARS array of median values, +err, -err (68%
confidence interval)
PDFNAME - A postscript plot of Probablity Distribution
Functions (PDFs)
COVARNAME - A postscript plot of covariances.
REVISION HISTORY:
2009/11/24 - Written: Jason Eastman - The Ohio State University
(See exofast_plotdist.pro)
EXOFAST_PREFITRV
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
EXOFAST_PREFITRV
PURPOSE:
Determines the global fit of RV data for a single planet.
CALLING SEQUENCE:
pars = exofast_prefitrv([/CHI2PERDOF,/CIRCULAR,CHI2=,])
OPTIONAL INPUTS:
MINPERIOD - The minimum period to consider (default=1 day)
MAXPERIOD - The maxmimum period to consider (default=range of
input times)
NMIN - The number of local minimums (from the Lomb-Scargle
periodogram) to fit fully. This should be increased
for highly eccentric planets. Default 5. Only the
parameters for the best fit (lowest chi2) are returned.
OPTIONAL OUTPUTS:
CHI2 - The chi2 of the fit (unaffected by /SCALEERR)
CHI2PERDOF - The chi2/dof of the fit (unaffected by /SCALEERR)
OPTIONAL KEYWORDS:
CIRCULAR - If set, the RV orbit will be assumed circulr (e=0)
NOSLOPE - If set, no slope will be fit to the RV data
SCALEERR - If set, the errors will be scaled such that the
probability of the chi^2 is 0.5.
PLOT - If set, plots of the best fit model and Lomb Scargle
periodograms will be displayed (useful for debgging)
RESULT:
The best-fit parameters for the RV fit of a single planet.
pars[0] = time of transit center
pars[1] = period
pars[2] = e*cos(omega)
pars[3] = e*sin(omega)
pars[4] = K
pars[5] = gamma
pars[6] = slope
COMMON BLOCKS:
RV_BLOCK - See exofast.pro for definition
MODIFICATION HISTORY
2012/06 -- Public release -- Jason Eastman (LCOGT)
2012/12 -- Add period, perr keywords to fix period
2013/01 -- Add error checking if no best fit is found
2013/02 -- Add degree of freedom when period is given
(See exofast_prefitrv.pro)
EXOFAST_RANDOM
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
EXOFAST_RANDOM
PURPOSE:
EXOFAST_RANDOM is a wrapper for pg_ran, the IDL translation
of ran from "Numerical Recipes Third Edition", page 342 by
Paolo C. Grigis, to make it a drop in replacement of RANDOMU
for normal or uniform deviates.
http://hea-www.cfa.harvard.edu/~pgrigis/idl_stuff/pg_ran.pro
While statistically, it is a far superior routine to IDL's
built-in, ran1 based, RANDOMU, it is also substantially (120x)
slower, as random number generators require serial
operations, which IDL does poorly.
If sequences larger than ~10^7 are required, IDL's built-in
routine is insufficient.
Call with a NEGATIVE or UNDEFINED SEED TO INITIALIZE
thereafter, do not alter seed between successive deviates in
a sequence.
CALLING SEQUENCE:
Result = EXOFAST_RANDOM(Seed [, D1, ..., Dn,/normal])
INPUTS:
Seed - A named variable containing the seed value for
random number generation. Seed is updated
once for each random number generated. The initial
value of Seed should be set to different NEGATIVE
values in order to obtain different random
sequences. An undefined variable name can be
specified and it will begin the seed with
-systime(/seconds). It is an error not to specify a
seed.
OPTIONAL INPUTS:
Di - The dimensions of the result. The dimension
parameters can be any scalar expression. Up to eight
dimensions can be specified. If no dimensions are
specified, RAN2 returns a scalar result
OPTIONAL KEYWORDS:
NORMAL - Set this keyword to generate random normal deviates
which uses an IDL translation of the Numerical
Recipes routine gasdev with the uniform deviates
generated from pg_ran.
OUTPUTS:
This function returns uniform random deviates between 0.0 and
1.0 (exclusive of the endpoint values) or normal deviates if
/NORMAL is set, with array dimensions given by the input
parameters, D1, ..., Dn.
MODIFICATION HISTORY:
Written by: Jason Eastman, OSU 3/2/2009
Based on Numerical Recipes Third Edition page 342
(See exofast_random.pro)
EXOFAST_ROSSITER
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
EXOFAST_ROSSITER
PURPOSE:
This routine implements the Ohta et al 2005 approximation for the
Rossiter-Mclaughlin effect, which assumes linear limb darkening
and is accurate to ~1 m/s for slowly rotating stars (see Hirano 2011).
CALLING SEQUENCE:
exofast_rossiter, x, y, u1, p, vsini, lambda, deltarv, z=z
INPUTS:
x - array of star-planet distances in the x direction
(stellar radii).
y - array or scalar of star-planet distance(s) in the y
direction. ie, the impact parameter (stellar radii)
u1 - linear limb-darkening coefficient -- epsilon in paper
p - occulting star size (stellar radii) -- gamma in paper
vsini - projected stellar surface velocity (same as output units)
lambda - projected spin-orbit angle (radians)
OPTIONAL INPUTS:
z - array of star-planet distances in the z direction.
-- Positive => in front of star (RM effect)
-- Negative => behind star (no effect)
-- assumed positive if not specified
OUTPUTS:
deltarv - Approximate Change in RV at each z due to the Rossiter
Mclaughlin effect (units same as vsini)
*** accurate to ~1 m/s ***
MODIFICATION HISTORY
2009/04/30 -- Jason Eastman (Ohio State University)
Loose translation from Scott Gaudi's (OSU) fortran code
2012/08 -- Jason Easmtan (LCOGT)
Make sign convention ahere to standard (now consistent with EXOFAST_GETB).
(See exofast_rossiter.pro)
EXOFAST_RV
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
exofast_rv
PURPOSE:
This function returns the radial velocity at each BJD. Optionally
includes the Rossiter McLaughlin effect, using the Ohta et al 2005
approximation good to ~1 m/s.
CALLING SEQUENCE:
result = exofast_rv(bjd, TPeriastron, period, V0, K, e, omega, $
i=i, a=a, u1=u1, p=p, vsini=vsini, lambda=lambda, $
r_star=r_star)
INPUTS:
bjd - Barycentric Julians dates for desired impact parameters
tperiastron - periastron passage time (BJD)
Period - Period of orbit (days)
V0 - systemic velocity (m/s)
K - radial velocity semi-amplitude (m/s)
OPTIONAL INPUTS:
e - eccentricity of orbit (0 if not specified)
omega - orbit's argument of periastron (radians)
- required if e is specified
- assumed to be pi/2 if e not specified
slope - the slope of the RV signal
t0 - the time at which V0 is referenced, if a slope is
given. If not given, (maxtime + mintime)/2 is used
OPTIONAL KEYWORDS:
rossiter - If set, calculates the RM effect. All values below
must also be specified.
Optional parameters (required for Rossiter-McLaughlin effect)
i - inclination of orbit (radians)
a - semi-major axis (R_star)
u1 - linear limb-darkening coefficient
p - occulting star size (stellar radii)
vsini - projected stellar surface velocity (m/s)
lambda - projected spin-orbit angle (radians)
OUTPUTS:
result - the radial velocity of the star at each point
MODIFICATION HISTORY
2009/05/01 -- Jason Eastman (Ohio State University)
(See exofast_rv.pro)
EXOFAST_TWILIGHT
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
exofast_twilight
PURPOSE:
Calculate the time of twilight (sunrise/sunset) to ~2 minute
accuracy (as compared to sunpos/eq2hor, which claims 2.2"/1"
accuracy, respectively).
Algorithm from http://williams.best.vwh.net/sunrise_sunset_algorithm.htm
Almanac for Computers, 1990
published by Nautical Almanac Office
United States Naval Observatory
Washington, DC 20392
INPUTS:
JD - The Julian date of days to calculate the twilight, in
UT. Scalar or vector.
OBSNAME - Name of the observatory (calls
observatory.pro). Either OBSNAME or
LATITUDE/LONGITUDE must be specified.
LATITUDE - Latitude of the observatory, in degrees.
LONGITUDE - Longitude of the observatory, in degrees West.
OPTIONAL INPUTS:
HORIZON - The Sun altitude that is considered twilight, in
degrees.
Civil Twilight = -6 degrees.
Nautical Twilight = -12 degrees. (Default)
Astronomical Twilight = -18 degrees.
TZ - The timezone of the observatory, in hours west of
Greenwich. Ignored unless /LOCAL is set.
OPTIONAL KEYWORDS:
SUNSET - Return the time of sunset instead of sunrise.
LOCAL - Return the result in local time. TZ must be specified
if this keyword is set.
NOTE: This does not take into account DST.
RESULT:
Returns the julian date of sunrise (or sunset, if /SUNSET is
set). If the sun is always down for the particular date, the
result will be -1. If the sun is always up for the particular
date, the result will be -2.
REVISION HISTORY:
2011/01/27 - Written: Jason Eastman - The Ohio State University
(See exofast_twilight.pro)
KTOM2
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
KTOM2
PURPOSE:
Determine the mass of the secondary given the velocity
semi-amplitude, eccentricity, inclination, Period, and mass of
the primary.
DESCRIPTION:
Determines the mass of the secondary by analytically solving the
cubic equation (only for the real root):
m2^3 - (m1+m2)^2*(K*sqrt(1d0-e^2)/sin(i))^3*Period/(2d0*!dpi*G) = 0
Inputs may be vectors or scalars, but must have the same number of
elements.
CALLING SEQUENCE:
m2 = ktom2(K, e, i, P, m1)
INPUTS:
K - The velocity semi-amplitude (m/s)
e - The eccentricity of the orbit
i - inclination of orbit (radians). To calculate the
minimum mass (msini, assuming M1 >> M2), use
i=!dpi/2d0.
Period - Period of orbit (days)
M1 - The mass of the Primary (m_sun)
OUTPUTS:
M2 - The mass of the secondary (m_sun)
MODIFICATION HISTORY
2011/08 -- Jason Eastman (Ohio State University)
(See ktom2.pro)
LINLD
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
LINLD
PURPOSE:
Interpolates the linear limb darkening tables of Claret and
Bloemen (2011). http://adsabs.harvard.edu/abs/2011A%26A...529A..75C
DESCRIPTION:
Loads an IDL save file found in $EXOFAST_PATH/quadld/, then
does a 3D linear interpolation of the table.
CALLING SEQUENCE:
coeffs = linld(logg, teff, feh, band [,MODEL=, METHOD=, VT=]);
INPUTS:
LOGG - The log of the stellar surface gravity
TEFF - The stellar effective temperature
FEH - The stellar metalicity
BAND - The observed bandpass. Allowed values are those defined in
Claret and Bloemen:
U,B,V,R,I,J,H,K, (Johnson/Cousins)
u',g',r',i',z', (Sloan)
Kepler, CoRoT,
Spitzer 3.6 um, Spitzer 4.5 um, Spitzer 5.8 um Spitzer 8.0 um,
u,b,v,y (Stromgren)
OPTIONAL INPUTS:
MODEL - The atmospheric model used to determine the limb
darkening values. Choose ATLAS or PHOENIX (default ATLAS).
METHOD - The method used. Choose L or F (default L)
VT - The microturbulent velocity (0,2,4,or 8, default 2)
RESULT:
The linear limb darkening parameter
COMMON BLOCKS:
LINLD_BLOCK - This is a self-contained block that stores the
contents of the IDL save files. This common block saves
the expensive step of restoring the same save files
for repeated calls (e.g., during and MCMC fit).
MODIFICATION HISTORY
2012/06 -- Public release -- Jason Eastman (LCOGT)
2013/01 -- Changed save filenames so they're not case sensitive
(now works with OSX)
Thanks Stefan Hippler
(See linld.pro)
MASSRADIUS_TORRES
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
MASSRADIUS_TORRES
PURPOSE:
Uses the Torres relation
(http://adsabs.harvard.edu/abs/2010A%26ARv..18...67T) to determine
the stellar mass and radius given the logg, Teff, and [Fe/H].
NOTE the errors in the empirical model of ulogm = 0.027d0, ulogr = 0.014d0
CALLING SEQUENCE:
massradius_torres, logg, teff, feh, mstar, rstar
INPUTS:
LOGG - The log of the stellar surface gravity
TEFF - The stellar effective temperature
FEH - The stellar metalicity
OUTPUTS:
MSTAR - The stellar mass, in solar masses
RSTAR - The stellar radius, in solar radii
MODIFICATION HISTORY
2012/06 -- Public release -- Jason Eastman (LCOGT)
(See massradius_torres.pro)
MYSETDIFF
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
MYSETDIFF
PURPOSE:
Finds the difference between two sets. Similar to
David Fanning's SETDIFFERENCE, but preserves the order and
duplicates at a substantial cost in speed.
CALLING SEQUENCE:
setdiff = mysetdiff(vec1, vec2)
INPUTS:
VEC1 - An array of elements
VEC2 - Another array of elements
RESULT:
The difference between VEC1 and VEC2.
EXAMPLE:
vec1 = indgen(5)
vec2 = indgen(5) + 3
print, mysetdiff(vec1,vec2)
0 1 2
MODIFICATION HISTORY
2012/06 -- Public release -- Jason Eastman (LCOGT)
(See mysetdiff.pro)
NONLINLD
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
NONLINLD
PURPOSE:
Interpolates the non-linear limb darkening tables of Claret and
Bloemen (2011). http://adsabs.harvard.edu/abs/2011A%26A...529A..75C
DESCRIPTION:
Loads an IDL save file found in $EXOFAST_PATH/quadld/, then
does a 3D linear interpolation of the table.
CALLING SEQUENCE:
coeffs = nonlinld(logg, teff, feh, band [,MODEL=, METHOD=, VT=]);
INPUTS:
LOGG - The log of the stellar surface gravity
TEFF - The stellar effective temperature
FEH - The stellar metalicity
BAND - The observed bandpass. Allowed values are those defined in
Claret and Bloemen:
U,B,V,R,I,J,H,K, (Johnson/Cousins)
u',g',r',i',z', (Sloan)
Kepler, CoRoT,
Spitzer 3.6 um, Spitzer 4.5 um, Spitzer 5.8 um Spitzer 8.0 um,
u,b,v,y (Stromgren)
OPTIONAL INPUTS:
MODEL - The atmospheric model used to determine the limb
darkening values. Choose ATLAS or PHOENIX (default ATLAS).
METHOD - The method used. Choose L or F (default L)
VT - The microturbulent velocity (0,2,4,or 8, default 2)
RESULT:
The non-linear limb darkening parameters
COMMON BLOCKS:
NONLINLD_BLOCK - This is a self-contained block that stores the
contents of the IDL save files. This common block saves
the expensive step of restoring the same save files
for repeated calls (e.g., during and MCMC fit).
MODIFICATION HISTORY
2012/06 -- Public release -- Jason Eastman (LCOGT)
2013/01 -- Changed save filenames so they're not case sensitive
(now works with OSX) -- thanks Stefan Hippler.
(See nonlinld.pro)
OCCULTQUAD_EXTERN
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
OCCULTQUAD_EXTERN
PURPOSE:
This routine computes the lightcurve for occultation of a
quadratically limb-darkened source without microlensing.
DESCRIPTION:
This is an IDL wrapper for the fortran version of the code, and is
a 2-5x faster, drop-in replacement for EXOFAST_OCCULTQUAD; see
complete documentation there. The shared libraries must be
compiled for your system before you can use this. See README file
in this directory.
MODIFICATION HISTORY
2012/06 -- Jason Eastman, Public release
(See other/occultquad_extern.pro)
QUADLD
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
QUADLD
PURPOSE:
Interpolates the quadratic limb darkening tables of Claret and
Bloemen (2011). http://adsabs.harvard.edu/abs/2011A%26A...529A..75C
DESCRIPTION:
Loads an IDL save file found in $EXOFAST_PATH/quadld/, then
does a 3D linear interpolation of the table.
CALLING SEQUENCE:
coeffs = quadld(logg, teff, feh, band [,MODEL=, METHOD=, VT=]);
INPUTS:
LOGG - The log of the stellar surface gravity
TEFF - The stellar effective temperature
FEH - The stellar metalicity
BAND - The observed bandpass. Allowed values are those defined in
Claret and Bloemen:
U,B,V,R,I,J,H,K, (Johnson/Cousins)
u',g',r',i',z', (Sloan)
Kepler, CoRoT,
Spitzer 3.6 um, Spitzer 4.5 um, Spitzer 5.8 um Spitzer 8.0 um,
u,b,v,y (Stromgren)
OPTIONAL INPUTS:
MODEL - The atmospheric model used to determine the limb
darkening values. Choose ATLAS or PHOENIX (default ATLAS).
METHOD - The method used. Choose L or F (default L)
VT - The microturbulent velocity (0,2,4,or 8, default 2)
RESULT:
The quadradic limb darkening parameters
COMMON BLOCKS:
LD_BLOCK - This is a self-contained block that stores the
contents of the IDL save files. This common block saves
the expensive step of restoring the same save files
for repeated calls (e.g., during and MCMC fit).
MODIFICATION HISTORY
2012/06 -- Public release -- Jason Eastman (LCOGT)
2013/01 -- Changed save filenames so they're not case sensitive
(now works with OSX) -- thanks Stefan Hippler.
(See quadld.pro)
READEXO
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
READEXO
PURPOSE:
Reads the http://exoplanets.org/csv-files/exoplanets.csv database into an
IDL structure for easy plotting and manipulation.
See http://exoplanets.org/help/common/data for the definitions
of each tag.
OPTIONAL KEYWORDS:
UPDATE - By default, if exoplanets.csv exists, it will not
retrieve a new copy. Set this keyword to wget the csv
file and overwrite exoplanets.csv. It will only overwrite
it if the copy on the server is newer. This program is
about 25% slower when this keyword is specified
unnecessarily and 2x slower when it actually needs to
update.
CALLING SEQUENCE:
result = readexo([/update])
EXAMPLE - Print all data for HAT-P-13 b:
planetdata = readexo(/update)
match = where(planetdata.name eq 'HAT-P-13 b')
names = tag_names(planetdata)
if match[0] ne -1 then for i=0, n_tags(planetdata)-1 do $
print, names[i], planetdata.(i)[match], format='(a15,a100)'
DEPENDENCIES:
wget
ENVIRONMENT VARIABLES:
EXOFAST_PATH - The environment variable specifying the directory
in which to save/look for the CSV file. If not set,
the program will use the current working directory.
NOTES:
All data types are strings for robustness and generality (no
update to this program is required if new rows or columns are
added to the database). Converting to numeric values is trivial,
but be sure to use the appropriate precision (eg, to convert ra
from hours to degrees, do radeg = planetdata.ra*15.d0, not radeg =
planetdata.ra*15).
Be aware that the order of tags or even the tagnames may change if
the database changes.
REVISION HISTORY
Created: 2011/01/19 Jason Eastman, OSU
2012/06/15 Updated for change to exoplanets.org.
2013/01/14 Updated for change to URL.
(See readexo.pro)
RECENTER
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
RECENTER
PURPOSE:
Recenters a distribution of periodic parameters to the domain
(mode-period/2,mode+period/2]
CALLING SEQUENCE:
par = recenter(par,period)
INPUTS:
PAR = An array of parameters
PERIOD = The period of the parameter distribution
EXAMPLE:
MODIFICATION HISTORY
2012/06 -- Jason Eastman (LCOGT)
(See recenter.pro)
ROUND5
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
ROUND5
PURPOSE:
Rounds a value to the nearest 5 in the second significant
digit. Useful for plotting.
CALLING SEQUENCE:
result = round5(value)
INPUTS:
VALUE - A number to be rounded
RESULT:
The rounded number
MODIFICATION HISTORY
2012/06 -- Public release -- Jason Eastman (LCOGT)
(See round5.pro)
TRANSITGIF
[Previous Routine]
[List of Routines]
NAME:
TRANSITGIF
PURPOSE:
Creates an animated gif of a transit from exoplanets.org.
CALLING SEQUENCE:
TRANSITGIF, ['pname',PARS=,BAND=,CADENCE=,PARS=,YMIN=,BAK=,
/YELLOW,/DISPLAY,/SAVEFRAMES,/UPDATE]
OPTIONAL INPUTS:
PNAME - The name of the planet from exoplanets.org (not case
sensitive). The output file will be the planet name
stripped of spaces + .gif.
PARS - A 10 element array of custom parameters with which to
make a transit animation:
PARS[0] = period (days)
PARS[1] = eccentricity
PARS[2] = omega_star (radians)
PARS[3] = inclination (radians)
PARS[4] = a/Rstar
PARS[5] = Rp/Rstar
PARS[6] = logg (cgs)
PARS[7] = Teff (K)
PARS[8] = [Fe/H]
PARS[9] = Rstar (R_sun)
NOTE: PNAME or PARS must be specified. If PARS is
specified and PNAME is not, the output will be "transit.gif".
BAND - The bandpass for the limb darkening calculation (see
quadld.pro for allowed values). Default = 'Sloanz'
CADENCE - A frame is generated every CADENCE minutes for a
duration of 2*T_FWHM*MAXSTAR/RSTAR. This duration
forces the X axis of the plot to be on the same scale as
the drawing, ignoring the curvature of the planet's orbit.
Default = 5.
NOTE: A constant cadence between animations means the
duration of the animation (number of frames)
corresponds to the duration of the transit.
YMIN - The minimum of transit plot. Default is 0.985
NOTE: A constant minimum makes the plots to scale from
animation to animation, but this default cuts off some
deeper planets.
BAK - A 3-element array containing the RGB values of the
background color. Default is [135,206,250]; (sky blue)
Use [255,255,255] for white.
MAXSTAR - The size of the maximum star, in R_sun. A star of this
size will be centered left/right and span the top half of
the image. The total image is 4*MAXSTAR
square. Keeping MAXSTAR the same (positive value)
for multiple animations means they will be to
scale. Default = 1.5.
Numbers between -1 and 0 are interpreted as a
fractions of RSTAR (e.g, -1 ensures the star will take
up the entire range; -0.9 leaves a nice margin). Negative
numbers are useful if you're making a single
animation and want it to fit nicely.
NOTE: If the Rstar is larger than MAXSTAR, it will
produce an error.
OPTIONAL KEYWORDS:
YELLOW - Makes the star yellow and the planet black, instead of
their color temperatures.
NOTE: Output file is then PNAME.yellow.gif
DISPLAY - Spawns firefox to display the animated gif.
SAVEFRAMES - The individual PNG frames are deleted upon
completion. Set this keyword to keep them.
UPDATE - Update the local copy of exoplanets.csv
OUTPUTS:
pname.gif - An animated gif of the transit.
EXAMPLES:
;; creates an animation of HAT-P-3b, called 'HAT-P-3b.gif'
transitgif, 'HAT-P-3b'
;; creates an animation Earth, called 'transit.gif'
pars = [365.242199d0,0.01671123d0,!dpi/2d0,!dpi/2d0,215.094177d0,$
0.0091705248d0,4.43834417d0,5778d0,0d0,1d0]
transitgif, pars=pars
;; Creates animations for all transiting exoplanets from exoplanets.org
planets = readexo(/update)
transit = where(planets.transit,ntransits)
for i=0, ntransits-1 do $
transitgif, strcompress(planets.name[transit[i]],/remove_all),$
maxstar=max(double(planets.rstar[transit])), band='Sloanz',$
ymin=1d0-max(double(planets.depth[transit]))
DEPENDENCIES:
convert. While IDL can write animated gifs directly,
they're only 8 bit. In color, that creates unsightly
ringing, so we use convert instead.
EXOFAST (http://astroutils.astronomy.ohio-state.edu/exofast/)
MODIFICATION HISTORY
2012/06 -- Public Release, Jason Eastman (LCOGT)
(See transitgif.pro)