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: Sun Dec 14 16:46:57 2014.


List of Routines


Routine Descriptions

ALPHA_BEAM

[Next Routine] [List of Routines]
 NAME:
   alpha_beam
     
 PURPOSE:
   Returns the alpha_beam parameter for the Doppler beaming effect
   (See Loeb & Gaudi [2003] ApJ 588, 117.)

 CALLING SEQUENCE:
   alpha_beam(Ts)  

 INPUTS:
   Ts -- stellar effective temperature (K)

 RESTRICTIONS:
   Assumes blackbody radiation from star

 REQUIRES:
   bbrad_kep_dop -- http://www.lpl.arizona.edu/~bjackson/idl_code/index.html
   push -- http://www.lpl.arizona.edu/~bjackson/idl_code/index.html

 EXAMPLE:
   Ts = 6350.
   ;Should be 1.0855685
   print, alpha_beam(Ts)

 MODIFICATION HISTORY:
   2013 Jul 15 -- Written by Brian Jackson (decaelus@gmail.com)

(See alpha_beam.pro)


BBRAD_KEP_DOP

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   bbrad_kep_dop
     
 PURPOSE:
   Returns stellar blackbody flux (arbitrary units) integrated over Kepler
     bandpass and possibly Doppler shifted

 CALLING SEQUENCE:
   bbrad_kep_dop(temp, vz)

 INPUTS:
   temp -- stellar effective temperature (K)
   vz -- stellar Doppler velocity in units of c (< 0 means TOWARD observer)

 EXAMPLE:
   Ts = 6350.
   vz = 0.
   ;Should be 4.2696805e-08
   print, bbrad_kep_dop(Ts, vz)

 MODIFICATION HISTORY:
   2011 Sep 15 - Returns black body radiation convolved with Kepler response
     function (http://keplergo.arc.nasa.gov/CalibrationResponse.shtml) and 
     Doppler shifts
   2012 Jan 11 - Using expression from Loeb & Gaudi (2003) ApJL 588, L117.

(See bbrad_kep_dop.pro)


CLARETLDC

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   claretldc
     
 PURPOSE:
   Returns limb-darkening coefficients determined via linear interpolation from the
     tables of Claret & Bloemen (2011) A&A 529, 75.

 CALLING SEQUENCE:
   ldc = claretldc(Teff, logg, M, ld_law)

 INPUTS:
   Teff - stellar effective temp
   logg - log(surf grav)
   M - metallicity
   LD_law - integer indicating which limb darkening law to use, 
     either linear (0), quadratic (1), or non-linear (2)

 REQUIRES:
   push.pro -- http://www.lpl.arizona.edu/~bjackson/idl_code/

 EXAMPLE:
   ;solar values
   Teff = 5780.
   logg = 4.44
   M = 0.
   print, claretldc(Teff, logg, M, 1)

 MODIFICATION HISTORY:
   2013 Sep 23 -- Written by Brian Jackson (decaelus@gmail.com)

(See claretldc.pro)


COS_F

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	cos_f

 PURPOSE:
	This function returns the cosine of the true anomaly, using relations given in 
	Murray & Dermott (1999), ch. 2.4. This routine uses keplereq.pro.

 CATEGORY:
	Celestial Mechanics.

 CALLING SEQUENCE:
	Result = cos_f(mean_anom, ecc)

 INPUTS:
	mean_anom: orbital mean anomaly
	ecc: orbital eccentricity

 OUTPUTS:
	Cosine of the orbital true anomaly	

 EXAMPLE:
       ;make mean anomaly array
	num = 101
	mean_anom = 2.*!pi*dindgen(num)
	;orbital eccentricity of 0.1
	ecc = 0.1
	cos_f = cos_f(mean_anom, ecc)

 MODIFICATION HISTORY:
 	Written by:	Brian Jackson (decaelus@gmail.com), 2011 Jul 25.

(See cos_f.pro)


EEBLS

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	eebls

 PURPOSE:
     >>>>>>>>>>>> This routine computes BLS spectrum <<<<<<<<<<<<<<

     [ see Kovacs, Zucker & Mazeh 2002, A&A, Vol. 391, 369 ]

     This is the slightly modified version of the original BLS routine 
     by considering Edge Effect (EE) as suggested by 
     Peter R. McCullough [ pmcc@stsci.edu ].

     This modification was motivated by considering the cases when 
     the low state (the transit event) happened to be devided between 
     the first and last bins. In these rare cases the original BLS 
     yields lower detection efficiency because of the lower number of 
     data points in the bin(s) covering the low state. 

     For further comments/tests see  www.konkoly.hu/staff/kovacs.html

	N.B. This is an IDL version of the code originally created by 
	Kovacs et al. I (Brian Jackson - decaelus@gmail.com) take no credit at all
	for developing this routine.

 CATEGORY:
	Astrophysical

 CALLING SEQUENCE:
 
	eebls,n,t,x,u,v,nf,fmin,df,nb,qmi,qma,$
	  p,bper,bpow,depth,qtran,in1,in2

 INPUTS:
     n    = number of data points
     t    = array {t(i)}, containing the time values of the time series
     x    = array {x(i)}, containing the data values of the time series
     u    = temporal/work/dummy array, must be dimensioned in the 
    calling program in the same way as  {t(i)}
     v    = the same as  {u(i)}
     nf   = number of frequency points in which the spectrum is computed
     fmin = minimum frequency (MUST be > 0)
     df   = frequency step
     nb   = number of bins in the folded time series at any test period   
     qmi  = minimum fractional transit length to be tested
     qma  = maximum fractional transit length to be tested

 OUTPUTS:
     p    = array {p(i)}, containing the values of the BLS spectrum
    at the i-th frequency value -- the frequency values are 
    computed as  f = fmin + (i-1)*df
     bper = period at the highest peak in the frequency spectrum
     bpow = value of {p(i)} at the highest peak
     depth= depth of the transit at   *bper*
     qtran= fractional transit length  [ T_transit/bper ]
     in1  = bin index at the start of the transit [ 0 < in1 < nb+1 ]
     in2  = bin index at the end   of the transit [ 0 < in2 < nb+1 ]

 RESTRICTIONS:
     -- *fmin* MUST be greater than  *1/total time span* 
     -- *nb*   MUST be lower than  *nbmax* 
     -- Dimensions of arrays {y(i)} and {ibi(i)} MUST be greater than 
    or equal to  *nbmax*. 
     -- The lowest number of points allowed in a single bin is equal 
    to   MAX(minbin,qmi*N),  where   *qmi*  is the minimum transit 
    length/trial period,   *N*  is the total number of data points,  
    *minbin*  is the preset minimum number of the data points per 
    bin.

 EXAMPLE:
 ;This is about 1 month of Kepler data
 n = 1440.
 ;Orbital period of ~5 hrs
 per = 20. + randomu(seed) - 0.5
 t = dindgen(n)
 x = replicate(1.0, n)
 u = dblarr(n)
 v = dblarr(n)
 ;According to Kovacs+ (2002), this is the relationship between nf, meant to 
 ;  represent the nber of independent frequencies to test, and n
 nf = fix(n^(0.83))
 fmin = 2./(max(t)-min(t))
 fmax = 2./(t[1]-t[0])
 df = (1./(0.05*per)-1/(5.*per))/(nf*10)
 ;must be less than 2000
 nb = 50
 folded_time = t mod per
 mid = randomu(seed)*per
 ;fake transit
 width = 0.03*per
 depth = 100d-6
 ind = where(abs(folded_time-mid) le width)
 fake_trans = replicate(1.0, n_elements(folded_time))
 fake_trans[ind] = 1.0-depth
 x = fake_trans
 qma = 0.05
 qmi = 0.01
 eebls,n,t,x,u,v,nf,fmin,df,nb,qmi,qma,$
   p,bper,bpow,depth,qtran,in1,in2
 freq = dindgen(nf)*df + fmin
 plot, 1./freq, p, /xlog, psym=2
 print, bper, per

 MODIFICATION HISTORY:
 	Original fortran code written by G. Kovacs, S. Zucker, & T. Mazeh
	Ported to IDL by Brian Jackson (decaelus@gmail.com), 2012 May

(See eebls.pro)


EVILMC

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	EVILMC

 PURPOSE:
	This function returns the ellipsoidal variation of a slowly-rotating star induced by a 
	low-mass companion. This code is publicly available with NO warranties whatsoever at 
	http://www.lpl.arizona.edu/~bjackson/code/idl.html. If you use the code, please cite the 
	following paper: Jackson et al. (2012) ApJ 750, 1. 

 CATEGORY:
	Astrophysics.

 CALLING SEQUENCE:
 
	Result = EVILMC(phs, params)

 INPUTS:
	phs:	orbital phase; phs = 0 corresponds to inferior conjunction (i.e. mid-transit)
		that positional parameters are shown with Initial Caps.
	params: array of ellipsoidal variation parameters, as follows --
  		params[0] - # of lat/long grid points on star, 
  		params[1] - mass ratio, q (= M_p/M_*)
		params[2] - K_z, stellar reflex velocity, in m/s
		params[3] - T_0, stellar effective temp in K
		params[4:6] - x, y, z of stellar rotation vector in units of mean motion
		params[7:8] - \gamma_i, stellar quadratic limb-darkening coefficients
		params[9] - \beta, stellar gravity darkening exponent
		params[10] - a, semi-major axis (in units of stellar radii)
		params[11] - orbital period (in days)
		params[12] - orbital inclination in degrees
		params[13] - eccentricity (assumed 0 for now)
		params[14] - longitude of ascending node (should probably be 0)
		params[15] - longitude of pericenter (while ecc = 0, this is assumed 0)

 OUTPUTS:
	Stellar disk emission in MKS units - Typically, the disk emission should be normalized
	by the disk emission at phs = 0.5, when a planet is eclipse by the star. See the 
	companion code, EVILMC_plphs.pro to see how this is done.

 RESTRICTIONS:
	No checking of parameters is done, so be sure everything is in the correct units, etc.
	Also, the code requires companion routines: orb_pos.pro. 
	These are available at http://www.lpl.arizona.edu/~bjackson/code/idl.html. Also, for now
	the code assumes an orbital eccentricity of 0 and blackbody radiation from the star (among
	other important assumptions).

 EXAMPLE:
	;See the companion code EVILMC_plphs.pro and example_EVILMC.pro for more help. Both are
	;available at http://www.lpl.arizona.edu/~bjackson/code/idl.html.

	q = 1.10d-3
	semi = 4.15 ;A/R_0
	per = 2.204733 ;days
	Omega = 4.73d-7 ;s^{-1}
	Ts = 6350. ;K
	Kz = 300. ;m/s
	logg = 4.07 ;log(cm/s^2)
	M = 0.26 ;[Fe/H]
	ecc = 0.0 ;orbital eccentricity
	asc_node = 0. ;longitude of planetary ascending node
	peri_long = 0. ;longitude of planetary pericenter
	inc = 83.1 ;orbital inclination in degrees
	bet = 0.0705696 ;gravity darkening exponent, (T/T_0) = (g/g_0)^\beta
	;limb-darkening coefficients, I(\mu)/I(1) = 1 - \gamma_1 (1-\mu) - \gamma_2 (1-\mu)^2
	gam1 = 0.314709
	gam2 = 0.312125
	gam = [gam1, gam2]
	;number of latitude or longitude grid points on stellar surface
	num_grid = 31
	;orientation (X, Y, Z) of stellar rotation axis
	Omegahat = [0., 0., 1.0]
	;stellar rotation axis
	Omega = Omega*Omegahat
	p = [num_grid, q, Kz, Ts, Omega, gam, bet, semi, per, inc, ecc, asc_node, peri_long]

	;orbital phase
	num_phs = 101
	phs = dindgen(num_phs)*1./(num_phs-1.)
	ret = EVILMC(phs, p)

	;calculate normalization
	norm_st_phs = EVILMC(0.5, p)
	ret /= norm_st_phs
	plot, phs, ret, yr=[min(ret), max(ret)]

 MODIFICATION HISTORY:
 	Written by:	Brian Jackson (decaelus@gmail.com, 2011 Oct 5.

(See EVILMC.pro)


EVILMC_PLPHS

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	EVILMC_plphs

 PURPOSE:
	Calls routines to calculate stellar ellipsoidal variation and planetary phase function. 
	Code is provided with no warranties whatsoever at 
	http://www.lpl.arizona.edu/~bjackson/code/idl.html. See also paper Jackson et al. (2012) 
	ApJ 750, 1 for more details.

 CATEGORY:
	Astrophysics.

 CALLING SEQUENCE:
	Result = EVILMC_plphs(phs, params)

 INPUTS:
       phs:    orbital phase; phs = 0 corresponds to inferior conjunction (i.e. mid-transit)
               that positional parameters are shown with Initial Caps.
       params: array of ellipsoidal variation parameters, as follows --
               params[0] - # of lat/long grid points on star, 
               params[1] - mass ratio, q (= M_p/M_*)
               params[2] - K_z, stellar reflex velocity, in m/s
               params[3] - T_0, stellar effective temp in K
               params[4:6] - x, y, z of stellar rotation vector in units of mean motion
               params[7:8] - \gamma_i, stellar quadratic limb-darkening coefficients
               params[9] - \beta, stellar gravity darkening exponent
               params[10] - a, semi-major axis (in units of stellar radii)
               params[11] - orbital period (in days)
               params[12] - orbital inclination in degrees
               params[13] - eccentricity (assumed 0 for now)
               params[14] - longitude of ascending node (should probably be 0)
               params[15] - longitude of pericenter (while ecc = 0, this is assumed 0)
               params[16] - F_0, planetary phase function parameter
               params[17] - F_1, planetary phase function parameter

 OUTPUTS:
	Normalized stellar ellipsoidal variation plus planetary phase function

 RESTRICTIONS:
	No checking of parameters is done, so be sure everything is in the correct units, etc.
	Also, the code requires a few companion routines: EVILMC.pro.
	These are available at http://www.lpl.arizona.edu/~bjackson/code/idl.html.

 EXAMPLE:

       See the companion code example_EVILMC.pro to see how to use this code.

 MODIFICATION HISTORY:
 	Written by:	Brian Jackson, 2012 April 18.

(See EVILMC_plphs.pro)


EXAMPLE_EVILMC

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	example_EVILMC

 PURPOSE:
	Provides an example system (HAT-P-7) to use to the EVIL-MC model. Code is provided with no
	warranties whatsoever at http://www.lpl.arizona.edu/~bjackson/code/idl.html. See also 
	paper Jackson et al. (2012) ApJ 750, 1 for more details.

 CATEGORY:
	Astrophysics.

 CALLING SEQUENCE:
	example_EVILMC

 INPUTS:
	none

 OUTPUTS:
	none

 SIDE EFFECTS:
	Makes a plot of the normalized ellipsoidal variation on the screen

 EXAMPLE:
	example_EVILMC

 MODIFICATION HISTORY:
 	Written by:	Brian Jackson (decaleus@gmail.com), 2012 April 18.

(See example_EVILMC.pro)


FIND_BETAS

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   find_betas
     
 PURPOSE:
   Estimates beta_mu parameters for MCMC analysis by determining acceptance
     fraction for trial betas spanning many orders of magnitude; 
     See Ford (2005) ApJ 129, 1706 for more

 CALLING SEQUENCE:
   find_betas, model_func, like_func, x, y, err, init_p, betamu[, $
     beta_pool=beta_pool, parinfo=parinfo, /verbose]

 INPUTS:
   model_func -- name of model function
   like_func -- function that returns log(likelihood) for jump
   x/y/err -- in/dependent variables and uncertainties
   init_p -- initial guesses for model parameters

 OUTPUTS:
   betamu -- the sought-for beta_mu's

 KEYWORD PARAMETERS:
   beta_pool -- beta's to try
   verbose -- prints out acceptance fraction for different trial betas; 
     If it's your first time running this routine for a problem, it's a good
     idea to activate this flag.
   num_tries -- number of links in chain to sample beta_mu (defaults to 100)

 RESTRICTIONS:
   This routine tries beta_mu's scaled by the initial guesses 
     from 1e-6 to 1.7 times the initial guesses. If the right betamu isn't in 
     that range, this routine won't find it.

 EXAMPLE:
   test_mh_gibbs provides an example

 MODIFICATION HISTORY:
   2011 Nov 27 -- Written by Brian Jackson (decaelus@gmail.com)
   2013 Apr 17 -- Modified by BJ to allow for like_func
   2013 Jun 4 -- Added print time

(See find_betas.pro)


GET_ASYMMETRIC_UNCERTAINTIES

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   get_asymmetric_uncertainties
     
 PURPOSE:
   Given a distribution of parameter values, this routine returns the 
     left/right uncertainties.

 CALLING SEQUENCE:
   errors = get_asymmetric_uncertainties(dst, num_sig=num_sig)

 INPUTS:
   dst -- distribution from which to extract uncertainties

 KEYWORD PARAMETERS:
   num_sig -- If the distribution is Gaussian about a mean value,
     this parameter represents the number of sigma away from that mean.
   sign -- If set, the routine will return a negative left error

 EXAMPLE:
  x = randomn(seed, 1001)
  errors = get_asymmetric_uncertainties(x, num_sig=1.)
  ;Should both be close to 1
  print, errors

 MODIFICATION HISTORY:
   2013 Jun 10 -- Written by Brian Jackson (decaelus@gmail.com)
   2013 Jul 31 -- Dealt with the case that there are NOT elements to the left or
     right of the mean value by returning min/max values
   2013 Jul 31 -- Added sign keyword

(See get_asymmetric_uncertainties.pro)


KEPCOTREND

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   kepcotrend
     
 PURPOSE:
   Cotrends Kepler lightcurves using the provided cotrend basis
     vectors in a way analogous to the kepcotrend python routine provided at 
     http://keplergo.arc.nasa.gov/PyKE.shtml. 

   The basis vector files required for cotrending are available here:
     http://archive.stsci.edu/kepler/cbv.html.

   Note that this routine doesn't provide for all the functionality that the
     Kepler team's kepcotrend.py routine does. In particular, this routine
     won't work for short-cadence data.

   The technique used here is discussed at length in 
     _Numerical Recipes in C_, 2nd ed., ch. 15. The CBV array below seems to be
     equivalent to the V array in NR.

 CALLING SEQUENCE:
   cotrended_flux = kepcotrend(lcfile, bvfile, listbv, model=model)

 INPUTS:
   lcfile -- long-cadence FITS file containing the data to be cotrended
   bvfile -- CBV file
   listbv -- array indicating which CBVs to use, starting at 0

 OUTPUT
   Cotrended flux, with bad data values (NaNs) still included

 OPTIONAL OUTPUT:
   model -- array containing the cotrending model used

 RESTRICTIONS:
   This routine isn't very clever; it only cotrends the data.

 EXAMPLE:
   lcfile = "kplr010666592-2009131105131_llc.fits"
   cbvfile = "kplr2009131105131-q00-d14_lcbv.fits"
   listbv = indgen(8)
   cotrended_flux = kepcotrend(lcfile, cbvfile, listbv)
   plot, indgen(n_elements(cotrended_flux)), cotrended_flux, psym=3, $
     yr=[min(cotrended_flux), max(cotrended_flux)], ystyle=2 

 MODIFICATION HISTORY:
   2012 Mar - Written by Brian Jackson (decaelus@gmail.com)
   2012 Dec 14 - Check for empty cbvs added by Nikole Thom (nklewis@mit.edu)

(See kepcotrend.pro)


KEPLEREQ

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
    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=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 keplereq.pro)


KEP_BBTEMP

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   Kep_BBtemp
     
 PURPOSE:
   Returns the blackbody temperature that corresponds to a Kepler eclipse

 CALLING SEQUENCE:
   kep_bbtemp(Ts, D, p)

 INPUTS:
   Ts -- stellar effective temperature (K)
   D -- eclipse depth
   p -- Rp/Rs radius ratio

 EXAMPLE:
   ;parameters from Jackson+ (2012), HAT-P-7 b -- Should be 2676 K
   print,kep_bbtemp(6350., 61d-6, 1./12.85)

 MODIFICATION HISTORY:
   2013 Jul 15 -- Written by Brian Jackson (decaelus@gmail.com)

(See kep_bbtemp.pro)


KEP_DOP_Q

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   kep_dop_q
     
 PURPOSE:
   Returns a planet-star mass ratio for a given Doppler beaming amplitude
   (See Loeb & Gaudi [2003] ApJ 588, 117.)

 CALLING SEQUENCE:
   kep_dop_q(Adop, Ts, Ms, P, sin_i)

 INPUTS:
   Adop -- Doppler beaming amplitude
   Ts -- stellar effective temperature (K)
   Ms -- stellar mass (solar masses)
   P -- orbital period (days)
   sin_i -- sine of the orbital inclination (probably close to 1)

 REQUIRES:
   alpha_beam() -- See http://www.lpl.arizona.edu/~bjackson/idl_code/index.html

 RESTRICTIONS:
   Assumes blackbody radiation from star

 EXAMPLE:
   Adop = 2.79d-6 & Ts = 5850 & Ms = 0.94 & P = 2.47 & sin_i = 0.994
   ;Should be 0.0011274410
   print, kep_dop_q(Adop, Ts, Ms, P, sin_i)

 MODIFICATION HISTORY:
   2013 Jul 15 -- Written by Brian Jackson (decaelus@gmail.com)

(See kep_dop_q.pro)


KEP_EV_Q

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   kep_ev_q
     
 PURPOSE:
   Returns the mass ratio q for observed ellipsoidal variations; 
   See Mazeh & Faigler (2010) A&A 521, L59.

 CALLING SEQUENCE:
   kep_ev_q(Aev, u, g, a, sin_i)

 INPUTS:
   Aev -- ellipsoidal variation amplitude
   u -- linear limb-darkening coefficient
   g -- gravity-darkening coefficient
   a -- semi-major axis (scaled to the stellar radius)
   sin_i -- sine of the orbital inclination (probably close to 1)

 EXAMPLE:
   Aev = 2.79d-6 & a = 7.931 & u = 0.579 & g = 4.*0.0884 & sin_i = 0.994
   ;Should be 0.00108
   print, kep_ev_q(Aev, u, g, a, sin_i)

 MODIFICATION HISTORY:
   2013 Jul 15 -- Written by Brian Jackson (decaelus@gmail.com)

(See kep_ev_q.pro)


LORENTZIAN

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   lorentzian
     
 PURPOSE:
   Returns a Lorentzian profile

 CALLING SEQUENCE:
   y = lorentzian(x, [x0, gamma])

 INPUTS:
   x -- time array
   p -- parameter array with p[0] = x0 and p[1] = gamma

 EXAMPLE:
   x = dindgen(1001)
   x0 = randomu(seed)*n_elements(x)
   gamma = 0.1
   plot, x, lorentzian(x, [x0, gamma])

 MODIFICATION HISTORY:
   2013 Jul 19 -- Written by Brian Jackson (decaelus@gmail.com)

(See lorentzian.pro)


MAD

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   mad
     
 PURPOSE:
   Returns the median absolute deviation

 CALLING SEQUENCE:
   mad = mad(data)

 INPUTS:
   data -- data

 EXAMPLE:
   data = randomn(seed, 101)
   print, mad(data)

 MODIFICATION HISTORY:
   2011 Nov 4 -- Written by Brian Jackson (decaelus@gmail.com)

(See mad.pro)


MH_GIBBS

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   mh_gibbs
     
 PURPOSE:
   Markov-chain Monte-Carlo routine using the Metropolis-Hastings algorithm
   and Gibbs sampler

 CALLING SEQUENCE:
   mh_gibbs, num_links, num_chains, model_func, like_func, x, y, err, init_p, $
     betamu, chain[, save_file=save_file, parinfo=parinfo, 
     restore_file=restore_file, /verbose]

 INPUTS:
   num_links -- number of links in each Markov chain; typically 1000s
   num_chains -- number of simultaneous Markov chains; typically ~5
   model_func -- the name of the model function to evaluate
   like_func -- the name of the function to evaluate log(likelihood) for jumping
   x/y/err -- in/dependent variables and associated uncertainties
   init_p -- initial guesses for model parameters
   betamu -- scaling params -- Ford [2005] ApJ 129, 1706 and find_betas.pro

 OUTPUTS:
   chain -- MCMC chain, shaped as [# fit params, num_chains, num_links], so,
     for example, chain[0,0,*] represents all values in the first chain for
     the first fit parameter.

 KEYWORD PARAMETERS:
   save_file -- name of save file to which to write chain
   verbose -- If set, the acceptance fraction for each chain will be printed
     every 100 links. If only one parameter is varied, acceptance fraction 
     should be ~0.44. If more than one, it should be ~0.25.
   parinfo -- structure containing information about parameters, exactly
     like parinfo for mpfitfun 
     (http://www.physics.wisc.edu/~craigm/idl/down/mpfitfun.pro), 
     except that only the .limited and .limits fields are used
   restore_file -- name of IDL save file that contains a previously-run chain

 RESTRICTIONS:
   This routine isn't too clever, so be sure you know what you're doing.   

 EXAMPLE:
   test_mh_gibb.pro provides a test case.

 MODIFICATION HISTORY:
   2011 Nov 27 -- Written by Brian Jackson (decaelus@gmail.com).
   2013 Apr 17 -- Modifed by BJ to allow for likelihood function
   2013 Jun 3 -- Added parinfo keyword
   2013 Jun 11 -- BJ changed the way in which initial values are chosen and
     in which variables are saved

(See mh_gibbs.pro)


MK_HTML_HELP2

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	MK_HTML_HELP2

 PURPOSE:
	Given a list of IDL procedure files (.PRO) or directories that
	contain such files, this procedure generates a file in the HTML
	format that contains the documentation for those routines that
	contain a DOC_LIBRARY style documentation template.  The output
	file is compatible with World Wide Web browsers.

 CATEGORY:
	Help, documentation.

 CALLING SEQUENCE:
	MK_HTML_HELP, Sources, Outfile

 INPUTS:
     Sources:  A string or string array containing the name(s) of the
		.pro files (or the names of directories containing such files)
		for which help is desired. If a source file is an IDL procedure,
		it must include the .PRO file extension. All other source files
		are assumed to be directories.

     Outfile:	The name of the output file which will be generated.

 KEYWORDS:
     TITLE:	If present, a string which supplies the name that
		should appear as the Document Title for the help.

     VERBOSE:	Normally, MK_HTML_HELP does its work silently.
		Setting this keyword to a non-zero value causes the procedure
		to issue informational messages that indicate what it
		is currently doing. !QUIET must be 0 for these messages
               to appear.

     STRICT: If this keyword is set to a non-zero value, MK_HTML_HELP
     will adhere strictly to the HTML format by scanning the document
     headers for characters that are reserved in HTML (<,>,&,").
     These are then converted to the appropriate HTML syntax in the
     output file. By default, this keyword is set to zero (to allow
     for faster processing).

 COMMON BLOCKS:
	None.

 SIDE EFFECTS:
	A help file with the name given by the Outfile argument is
	created.

 RESTRICTIONS:
	The following rules must be followed in formatting the .pro
	files that are to be searched.
		(a) The first line of the documentation block contains
		    only the characters ";+", starting in column 1.
     (b) There must be a line which contains the string "NAME:",
         which is immediately followed by a line containing the
         name of the procedure or function being described in
         that documentation block.  If this NAME field is not
         present, the name of the source file will be used.
		(c) The last line of the documentation block contains
		    only the characters ";-", starting in column 1.
		(d) Every other line in the documentation block contains
		    a ";" in column 1.

  Note that a single .pro file can contain multiple procedures and/or
  functions, each with their own documentation blocks. If it is desired
  to have "invisible" routines in a file, i.e. routines which are only
  for internal use and should not appear in the help file, simply leave
  out the ";+" and ";-" lines in the documentation block for those
  routines.

	No reformatting of the documentation is done.

 MODIFICATION HISTORY:
  July 5, 1995, DD, RSI. Original version.
  July 13, 1995, Mark Rivers, University of Chicago. Added support for
          multiple source directories and multiple documentation
          headers per .pro file.
  July 17, 1995, DD, RSI. Added code to alphabetize the subjects;
          At the end of each description block in the HTML file,
          added a reference to the source .pro file.
  July 18, 1995, DD, RSI. Added STRICT keyword to handle angle brackets.
  July 19, 1995, DD, RSI. Updated STRICT to handle & and ".
          Changed calling sequence to accept .pro filenames, .tlb
          text library names, and/or directory names.
          Added code to set default subject to name of file if NAME
          field is not present in the doc header.
  27 August 2003, AB, Remove VMS support. Fix bug with leaking
          file units.
  8 March 2005, DD, RSI. If none of the input .pro files contain
          documentation headers, display an informative error
          message, clean up, and exit. Also report the full path
          to the output file after creation.
  2013 Jul 16 -- Brian Jackson (decaelus@gmail.com) added linking to source 

(See mk_html_help2.pro)


MODIFIED_LORENTZIAN

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   modified_lorentzian
     
 PURPOSE:

 CALLING SEQUENCE:
   lor = modified_lorentzian(time, [t0, width, depth])

 INPUTS:
   time -- sampled times
   t0 -- center time for Lorentzian
   width -- full-width-half-max for Lorentzian profile
   depth -- depth of Lorentzian profile (should be positive)

 REQUIRES:
   lorentzian -- http://www.lpl.arizona.edu/~bjackson/idl_code/index.html

 EXAMPLE:
   time = dindgen(1001)
   t0 = randomu(seed)*(max(time) - min(time)) + min(time)
   width = 300.
   depth = 0.1
   lor = modified_lorentzian(time, [t0, width, depth])
   plot, time, lor

 MODIFICATION HISTORY:
   2013 Aug 4 -- Written by Brian Jackson (decaelus@gmail.com)

(See modified_lorentzian.pro)


ORB_POS

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	orb_pos

 PURPOSE:
	This routine returns the 3-D position vector of an orbital body, given the orbital 
	elements. Uses equations worked out in Murray & Dermott (1999) _Solar System Dynamics_,
	ch. 2. Available with no warranties whatsoever at
	http://www.lpl.arizona.edu/~bjackson/code/idl.html.

 CATEGORY:
	Astrophysics.

 CALLING SEQUENCE:

	Result = orb_pos(semi, ecc, asc_node, peri_long, inc, mean_anom)

 INPUTS:
	semi:		orbital semi-major axis
	ecc:		eccentricity
	asc_node:	longitude of ascending node, in degrees
	peri_long: 	longitude of pericenter, in degrees
	inc: 		orbital inclination, in degrees
	mean_anom: 	orbital mean anomaly

 OUTPUTS:
	3 x n_elements(mean_anom) array with x, y, z position of orbital body

 RESTRICTIONS:
	Code doesn't do any checking of parameters, so do them yourself. The code requires a 
	companion routine keplereq.pro (written by Marc Buie and Joern Wilms, among others) 
	available here: http://www.lpl.arizona.edu/~bjackson/code/idl.html.

 EXAMPLE:
	semi = 1.0
	ecc = 0.0 ;orbital eccentricity
	asc_node = 0. ;longitude of planetary ascending node
	peri_long = 0. ;longitude of planetary pericenter
	inc = 83.1 ;orbital inclination in degrees
	mean_anom = 2.*!pi*dindgen(101)
	r = orb_pos(semi, ecc, asc_node, peri_long, inc, mean_anom)

 MODIFICATION HISTORY:
 	Written by:	Brian Jackson (decaelus@gmail.com), 2011 July 25.

(See orb_pos.pro)


P3

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   p3
     
 PURPOSE:
   Makes a scatter plot with a couple default options I usually invoke

 CALLING SEQUENCE:
   p3, [x,] y, _extra=_extra

 INPUTS:
   x/y -- data points

 KEYWORD PARAMETERS:
   All keyword parameters are passed to PLOT.

 EXAMPLE:
   p3, findgen(10)

 MODIFICATION HISTORY:
   2011 Mar 8 -- WRitten by Brian jackson (decaelus@mgail.com)

(See p3.pro)


PL_PHS_CRV

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	pl_phs_crv

 PURPOSE:
	Returns a sinusoidal planetary phase curve for a planet in 
 	a circular orbit

 CATEGORY:
	Astrophysics

 CALLING SEQUENCE:

	pl_phs_crv = pl_phs_crv(phs, F0, F1)

 INPUTS:
	phs: 	orbital phase (0 = mid-transit, 0.5 = mid-eclipse)
	F0/F1:	phase curve = F0 - F1*cos(2*!pi*phs)
	
 OUTPUTS:
	planetary phase curve

 MODIFICATION HISTORY:
 	2012 Sep 21 -- Written by Brian Jackson (decaelus@gmail.com)

(See pl_phs_crv.pro)


PUSH

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   push
     
 PURPOSE:
   Adds an element to a vector

 CALLING SEQUENCE:
   push, array, value

 INPUTS:
   array -- the vector to which to add a new element
   value -- the new element

 RESTRICTIONS:
   This routine isn't very smart, so don't try to do anything clever, or 
     you'll get what's coming to you.

 EXAMPLE:
   push, new_vector, [1.]

 MODIFICATION HISTORY:
   Taken from some forgotten website long ago by 
     Brian Jackson (decaelus@gmail.com)

(See push.pro)


RANDOM_NON_UNI

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
      random_non_uni
     
 PURPOSE:
      Returns a random number with a given probability distribution

 CALLING SEQUENCE:
      result = random_non_uni(pop, num_wanted)

 INPUTS:
	pop: the population of values that gives the probability distribution
	num_wanted: the number of random values to return

 RESTRICTIONS:
	Does NOT return multi-dimensional arrays, like randomu and randomn.
	Also, this routine doesn't do any checking of the input values, so make
	  sure you know what you're doing.

 EXAMPLE:
 ;Some arbitrary, random population of numbers
 pop = exp(-randomu(seed, 101))
 res = random_non_uni(pop, 10000)

 MODIFICATION HISTORY:
	Written by Brian Jackson (decaelus@gmail.com), 2012 June 19

(See random_non_uni.pro)


ROBUST_STDDEV

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   robust_stddev
     
 PURPOSE:
   Returns outlier-resistant standard deviation, 1.4826 x MAD

 CALLING SEQUENCE:
   std = robust_stddev(data)

 INPUTS:
   data -- data

 KEYWORD PARAMETERS:
   iterate -- iteratively removes outliers to calculate stddev

 REQUIRES:
   mad -- http://www.lpl.arizona.edu/~bjackson/idl_code/index.html#MAD

 EXAMPLE:
   data = randomn(seed, 1001)
   print, robust_stddev(data)

 MODIFICATION HISTORY:
   2013 Jun 28 -- Written by Brian Jackson (decaelus@gmail.com)

(See robust_stddev.pro)


SIN_F

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	sin_f

 PURPOSE:
	This function returns the sine of the true anomaly, using relations given in 
	Murray & Dermott (1999), ch. 2.4. This routine uses keplereq.pro.

 CATEGORY:
	Celestial Mechanics.

 CALLING SEQUENCE:
	Result = sin_f(mean_anom, ecc)

 INPUTS:
	mean_anom: orbital mean anomaly
	ecc: orbital eccentricity

 OUTPUTS:
	Sine of the orbital true anomaly	

 EXAMPLE:
       ;make mean anomaly array
	num = 101
	mean_anom = 2.*!pi*dindgen(num)
	;orbital eccentricity of 0.1
	ecc = 0.1
	sin_f = sin_f(mean_anom, ecc)

 MODIFICATION HISTORY:
 	Written by:	Brian Jackson (decaelus@gmail.com), 2011 Jul 25.

(See sin_f.pro)


TEST_KEPCOTREND

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   test_kepcotrend
     
 PURPOSE:
   Provides an example for how to use kepcotrend.pro

 CALLING SEQUENCE:
   test_kepcotrend

 EXAMPLE:
   test_kepcotrend

 MODIFICATION HISTORY:
   2012 Dec 13 -- Written by Brian Jackson (decaelus@gmail.com)

(See test_kepcotrend.pro)


TEST_MH_GIBBS

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   test_mh_gibbs
     
 PURPOSE:
   Provides an example showing how to use find_betas and mh_gibbs

 CALLING SEQUENCE:
   test_mh_gibbs

 INPUTS:

 REQUIRES:
   mh_gibbs, find_betas, log_like_chisq -- 
     available at http://www.lpl.arizona.edu/~bjackson/idl_code/index.html

 SIDE EFFECTS:
   Plots the Markov Chain to the screen and shows the correct value

 EXAMPLE:
   test_mh_gibbs ;That's it.

 MODIFICATION HISTORY:
   2013 Jun 10 -- Written by Brian Jackson (decaelus@gmail.com)
   2014 Jan 12 -- Added complete comments

(See test_mh_gibbs.pro)


TRIM_ZEROS

[Previous Routine] [List of Routines]
 Name        : 
	TRIM_ZEROS
 Purpose     : 
	Converts numbers to strings, without trailing zeros.
 Explanation : 
	Converts numbers into a string representation, and trims off leading
	and/or trailing blanks.  Differs from STRTRIM in that trailing zeros
	after the period are also trimmed off, unless NUMBER is already a
	string, or an explicit format is passed.
 Use         : 
	Result = TRIM_ZEROS( NUMBER  [, FORMAT ]  [, FLAG ] )
 Inputs      : 
	NUMBER	= Variable or constant.  May be of any ordinary including
		  string.  However, structures are not allowed.
 Opt. Inputs : 
	FORMAT	- Format specification for STRING function.  Must be a string
		  variable, start with the "(" character, end with the ")"
		  character, and be a valid FORTRAN format specification.  If
		  NUMBER is complex, then FORMAT will be applied separately to
		  the real and imaginary parts.

	FLAG	- Flag passed to STRTRIM to control the type of trimming:

			FLAG = 0	Trim trailing blanks.
			FLAG = 1	Trim leading blanks.
			FLAG = 2	Trim both leading and trailing blanks.

		  The default value is 2.  If NUMBER is complex, then FORMAT
		  will be applied separately to the real and imaginary parts.

 Outputs     : 
	Function returns as a string variable representing the value NUMBER.
 Opt. Outputs: 
	None.
 Keywords    : 
	None.
 Calls       : 
	None.
 Common      : 
	None.
 Restrictions: 
	NUMBER must not be a structure.
	FORMAT must be a valid format specification, and must not be passed
		if NUMBER is of type string.
	FLAG must not be of string type, or an array.
 Side effects: 
	None.
 Category    : 
	Utilities, Strings.
 Prev. Hist. : 
	William Thompson	Applied Research Corporation
	May, 1987		8201 Corporate Drive
				Landover, MD  20785

	William Thompson, Feb. 1992, added support for complex numbers, and
				     fixed Unix problem with lowercase "e".
 Written     : 
	William Thompson, GSFC, May 1987.
 Modified    : 
	Version 1, William Thompson, GSFC, 9 April 1993.
		Incorporated into CDS library.
	Version 2, Zarro (SAC/GSFC), 3-Jun-98
		Added check for undefined input
       Version 3, Zarro (SM&A/GSFC), 1-Dec-99
               Returned invalid input as blank string
               to avoid downstream problems.
       Version 4, Zarro (SM&A/GSFC), 4-Jan-00
               Added /QUIET
       Version 5, Zarro (SM&A/GSFC), 20-Jan-00
               Vectorized
	Version 6, 24-Jan-2000, William Thompson, GSFC
		Fixed bug introduced in version 5.
	Version 7, 14-Mar-2000, Zarro (SM&A/GSFC)
		Moved check for unsupported type ahead of recursion

(See trim_zeros.pro)