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: **Thu Sep 18 09:32:31 2014.

- BJD2TARGET
- ELLKE
- ELLPIC_BULIRSCH
- EXOFAST
- EXOFAST_AMOEBA
- EXOFAST_CHI2
- EXOFAST_DEMC
- EXOFAST_ERRELL
- EXOFAST_GELMANRUBIN
- EXOFAST_GETB
- EXOFAST_GETCHI2_RV_FIT
- EXOFAST_GETCHI2_RV_FITCIR
- EXOFAST_GETCHI2_RV_FITCIRNOM
- EXOFAST_GETCHI2_RV_FITNOM
- EXOFAST_GETMCMCSCALE
- EXOFAST_GETPHASE
- EXOFAST_KEPLEREQ
- EXOFAST_LATEXTAB
- EXOFAST_LOMBSCARGLE
- EXOFAST_OCCULTQUAD
- EXOFAST_PLOTDIST
- EXOFAST_PREFITRV
- EXOFAST_RANDOM
- EXOFAST_ROSSITER
- EXOFAST_RV
- EXOFAST_TWILIGHT
- KTOM2
- LINLD
- MASSRADIUS_TORRES
- MYSETDIFF
- NONLINLD
- OCCULTQUAD_EXTERN
- QUADLD
- READEXO
- RECENTER
- ROUND5
- TRANSITGIF

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)