Extended IDL Help

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.


List of Routines


Routine Descriptions

BJD2TARGET

[Next Routine] [List of Routines]
 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)


ELLKE

[Previous Routine] [Next Routine] [List of Routines]
 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)


ELLPIC_BULIRSCH

[Previous Routine] [Next Routine] [List of Routines]
 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)


EXOFAST

[Previous Routine] [Next Routine] [List of Routines]
 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)


EXOFAST_AMOEBA

[Previous Routine] [Next Routine] [List of Routines]
 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)


EXOFAST_CHI2

[Previous Routine] [Next Routine] [List of Routines]
 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)


EXOFAST_DEMC

[Previous Routine] [Next Routine] [List of Routines]
 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)


EXOFAST_ERRELL

[Previous Routine] [Next Routine] [List of Routines]
 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)


EXOFAST_GELMANRUBIN

[Previous Routine] [Next Routine] [List of Routines]
 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)


EXOFAST_GETB

[Previous Routine] [Next Routine] [List of Routines]
 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)


EXOFAST_GETCHI2_RV_FIT

[Previous Routine] [Next Routine] [List of Routines]
 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)


EXOFAST_GETCHI2_RV_FITCIR

[Previous Routine] [Next Routine] [List of Routines]
 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)


EXOFAST_GETCHI2_RV_FITCIRNOM

[Previous Routine] [Next Routine] [List of Routines]
 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)


EXOFAST_GETCHI2_RV_FITNOM

[Previous Routine] [Next Routine] [List of Routines]
 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)


EXOFAST_GETMCMCSCALE

[Previous Routine] [Next Routine] [List of Routines]
 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)


EXOFAST_GETPHASE

[Previous Routine] [Next Routine] [List of Routines]
 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)


EXOFAST_KEPLEREQ

[Previous Routine] [Next Routine] [List of Routines]
 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)