single_component.py 文件源码

def deblend_cube(region='OrionA',vmin=0.,vmax=20.):
  tex = fits.getdata('{0}/parameterMaps/{0}_Tex_DR1_rebase3_flag.fits'.format(region))
  mom0 = fits.getdata('{0}/{0}_NH3_11_DR1_rebase3_mom0_QA_trim.fits'.format(region))
  vlsr = fits.getdata('{0}/parameterMaps/{0}_Vlsr_DR1_rebase3_flag.fits'.format(region))
  sigv = fits.getdata('{0}/parameterMaps/{0}_Sigma_DR1_rebase3_flag.fits'.format(region))
  nnh3 = fits.getdata('{0}/parameterMaps/{0}_N_NH3_DR1_rebase3_flag.fits'.format(region))
  cube = SpectralCube.read('{0}/{0}_NH3_11_DR1_rebase3_trim.fits'.format(region))
  cube  = cube.with_spectral_unit(u.km/u.s,velocity_convention='radio')

  tpeak = mom0/(np.sqrt(2*np.pi)*sigv)

  vlsr[vlsr==0]=np.nan
  sigv[sigv==0]=np.nan

  deblend = np.zeros(cube.shape)
  hdr = cube.wcs.to_header()
  spaxis = cube.spectral_axis.value

  for plane in np.arange(deblend.shape[0]):
      deblend[plane,:,:] = tpeak*np.exp(-(spaxis[plane]-vlsr)**2/(2*sigv**2))

  newcube = SpectralCube(deblend,cube.wcs,header=hdr)
  slab = newcube.spectral_slab(vmin*u.km/u.s,vmax*u.km/u.s)
  slab.write('{0}/{0}_singlecomp.fits'.format(region),overwrite=True)