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... hmm (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

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 TKIN fitter translated from Erik Rosolowsky’s http://svn.ok.ubc.ca/svn/signals/nh3fit/

Module API

pyspeckit.spectrum.models.ammonia.ammonia(xarr, tkin=20, tex=None, ntot=100000000000000.0, width=1, xoff_v=0.0, fortho=0.0, tau=None, fillingfraction=None, return_tau=False, background_tb=2.7315, thin=False, verbose=False, return_components=False, debug=False)

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

Parameters:

xarr: `pyspeckit.spectrum.units.SpectroscopicAxis` :

Array of wavelength/frequency values

ntot: float :

can be specified as a column density (e.g., 10^15) or a log-column-density (e.g., 15)

tex: float or None :

Excitation temperature. Assumed LTE if unspecified (None), if tex>tkin, or if thin is specified.

ntot: float :

Total column density of NH3. Can be specified as a float in the range 5-25 or an exponential (1e5-1e25)

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

thin: bool :

uses a different parametetrization and requires only the optical depth, width, offset, and tkin to be specified. In the ‘thin’ approximation, tex is not used in computation of the partition function - LTE is implicitly assumed

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)

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(x, A, dx, w, return_components=False, normalized=False, return_hyperfine_components=False)

Returns a 1-dimensional gaussian of form A*numpy.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:

x : np.ndarray

array of x values

A : 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

w : 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 rest of this needs to be translated from C

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_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

The simplest and most useful model.

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 at 0x7f1c7c780f50>)

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)

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

Model Documentation Table of Contents

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

[edit this page on github]

[edit this page on bitbucket]