fit_sed_mpfit_hz

dust_emissivity.fit_sed.fit_sed_mpfit_hz(xdata, flux, guesses=(0, 0), err=None, blackbody_function='blackbody', quiet=True, sc=1e+20, **kwargs)

Parameters:

xdata : array

Array of the frequencies of the data

flux : array

The fluxes corresponding to the xdata values. Should be in erg/s/cm^2/Hz

guesses : (Temperature,Column) or (Temperature,Beta,Column)

The input guesses. 3 parameters are used for modified blackbody fitting, two for temperature fitting.

blackbody_function: str

The blackbody function to fit, either ‘blackbody’, ‘modified’, or ‘modified_blackbody’

quiet : bool

quiet flag passed to mpfit

sc : float

A numerical parameter to enable the fitter to function properly. It is unclear what values this needs to take, 1e20 seems to work by bringing the units from erg/s/cm^2/Hz to Jy, i.e. bringing them into the “of order 1” regime. This does NOT affect the output units, though it may affect the quality of the fit.

Returns:

mp : mpfit structure

An mpfit structure. Access parameters and errors via mp.params and mp.perror. The covariance matrix is in mp.covar.

Examples

>>> from astropy import units as u
>>> import numpy as np
>>> wavelengths = np.array([20,70,160,250,350,500,850,1100]) * u.um
>>> frequencies = wavelengths.to(u.Hz, u.spectral())
>>> temperature = 15 * u.K
>>> column = 1e22 * u.cm**-2
>>> flux = modified_blackbody(frequencies, temperature, beta=1.75,
...                           column=column)
>>> err = 0.1 * flux
>>> np.random.seed(0)
>>> noise = np.random.randn(frequencies.size) * err
>>> tguess, bguess, nguess = 20.,2.,21.5
>>> bbunit = u.erg/u.s/u.cm**2/u.Hz
>>> mp = fit_sed_mpfit_hz(frequencies.to(u.Hz).value,
...                       (flux+noise).to(bbunit).value, err=err.to(bbunit).value,
...                       blackbody_function='modified',
...                       guesses=(tguess, bguess, nguess))
>>> print(mp.params)
[ 14.99095224   1.78620237  22.05271119]
>>> # T~14.9 K, beta ~1.79, column ~10^22