In [5]:
%matplotlib inline
import matplotlib 
#matplotlib.use('TkAgg')
import matplotlib.pyplot as plt

###Hill5 model in the pyspeckit package
import pyspeckit.spectrum.models.hill5infall as hill5infall
from pyspeckit.spectrum.classes import Spectrum
from pyspeckit.spectrum.units import SpectroscopicAxis
import numpy as np
from astropy import units as u

from astropy.io import fits

####an example of my own HNC line
hd = fits.getheader('myhnc.fits')
data = fits.getdata('myhnc.fits')

velo = hd['VELO']

dv = hd['deltav']

refc = hd['CRPIX1']

xaxis = np.divide(np.add(np.multiply(np.subtract(np.arange(1,hd['NAXIS1']+1), refc), dv), velo), 1000)
yaxis = data[0,0,0,:]
x =xaxis*u.km/u.s
xarr = SpectroscopicAxis(x, refX=hd['restfreq']*u.Hz, velocity_convention='radio')
sp = Spectrum(xarr=xarr, data=yaxis)
sp.plotter()
plt.xlim(5,12)
plt.ylabel('Main beam temperature (K)')
#plt.plot([8,8],[-0.5,2],'r')
Out[5]:
<matplotlib.text.Text at 0x1180267d0>
In [6]:
##register hill5 model for fit
sp.Registry.add_fitter('hill5_fitter', hill5infall.hill5_fitter, 5)

##fitting, 5 pars correspond to tau, v_lsr, v_in, delta_v, T_ex
sp.specfit(fittype='hill5_fitter', guesses = [8, 8.2, 0.03, 0.95/2.35, 7])

invalid value encountered in divide
invalid value encountered in divide
In [7]:
##fitting results
sp.specfit.parinfo
Out[7]:
[Param #0         TAU0 =      4.60516 +/-          2.4544   Range:   [0,inf),
 Param #1       V_LSR0 =      8.07115 +/-       0.0727654 ,
 Param #2    V_INFALL0 =     0.118602 +/-        0.103928   Range:   [0,inf),
 Param #3       SIGMA0 =     0.265835 +/-       0.0482812   Range:   [0,inf),
 Param #4       TPEAK0 =      5.98288 +/-         0.85056   Range:   [0,inf)]
In [8]:
sp.plotter()
sp.specfit.plot_fit()
plt.xlim(2,15)
plt.ylabel('Main beam temperature (K)')
Out[8]:
<matplotlib.text.Text at 0x1185b1490>
In [ ]: