Models

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'),
        multisingle=multisingle,
        )

Then you can register these fitters.

Generic SpectralModel wrapper

class pyspeckit.spectrum.models.model.SpectralModel(modelfunc, npars, shortvarnames=('A', '\Delta x', '\sigma'), multisingle='multi', **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.

modelfunc: the model function to be fitted. Should take an X-axis (spectroscopic axis) as an input, followed by input parameters. npars - number of parameters required by the model

parnames - a list or tuple of the parameter names

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

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

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

parerror - technically an output parameter... hmm (defaults to ZEROS)

partied - 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 - convert X-axis to these units before passing to model

parsteps - minimum step size for each paremeter (defaults to ZEROS)

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

shortvarnames - TeX names of the variables to use when annotating

multisingle - Are there multiple peaks (no background will be fit) or

just a single peak (a background may/will be fit)

annotations(shortvarnames=None)

Return a list of TeX-formatted labels

components(xarr, pars, **kwargs)

Return a numpy ndarray of the independent components of the fits

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

Run the fitter. Must pass the x-axis and data. Can include error, parameter guesses, and a number of verbosity parameters.

quiet - pass to mpfit. If False, will print out the parameter values

for each iteration of the fitter

veryverbose - print out a variety of mpfit output parameters

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

accepts tied, limits, limited, and fixed as keyword arguments.

They must be lists of length len(params)

parinfo - You can override the class parinfo dict with this, though

that largely defeats the point of having the wrapper class. This class does NO checking for whether the parinfo dict is valid.

kwargs are passed to mpfit after going through _make_parinfo to strip out things used by this class

integral(modelpars, **kwargs)

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

mpfitfun(x, y, err=None)

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

n_modelfunc(pars, **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

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, thin=False, verbose=False, return_components=False, debug=False)

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

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

tex can be specified or can be assumed LTE if unspecified, if tex>tkin, or if “thin”

is specified

“thin” 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

If tau is specified, ntot is NOT fit but is set to a fixed value fillingfraction is an arbitrary scaling factor to apply to the model fortho is the ortho/(ortho+para) fraction. The default is to assume all ortho. xoff_v is the velocity offset in km/s

tau refers to the optical depth of the 1-1 line. The optical depths of the other lines are fixed relative to tau_oneone

(not implemented) if tau is specified, ntot is ignored

SimpleFitter wrapper

Adds a variable height (background) component to any model

Formaldehyde cm-line fitter

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_components=False)

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

Gaussian Fitter

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)

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

[height,amplitude,center,width]

return_components does nothing but is required by all fitters

pyspeckit.spectrum.models.inherited_gaussfitter.gaussian_fitter(multisingle='multi')

Generator for Gaussian fitter class

HCN Hyperfine Fitter

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.4074999999999999e-05, '12-01': 2.4074999999999999e-05, '11-01': 2.4074999999999999e-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

Hill5 analytic infall model

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?

Generalized hyperfine component fitter

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

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

Initialize the various parameters defining the hyperfine transitions

line_names is a LIST of the line names to be used as indices for the dictionaries

voff_lines_dict is a linename:v_off dictionary of velocity offsets for the hyperfine components. Technically,

this is redundant with freq_dict

freq_dict - frequencies of the indvidual transitions

line_strength_dict - Relative strengths of the hyperfine components, usually determined by their degeneracy and

Einstein A coefficients

hyperfine(xarr, Tex=5.0, tau=0.10000000000000001, xoff_v=0.0, width=1.0, return_components=False, Tbackground=2.73, amp=None)

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

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

wrapper of self.hyperfine with order of arguments changed

Lorentzian Fitter

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(multisingle='multi')

Generator for lorentzian fitter class

Model Grid

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 0x106e77cf8>)

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

pyspeckit.spectrum.models.n2hp.line_strength_dict = {'121-011': 0.41699999999999998, '121-010': 0.55600000000000005, '121-012': 0.028000000000000001, '111-010': 0.33300000000000002, '111-011': 0.25, '111-012': 0.41699999999999998, '122-012': 0.41699999999999998, '122-011': 1.25, '112-012': 1.25, '112-011': 0.41699999999999998, '110-011': 0.33300000000000002, '101-012': 0.55000000000000004, '101-011': 0.33300000000000002, '101-010': 0.111, '123-012': 2.3300000000000001}

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

Voigt Profile Fitter

pyspeckit.spectrum.models.inherited_voigtfitter.voigt(xarr, amp, xcen, Gfwhm, Lfwhm)

voigt profile

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

Converted from http://mail.scipy.org/pipermail/scipy-user/2011-January/028327.html

pyspeckit.spectrum.models.inherited_voigtfitter.voigt_fitter(multisingle='multi')

Generator for voigt fitter class

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 [HummerStorey1987]. They are also computed in [Hummer1994] and tabulated at a [wiki]. I had to OCR and pull out by hand some of the coefficients.

[HummerStorey1987]

: Recombination-line intensities for hydrogenic ions. I - Case B calculations for H I and He II

[Hummer1994]

: http://adsabs.harvard.edu/abs/1994MNRAS.268..109H

[wiki]

: http://wiki.hmet.net/index.php/HII_Case_B_Recombination_Coefficients

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.rrl(n, dn=1, amu=1.007825)

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