Models

See Parameters for information on how to restrict/modify model parameters.

The generic SpectralModel class is a wrapper for model functions. A model should take in an X-axis and some number of parameters. In order to declare a SpectralModel, you give SpectralModel the function name and the number of parameters it requires. The rest of the options are optional, though parnames & shortvarnames are strongly recommended. If you do not specify fitunits, your fitting code must deal with units internally.

Here are some examples of how to make your own fitters:

hill5_fitter = model.SpectralModel(hill5_model, 5,
        parnames=['tau', 'v_lsr',  'v_infall',  'sigma', 'tpeak'],
        parlimited=[(True,False),(False,False),(True,False),(True,False), (True,False)],
        parlimits=[(0,0), (0,0), (0,0), (0,0), (0,0)],
        # specify the parameter names (TeX is OK)
        shortvarnames=("\\tau","v_{lsr}","v_{infall}","\\sigma","T_{peak}"),
        fitunits='Hz' )

gaussfitter = model.SpectralModel(gaussian, 3,
        parnames=['amplitude','shift','width'],
        parlimited=[(False,False),(False,False),(True,False)],
        parlimits=[(0,0), (0,0), (0,0)],
        shortvarnames=('A',r'\Delta x',r'\sigma'))

Then you can register these fitters.

API Documentation for Models

class pyspeckit.spectrum.models.model.SpectralModel(modelfunc, npars, shortvarnames=('A', '\Delta x', '\sigma'), fitunits=None, centroid_par=None, fwhm_func=None, fwhm_pars=None, integral_func=None, use_lmfit=False, **kwargs)

A wrapper class for a spectra model. Includes internal functions to generate multi-component models, annotations, integrals, and individual components. The declaration can be complex, since you should name individual variables, set limits on them, set the units the fit will be performed in, and set the annotations to be used. Check out some of the hyperfine codes (hcn, n2hp) for examples.

Spectral Model Initialization

Create a Spectral Model class for data fitting

Parameters:

modelfunc : function

the model function to be fitted. Should take an X-axis (spectroscopic axis) as an input followed by input parameters. Returns an array with the same shape as the input X-axis

npars : int

number of parameters required by the model

parnames : list (optional)

a list or tuple of the parameter names

parvalues : list (optional)

the initial guesses for the input parameters (defaults to ZEROS)

parlimits : list (optional)

the upper/lower limits for each variable (defaults to ZEROS)

parfixed : list (optional)

Can declare any variables to be fixed (defaults to ZEROS)

parerror : list (optional)

technically an output parameter. Specifying it here will have no effect. (defaults to ZEROS)

partied : list (optional)

not the past tense of party. Can declare, via text, that some parameters are tied to each other. Defaults to zeros like the others, but it’s not clear if that’s a sensible default

fitunits : str (optional)

convert X-axis to these units before passing to model

parsteps : list (optional)

minimum step size for each paremeter (defaults to ZEROS)

npeaks : list (optional)

default number of peaks to assume when fitting (can be overridden)

shortvarnames : list (optional)

TeX names of the variables to use when annotating

Returns:

A tuple containing (model best-fit parameters, the model, parameter

errors, chi^2 value)

analytic_centroids(centroidpar=None)

Return the analytic centroids of the model components

Parameters:

centroidpar : None or string

The name of the parameter in the fit that represents the centroid some models have default centroid parameters - these will be used if centroidpar is unspecified

Returns:

List of the centroid values (even if there’s only 1)

analytic_fwhm(parinfo=None)

Return the FWHMa of the model components if a fwhm_func has been defined

Done with incomprehensible list comprehensions instead of nested for loops... readability sacrificed for speed and simplicity. This is unpythonic.

analytic_integral(modelpars=None, npeaks=None, npars=None)

Placeholder for analyic integrals; these must be defined for individual models

annotations(shortvarnames=None, debug=False)

Return a list of TeX-formatted labels

The values and errors are formatted so that only the significant digits are displayed. Rounding is performed using the decimal package.

Parameters:

shortvarnames : list

A list of variable names (tex is allowed) to include in the annotations. Defaults to self.shortvarnames

Examples

>>> # Annotate a Gaussian
>>> sp.specfit.annotate(shortvarnames=['A','\Delta x','\sigma'])

component_integrals(xarr, dx=None)

Compute the integrals of each component

components(xarr, pars, **kwargs)

Return a numpy ndarray of shape [npeaks x modelshape] of the independent components of the fits

computed_centroid(xarr=None)

Return the computed centroid of the model

Parameters:

xarr : None or np.ndarray

The X coordinates of the model over which the centroid should be computed. If unspecified, the centroid will be in pixel units

fitter(xax, data, err=None, quiet=True, veryverbose=False, debug=False, parinfo=None, **kwargs)

Run the fitter using mpfit.

kwargs will be passed to _make_parinfo and mpfit.

Parameters:

xax : SpectroscopicAxis

The X-axis of the spectrum

data : ndarray

The data to fit

err : ndarray (optional)

The error on the data. If unspecified, will be uniform unity

parinfo : ParinfoList

The guesses, parameter limits, etc. See pyspeckit.spectrum.parinfo for details

quiet : bool

pass to mpfit. If False, will print out the parameter values for each iteration of the fitter

veryverbose : bool

print out a variety of mpfit output parameters

debug : bool

raise an exception (rather than a warning) if chi^2 is nan

get_emcee_ensemblesampler(xarr, data, error, nwalkers, **kwargs)

Get an emcee walker ensemble for the data & model

Parameters:

data : np.ndarray

error : np.ndarray

nwalkers : int

Number of walkers to use

Examples

>>> import pyspeckit
>>> x = pyspeckit.units.SpectroscopicAxis(np.linspace(-10,10,50), unit='km/s')
>>> e = np.random.randn(50)
>>> d = np.exp(-np.asarray(x)**2/2.)*5 + e
>>> sp = pyspeckit.Spectrum(data=d, xarr=x, error=np.ones(50)*e.std())
>>> sp.specfit(fittype='gaussian')
>>> nwalkers = sp.specfit.fitter.npars * 2
>>> emcee_ensemble = sp.specfit.fitter.get_emcee_ensemblesampler(sp.xarr, sp.data, sp.error, nwalkers)
>>> p0 = np.array([sp.specfit.parinfo.values] * nwalkers)
>>> p0 *= np.random.randn(*p0.shape) / 10. + 1.0
>>> pos,logprob,state = emcee_ensemble.run_mcmc(p0,100)

get_emcee_sampler(xarr, data, error, **kwargs)

Get an emcee walker for the data & model

Parameters:

xarr : pyspeckit.units.SpectroscopicAxis

data : np.ndarray

error : np.ndarray

Examples

>>> import pyspeckit
>>> x = pyspeckit.units.SpectroscopicAxis(np.linspace(-10,10,50), unit='km/s')
>>> e = np.random.randn(50)
>>> d = np.exp(-np.asarray(x)**2/2.)*5 + e
>>> sp = pyspeckit.Spectrum(data=d, xarr=x, error=np.ones(50)*e.std())
>>> sp.specfit(fittype='gaussian')
>>> emcee_sampler = sp.specfit.fitter.get_emcee_sampler(sp.xarr, sp.data, sp.error)
>>> p0 = sp.specfit.parinfo
>>> emcee_sampler.run_mcmc(p0,100)

get_pymc(xarr, data, error, use_fitted_values=False, inf=inf, use_adaptive=False, return_dict=False, **kwargs)

Create a pymc MCMC sampler. Defaults to ‘uninformative’ priors

Parameters:

data : np.ndarray

error : np.ndarray

use_fitted_values : bool

Each parameter with a measured error will have a prior defined by the Normal distribution with sigma = par.error and mean = par.value

Examples

>>> x = pyspeckit.units.SpectroscopicAxis(np.linspace(-10,10,50), unit='km/s')
>>> e = np.random.randn(50)
>>> d = np.exp(-np.asarray(x)**2/2.)*5 + e
>>> sp = pyspeckit.Spectrum(data=d, xarr=x, error=np.ones(50)*e.std())
>>> sp.specfit(fittype='gaussian')
>>> MCuninformed = sp.specfit.fitter.get_pymc(sp.xarr, sp.data, sp.error)
>>> MCwithpriors = sp.specfit.fitter.get_pymc(sp.xarr, sp.data, sp.error, use_fitted_values=True)
>>> MCuninformed.sample(1000)
>>> MCuninformed.stats()['AMPLITUDE0']
>>> # WARNING: This will fail because width cannot be set <0, but it may randomly reach that...
>>> # How do you define a likelihood distribution with a lower limit?!
>>> MCwithpriors.sample(1000)
>>> MCwithpriors.stats()['AMPLITUDE0']

integral(modelpars, dx=None, **kwargs)

Extremely simple integrator: IGNORES modelpars; just sums self.model

lmfitfun(x, y, err=None, debug=False)

Wrapper function to compute the fit residuals in an lmfit-friendly format

lmfitter(xax, data, err=None, parinfo=None, quiet=True, debug=False, **kwargs)

Use lmfit instead of mpfit to do the fitting

Parameters:

xax : SpectroscopicAxis

The X-axis of the spectrum

data : ndarray

The data to fit

err : ndarray (optional)

The error on the data. If unspecified, will be uniform unity

parinfo : ParinfoList

The guesses, parameter limits, etc. See pyspeckit.spectrum.parinfo for details

quiet : bool

If false, print out some messages about the fitting

logp(xarr, data, error, pars=None)

Return the log probability of the model. If the parameter is out of range, return -inf

mpfitfun(x, y, err=None)

Wrapper function to compute the fit residuals in an mpfit-friendly format

n_modelfunc(pars=None, debug=False, **kwargs)

Simple wrapper to deal with N independent peaks for a given spectral model

slope(xinp)

Find the local slope of the model at location x (x must be in xax’s units)

Ammonia inversion transition TROT fitter translated from Erik Rosolowsky’s http://svn.ok.ubc.ca/svn/signals/nh3fit/

Module API

pyspeckit.spectrum.models.ammonia.ammonia(xarr, trot=20, tex=None, ntot=14, width=1, xoff_v=0.0, fortho=0.0, tau=None, fillingfraction=None, return_tau=False, background_tb=2.7315, verbose=False, return_components=False, debug=False, line_names=['oneone', 'twotwo', 'threethree', 'fourfour', 'fivefive', 'sixsix', 'sevenseven', 'eighteight'])

Generate a model Ammonia spectrum based on input temperatures, column, and gaussian parameters

Parameters:

xarr: `pyspeckit.spectrum.units.SpectroscopicAxis`

Array of wavelength/frequency values

trot: float

The rotational temperature of the lines. This is the excitation temperature that governs the relative populations of the rotational states.

tex: float or None

Excitation temperature. Assumed LTE if unspecified (None) or if tex>trot. This is the excitation temperature for all of the modeled lines, which means we are explicitly assuming T_ex is the same for all lines.

ntot: float

Total log column density of NH3. Can be specified as a float in the range 5-25

width: float

Line width in km/s

xoff_v: float

Line offset in km/s

fortho: float

Fraction of NH3 molecules in ortho state. Default assumes all para (fortho=0).

tau: None or float

If tau (optical depth in the 1-1 line) is specified, ntot is NOT fit but is set to a fixed value. The optical depths of the other lines are fixed relative to tau_oneone

fillingfraction: None or float

fillingfraction is an arbitrary scaling factor to apply to the model

return_tau: bool

Return a dictionary of the optical depths in each line instead of a synthetic spectrum

return_components: bool

Return a list of arrays, one for each hyperfine component, instead of just one array

background_tb : float

The background brightness temperature. Defaults to TCMB.

verbose: bool

More messages

debug: bool

For debugging.

Returns:

spectrum: numpy.ndarray

Synthetic spectrum with same shape as xarr

component_list: list

List of numpy.ndarray‘s, one for each hyperfine component

tau_dict: dict

Dictionary of optical depth values for the various lines (if return_tau is set)

pyspeckit.spectrum.models.ammonia.ammonia_thin(xarr, tkin=20, tex=None, ntot=14, width=1, xoff_v=0.0, fortho=0.0, tau=None, return_tau=False, **kwargs)

Use optical depth in the 1-1 line as a free parameter The optical depths of the other lines are then set by the kinetic temperature

tkin is used to compute trot assuming a 3-level system consisting of (1,1), (2,1), and (2,2) as in Swift et al, 2005 [2005ApJ...620..823S]

pyspeckit.spectrum.models.ammonia.cold_ammonia(xarr, tkin, **kwargs)

Generate a model Ammonia spectrum based on input temperatures, column, and gaussian parameters

Parameters:

xarr: `pyspeckit.spectrum.units.SpectroscopicAxis`

Array of wavelength/frequency values

tkin: float

The kinetic temperature of the lines in K. Will be converted to rotational temperature following the scheme of Swift et al 2005 (http://esoads.eso.org/abs/2005ApJ...620..823S, eqn A6) and further discussed in Equation 7 of Rosolowsky et al 2008 (http://adsabs.harvard.edu/abs/2008ApJS..175..509R)

Adds a variable height (background) component to any model

This is a formaldehyde 1_11-1_10 / 2_12-2_11 fitter. It includes hyperfine components of the formaldehyde lines and has both LTE and RADEX LVG based models

pyspeckit.spectrum.models.formaldehyde.formaldehyde(xarr, amp=1.0, xoff_v=0.0, width=1.0, return_hyperfine_components=False, texscale=0.01, tau=0.01, **kwargs)

Generate a model Formaldehyde spectrum based on simple gaussian parameters

the “amplitude” is an essentially arbitrary parameter; we therefore define it to be Tex given tau=0.01 when passing to the fitter The final spectrum is then rescaled to that value

pyspeckit.spectrum.models.formaldehyde.formaldehyde_pyradex(xarr, density=4, column=13, temperature=20, xoff_v=0.0, opr=1.0, width=1.0, tbackground=2.73, grid_vwidth=1.0, debug=False, verbose=False, **kwargs)

Use a grid of RADEX-computed models to make a model line spectrum

The RADEX models have to be available somewhere. OR they can be passed as arrays. If as arrays, the form should be: texgrid = ((minfreq1,maxfreq1,texgrid1),(minfreq2,maxfreq2,texgrid2))

xarr must be a SpectroscopicAxis instance xoff_v, width are both in km/s

grid_vwidth is the velocity assumed when computing the grid in km/s this is important because tau = modeltau / width (see, e.g., Draine 2011 textbook pgs 219-230)

pyspeckit.spectrum.models.formaldehyde.formaldehyde_radex(xarr, density=4, column=13, xoff_v=0.0, width=1.0, grid_vwidth=1.0, grid_vwidth_scale=False, texgrid=None, taugrid=None, hdr=None, path_to_texgrid='', path_to_taugrid='', temperature_gridnumber=3, debug=False, verbose=False, **kwargs)

Use a grid of RADEX-computed models to make a model line spectrum

The RADEX models have to be available somewhere. OR they can be passed as arrays. If as arrays, the form should be: texgrid = ((minfreq1,maxfreq1,texgrid1),(minfreq2,maxfreq2,texgrid2))

xarr must be a SpectroscopicAxis instance xoff_v, width are both in km/s

grid_vwidth is the velocity assumed when computing the grid in km/s this is important because tau = modeltau / width (see, e.g., Draine 2011 textbook pgs 219-230) grid_vwidth_scale is True or False: False for LVG, True for Sphere

pyspeckit.spectrum.models.formaldehyde.formaldehyde_radex_orthopara_temp(xarr, density=4, column=13, orthopara=1.0, temperature=15.0, xoff_v=0.0, width=1.0, Tbackground1=2.73, Tbackground2=2.73, grid_vwidth=1.0, grid_vwidth_scale=False, texgrid=None, taugrid=None, hdr=None, path_to_texgrid='', path_to_taugrid='', debug=False, verbose=False, getpars=False, **kwargs)

Use a grid of RADEX-computed models to make a model line spectrum

The RADEX models have to be available somewhere. OR they can be passed as arrays. If as arrays, the form should be: texgrid = ((minfreq1,maxfreq1,texgrid1),(minfreq2,maxfreq2,texgrid2))

xarr must be a SpectroscopicAxis instance xoff_v, width are both in km/s

grid_vwidth is the velocity assumed when computing the grid in km/s this is important because tau = modeltau / width (see, e.g., Draine 2011 textbook pgs 219-230) grid_vwidth_scale is True or False: False for LVG, True for Sphere

pyspeckit.spectrum.models.formaldehyde.formaldehyde_radex_tau(xarr, density=4, column=13, xoff_v=0.0, width=1.0, grid_vwidth=1.0, grid_vwidth_scale=False, taugrid=None, hdr=None, path_to_taugrid='', temperature_gridnumber=3, debug=False, verbose=False, return_hyperfine_components=False, **kwargs)

Use a grid of RADEX-computed models to make a model line spectrum

• uses hyperfine components
• assumes tau varies but tex does not!

The RADEX models have to be available somewhere. OR they can be passed as arrays. If as arrays, the form should be: texgrid = ((minfreq1,maxfreq1,texgrid1),(minfreq2,maxfreq2,texgrid2))

xarr must be a SpectroscopicAxis instance xoff_v, width are both in km/s

grid_vwidth is the velocity assumed when computing the grid in km/s this is important because tau = modeltau / width (see, e.g., Draine 2011 textbook pgs 219-230) grid_vwidth_scale is True or False: False for LVG, True for Sphere

The simplest and most useful model.

Until 12/23/2011, gaussian fitting used the complicated and somewhat bloated gaussfitter.py code. Now, this is a great example of how to make your own model! Just make a function like gaussian and plug it into the SpectralModel class.

pyspeckit.spectrum.models.inherited_gaussfitter.gaussian(xarr, amplitude, dx, width, return_components=False, normalized=False, return_hyperfine_components=False)

Returns a 1-dimensional gaussian of form A*np.exp(-(x-dx)**2/(2*w**2))

Area is sqrt(2*pi*sigma^2)*amplitude - i.e., this is NOT a normalized gaussian, unless normalized=True in which case A = Area

Parameters:

xarr : np.ndarray

array of x values

amplitude : float

Amplitude of the Gaussian, i.e. its peak value, unless normalized=True then A is the area of the gaussian

dx : float

Center or “shift” of the gaussian

width : float

Width of the gaussian (sigma)

return_components : bool

dummy variable; return_components does nothing but is required by all fitters

return_hyperfine_components : bool

dummy variable; does nothing but is required by all fitters

normalized : bool

Return a normalized Gaussian?

pyspeckit.spectrum.models.inherited_gaussfitter.gaussian_fitter()

Generator for Gaussian fitter class

pyspeckit.spectrum.models.inherited_gaussfitter.gaussian_integral(amplitude, sigma)

Integral of a Gaussian

pyspeckit.spectrum.models.inherited_gaussfitter.gaussian_vheight_fitter()

Generator for Gaussian fitter class

This is an HCN fitter... ref for line params: http://www.strw.leidenuniv.nl/~moldata/datafiles/hcn@hfs.dat

pyspeckit.spectrum.models.hcn.aval_dict = {'10-01': 2.4075e-05, '12-01': 2.4075e-05, '11-01': 2.4075e-05}

Line strengths of the 15 hyperfine components in J = 1 - 0 transition. The thickness of the lines indicates their relative weight compared to the others. Line strengths are normalized in such a way that summing over all initial J = 1 levels gives the degeneracy of the J = 0 levels, i.e., for JF1F = 012, three for JF1F = 011, and one for JF1F = 010. Thus, the sum over all 15 transitions gives the total spin degeneracy

pyspeckit.spectrum.models.hcn.hcn_radex(xarr, density=4, column=13, xoff_v=0.0, width=1.0, grid_vwidth=1.0, grid_vwidth_scale=False, texgrid=None, taugrid=None, hdr=None, path_to_texgrid='', path_to_taugrid='', temperature_gridnumber=3, debug=False, verbose=False, **kwargs)

Use a grid of RADEX-computed models to make a model line spectrum

The RADEX models have to be available somewhere. OR they can be passed as arrays. If as arrays, the form should be: texgrid = ((minfreq1,maxfreq1,texgrid1),(minfreq2,maxfreq2,texgrid2))

xarr must be a SpectroscopicAxis instance xoff_v, width are both in km/s

grid_vwidth is the velocity assumed when computing the grid in km/s this is important because tau = modeltau / width (see, e.g., Draine 2011 textbook pgs 219-230) grid_vwidth_scale is True or False: False for LVG, True for Sphere

Code translated from: https://bitbucket.org/devries/analytic_infall/overview

Original source: http://adsabs.harvard.edu/abs/2005ApJ...620..800D

pyspeckit.spectrum.models.hill5infall.hill5_model(xarr, tau, v_lsr, v_infall, sigma, tpeak, TBG=2.73)

The hill5 from de Vries and Myers 2005. This model implicitly has zero optical depth in the envelope, no envelope velocity, and has a fixed background radiation temperature (see Table 2 in the paper).

Parameters:

xarr : np.ndarray

array of x values

tau : float

tau_c, the core-center optical depth

v_lsr : float

The centroid velocity in km/s

v_infall : float

The infall velocity

sigma : float

The line width

tpeak : float

The peak brightness temperature

TBG : float

The background temperature

pyspeckit.spectrum.models.hill5infall.jfunc(t, nu)

t- kelvin nu - Hz?

class pyspeckit.spectrum.models.hyperfine.hyperfinemodel(line_names, voff_lines_dict, freq_dict, line_strength_dict, relative_strength_total_degeneracy)

Wrapper for the hyperfine model class. Specify the offsets and relative strengths when initializing, then you’ve got yourself a hyperfine modeler.

There are a wide variety of different fitter attributes, each designed to free a different subset of the parameters. Their purposes should be evident from their names.

Initialize the various parameters defining the hyperfine transitions

Parameters:

line_names: list

list of the line names to be used as indices for the dictionaries

voff_lines_dict: dict

a linename:v_off dictionary of velocity offsets for the hyperfine components. Technically, this is redundant with freq_dict

freq_dict: dict

frequencies of the indvidual transitions

line_strength_dict: dict

Relative strengths of the hyperfine components, usually determined by their degeneracy and Einstein A coefficients

hyperfine(xarr, Tex=5.0, tau=0.1, xoff_v=0.0, width=1.0, return_hyperfine_components=False, Tbackground=2.73, amp=None, return_tau=False, tau_total=None, vary_hyperfine_tau=False, vary_hyperfine_width=False)

Generate a model spectrum given an excitation temperature, optical depth, offset velocity, and velocity width.

Parameters:

return_tau : bool

If specified, return just the tau spectrum, ignoring Tex

tau_total : bool

If specified, use this instead of tau, and it tries to normalize to the peak of the line

vary_hyperfine_tau : bool

If set to true, allows the hyperfine transition amplitudes to vary and does not use the line_strength_dict. If set, tau must be a dict

hyperfine_addbackground(xarr, Tbackground=2.73, Tex=5.0, tau=0.1, xoff_v=0.0, width=1.0, return_tau=False, **kwargs)

Identical to hyperfine, but adds Tbackground as a constant continuum level

hyperfine_amp(xarr, amp=None, xoff_v=0.0, width=1.0, return_hyperfine_components=False, Tbackground=2.73, Tex=5.0, tau=0.1)

wrapper of self.hyperfine with order of arguments changed

hyperfine_background(xarr, Tbackground=2.73, Tex=5.0, tau=0.1, xoff_v=0.0, width=1.0, return_tau=False, **kwargs)

Identical to hyperfine, but with Tbackground free. Assumes already background-subtracted

hyperfine_tau(xarr, tau, xoff_v, width, **kwargs)

same as hyperfine, but with arguments in a different order, AND tau is returned instead of exp(-tau)

hyperfine_tau_total(xarr, tau_total, xoff_v, width, **kwargs)

same as hyperfine, but with arguments in a different order, AND tau is returned instead of exp(-tau), AND the peak tau is used

hyperfine_varyhf(xarr, Tex, xoff_v, width, *args, **kwargs)

Wrapper of hyperfine for using a variable number of peaks with specified tau

hyperfine_varyhf_amp(xarr, xoff_v, width, *args, **kwargs)

Wrapper of hyperfine for using a variable number of peaks with specified amplitude (rather than tau). Uses some opaque tricks: Tex is basically ignored, and return_tau means you’re actually returning the amplitude, which is just passed in as tau

hyperfine_varyhf_amp_width(xarr, xoff_v, *args, **kwargs)

Wrapper of hyperfine for using a variable number of peaks with specified amplitude (rather than tau). Uses some opaque tricks: Tex is basically ignored, and return_tau means you’re actually returning the amplitude, which is just passed in as tau

Until 12/23/2011, lorentzian fitting used the complicated and somewhat bloated gaussfitter.py code. Now, this is a great example of how to make your own model!

pyspeckit.spectrum.models.inherited_lorentzian.lorentzian(x, A, dx, w, return_components=False)

Returns a 1-dimensional lorentzian of form A/(2*pi)*w/((x-dx)**2 + ((w/2)**2))

[amplitude,center,width]

return_components does nothing but is required by all fitters

pyspeckit.spectrum.models.inherited_lorentzian.lorentzian_fitter()

Generator for lorentzian fitter class

Fit a line based on parameters output from a grid of models

pyspeckit.spectrum.models.modelgrid.gaussian_line(xax, maxamp, tau, offset, width)

A Gaussian line function in which the

pyspeckit.spectrum.models.modelgrid.line_model_2par(xax, center, width, gridval1, gridval2, griddim1, griddim2, maxampgrid, taugrid, linefunction=<function gaussian_line>)

Returns the spectral line that matches the given x-axis

xax, center, width must be in the same units!

pyspeckit.spectrum.models.modelgrid.line_params_2D(gridval1, gridval2, griddim1, griddim2, valuegrid)

Given a 2D grid of modeled line values - the amplitude, e.g. excitation temperature, and the optical depth, tau - return the model spectrum

griddims contains the names of the axes and their values... it should have the same number of entries as gridpars

N2H+ fitter

Reference for line params: Daniel, F., Dubernet, M.-L., Meuwly, M., Cernicharo, J., Pagani, L. 2005, MNRAS 363, 1083

http://www.strw.leidenuniv.nl/~moldata/N2H+.html

http://adsabs.harvard.edu/abs/2005MNRAS.363.1083D

Does not yet implement: http://adsabs.harvard.edu/abs/2010ApJ...716.1315K

pyspeckit.spectrum.models.n2hp.n2hp_radex(xarr, density=4, column=13, xoff_v=0.0, width=1.0, grid_vwidth=1.0, grid_vwidth_scale=False, texgrid=None, taugrid=None, hdr=None, path_to_texgrid='', path_to_taugrid='', temperature_gridnumber=3, debug=False, verbose=False, **kwargs)

Use a grid of RADEX-computed models to make a model line spectrum

The RADEX models have to be available somewhere. OR they can be passed as arrays. If as arrays, the form should be: texgrid = ((minfreq1,maxfreq1,texgrid1),(minfreq2,maxfreq2,texgrid2))

xarr must be a SpectroscopicAxis instance xoff_v, width are both in km/s

grid_vwidth is the velocity assumed when computing the grid in km/s this is important because tau = modeltau / width (see, e.g., Draine 2011 textbook pgs 219-230) grid_vwidth_scale is True or False: False for LVG, True for Sphere

pyspeckit.spectrum.models.n2hp.relative_strength_total_degeneracy = {'121-011': 9.0, '121-010': 9.0, '121-012': 9.0, '111-010': 9.0, '111-011': 9.0, '111-012': 9.0, '122-012': 9.0, '122-011': 9.0, '112-012': 9.0, '112-011': 9.0, '110-011': 9.0, '101-012': 9.0, '101-011': 9.0, '101-010': 9.0, '123-012': 9.0}

Line strengths of the 15 hyperfine components in J=1-0 transition. The thickness of the lines indicates their relative weight compared to the others. Line strengths are normalized in such a way that summing over all initial J = 1 levels gives the degeneracy of the J = 0 levels, i.e., for JF1F 012, three for JF1F 011, and one for JF1F 010. Thus, the sum over all 15 transitions gives the total spin degeneracy

pyspeckit.spectrum.models.inherited_voigtfitter.voigt(xarr, amp, xcen, sigma, gamma, normalized=False)

Normalized Voigt profile

z = (x+i*gam)/(sig*sqrt(2)) V(x,sig,gam) = Re(w(z))/(sig*sqrt(2*pi))

The area of V in this definition is 1. If normalized=False, then you can divide the integral of V by sigma*sqrt(2*pi) to get the area.

Original implementation converted from http://mail.scipy.org/pipermail/scipy-user/2011-January/028327.html (had an incorrect normalization and strange treatment of the input parameters)

Modified implementation taken from wikipedia, using the definition. http://en.wikipedia.org/wiki/Voigt_profile

Parameters:

xarr : np.ndarray

The X values over which to compute the Voigt profile

amp : float

Amplitude of the voigt profile if normalized = True, amp is the AREA

xcen : float

The X-offset of the profile

sigma : float

The width / sigma parameter of the Gaussian distribution

gamma : float

The width / shape parameter of the Lorentzian distribution

normalized : bool

Determines whether “amp” refers to the area or the peak of the voigt profile

pyspeckit.spectrum.models.inherited_voigtfitter.voigt_fitter()

Generator for voigt fitter class

pyspeckit.spectrum.models.inherited_voigtfitter.voigt_fwhm(sigma, gamma)

Approximation to the Voigt FWHM from wikipedia

http://en.wikipedia.org/wiki/Voigt_profile

Parameters:

sigma : float

The width / sigma parameter of the Gaussian distribution

gamma : float

The width / shape parameter of the Lorentzian distribution

pyspeckit.spectrum.models.inherited_voigtfitter.voigt_moments(self, *args, **kwargs)

Get the spectral moments from the moments package. Use the gaussian width for the lorentzian width (not a great guess!)

Hydrogen Models

Hydrogen in HII regions is typically assumed to follow Case B recombination theory.

The values for the Case B recombination coefficients are given by Hummer & Storey (1987). They are also computed in Hummer (1994) and tabulated at a wiki. I had to OCR and pull out by hand some of the coefficients.

pyspeckit.spectrum.models.hydrogen.add_to_registry(sp)

Add the Hydrogen model to the Spectrum’s fitter registry

pyspeckit.spectrum.models.hydrogen.find_lines(xarr)

Given a pyspeckit.units.SpectrosopicAxis instance, finds all the lines that are in bounds. Returns a list of line names.

pyspeckit.spectrum.models.hydrogen.hydrogen_fitter(sp, temperature=10000, tiedwidth=False)

Generate a set of parameters identifying the hydrogen lines in your spectrum. These come in groups of 3 assuming you’re fitting a gaussian to each. You can tie the widths or choose not to.

temperature [ 5000, 10000, 20000 ]

The case B coefficients are computed for 3 temperatures

tiedwidth [ bool ]

Should the widths be tied?

Returns a list of tied and guesses in the xarr’s units

pyspeckit.spectrum.models.hydrogen.hydrogen_model(xarr, amplitude=1.0, width=0.0, velocity=0.0, a_k=0.0, temperature=10000)

Generate a set of parameters identifying the hydrogen lines in your spectrum. These come in groups of 3 assuming you’re fitting a gaussian to each. You can tie the widths or choose not to.

Parameters:

sp : pyspeckit.Spectrum

The spectrum to fit

temperature : [ 5000, 10000, 20000 ]

The case B coefficients are computed for 3 temperatures

a_k : float

The K-band extinction normalized to 2.2 microns. Simple exponential.

width : float

Line width in km/s

velocity : float

Line center in km/s

amplitude : float

arbitrary amplitude of the first line (all other lines will be scaled to this value)

Returns:

np.ndarray with same shape as sp.xarr

pyspeckit.spectrum.models.hydrogen.rrl(n, dn=1, amu=1.007825, Z=1)

compute Radio Recomb Line freqs in GHz from Brown, Lockman & Knapp ARAA 1978 16 445

UPDATED:

Gordon & Sorochenko 2009, eqn A6

Parameters:

n : int

The number of the lower level of the recombination line (H1a is Lyman alpha, for example)

dn : int

The delta-N of the transition. alpha=1, beta=2, etc.

amu : float

The mass of the central atom

Z : int

The ionization parameter for the atom. Z=1 is neutral, Z=2 is singly ionized, etc. For hydrogen, only z=1 makes sense, since ionized hydrogen has no electrons and therefore cannot have recombination lines.

Returns:

frequency in GHz

Model Documentation Table of Contents

• Parameters
  • Simple Example
  • Tying Parameters
  • Making your own parinfo

[edit this page on github]

[edit this page on bitbucket]