#!/usr/bin/env """ make line datacubes and integrated maps and compare applying continuum selfcal NOTE: this is intended to be an interactive, iterative process so this is more a log that should be run by cutting and pasting into casa rather than as an executable script search "CHANGEME" for variables to be changed 10/9/15 MCA """ # ================================== Setup ===================================== # III_18 M5.5 F # 16:07:03.84 -39:11:11.8 field = 48 # CHANGEME fitspw = '2,3,4,7,8,9' # line-free channels for fitting continuum linespw = '0,1,3,5,6,8' # line spectral windows (C18O, 13CO, CN) robust = 0.5 # CHANGEME imsize = [640,640] cell = '0.03arcsec' imagermode = 'csclean' outframe = 'lsrk' veltype = 'radio' width = '1.0km/s' start = '0km/s' nchan = 8 xc = 330 # CHANGEME yc = 314 # CHANGEME in_a = 80 out_a = 120 aper = 1.25 boxwidth = 300. box = rg.box([xc-boxwidth,yc-boxwidth],[xc+boxwidth,yc+boxwidth]) # ================ Create continuum subtracted line datasets =================== uvcontsub(vis = 'f'+str(field)+'.vis', # full vis file for this field spw = linespw, # line spw (for cont subtraction) fitspw = fitspw, # cont spw combine = 'spw', solint = 'int', fitorder = 1, want_cont = False) # should not be changed. # ============================== 13CO line ===================================== # first try on 13CO line as thats going to be the brightest # os.system('rm -rf test_f'+str(field)+'_13co32*') clean(vis = 'f'+str(field)+'.vis.contsub', imagename = 'test_f'+str(field)+'_13co32', mode = 'velocity', start = start, width = width, nchan = nchan, outframe = outframe, veltype = veltype, restfreq = '330.58797GHz', niter = 2000, threshold = 0, interactive = True, imsize = imsize, cell = cell, weighting ='briggs', robust = robust, imagermode = imagermode) # use viewer to check channel maps and spectrum # make sure that velocity range is adequate and continuum subtraction ok imview(raster = [{'file':'test_f'+str(field)+'_13co32.image'}], contour = [{'file':'f'+str(field)+'_cont.fits'}]) # export cube to fits file fbase = 'test_f'+str(field)+'_13co32' os.system('rm -rf '+fbase+'.cube.fits') exportfits(imagename=fbase+'.image',fitsimage=fbase+'.cube.fits') # redo moment0 maps (now excluding noisy first channel!!) # os.system('rm -rf f'+str(field)+'_13co32.mom0*') immoments(imagename = 'f'+str(field)+'_13co32.image', # CHANGEME (based on above analysis!!) outfile = 'f'+str(field)+'_13co32.mom0', moments = [0], includepix = [-10.0,100.0], chans = ('range=[2km/s,5km/s]')) # CHANGEME (based on above analysis!!) # can see emission; best snr # export to fits file fbase = 'f'+str(field)+'_13co32.mom0' os.system('rm -rf '+fbase+'.fits') exportfits(imagename=fbase,fitsimage=fbase+'.fits') #### FIRST MOMENT MAP sigma = 15e-3 # Jy/beam in peak velocity channel os.system('rm -rf test_f48_13co32_mom1.image') immoments(imagename = 'test_f48_13co32.cube.fits', outfile = 'test_f48_13co32_mom1.image', moments = [1], includepix = [3.0*sigma,100.0], chans = ('range=[2km/s,5km/s]')) os.system('rm -rf test_f48_13co32_mom1.fits') exportfits(imagename='test_f48_13co32_mom1.image',fitsimage='test_f48_13co32_mom1.fits') os.system('rm -rf test_f48_13co32_mom1.image') # measure flux im_max = imstat(imagename = 'f'+str(field)+'_13co32.mom0')['max'][0] im_rms = imstat(imagename = 'f'+str(field)+'_13co32.mom0', region='annulus[['+str(xc)+'pix,'+str(yc)+'pix],['+str(in_a)+'pix,'+str(out_a)+'pix]]')['rms'][0] im_flux = imstat(imagename = 'f'+str(field)+'_13co32.mom0', region='circle[['+str(xc)+'pix,'+str(yc)+'pix],'+str(aper)+'arcsec]')['flux'][0] print 'Flux = {0:.2f} mJy, rms = {1:.2f} mJy, S/N = {2:.1f}'.format(1000*im_flux, 1000*im_rms, im_flux/im_rms) # Flux = 1236.34 mJy, rms = 27.23 mJy, S/N = 45.4 # view continuum and gas imview(raster=[{'file':'f'+str(field)+'_13co32.mom0'}], contour = [{'file':'f'+str(field)+'_cont.fits'}]) # re-center image on source and use measure.py to get COG flux # os.system('rm -rf f'+str(field)+'_13co32.mom0_cropped*') ia.fromimage(outfile = 'f'+str(field)+'_13co32.mom0_cropped.image', infile = 'f'+str(field)+'_13co32.mom0.fits', region = box ) ia.close() exportfits(imagename = 'f'+str(field)+'_13co32.mom0_cropped.image', fitsimage = 'f'+str(field)+'_13co32.mom0_cropped.fits') ''' Measuring COG for M/f48_13co32.mom0_cropped.fits Assuming object center (300.0,300.0) Background: -0.68 mJy/beam km/s RMS in annulus 4.0-6.0 arcsec = 28.19 mJy/beam km/s i radius flux err snr (asec) (mJy) (mJy) 0 0.10 41.41 7.59 5.5 1 0.20 165.46 16.90 9.8 2 0.30 385.02 29.36 13.1 3 0.40 615.18 39.80 15.5 4 0.50 822.80 50.40 16.3 5 0.60 965.70 59.91 16.1 6 0.70 1067.25 74.55 14.3 7 0.80 1143.76 81.16 14.1 8 0.90 1200.00 96.39 12.4 9 1.00 1242.11 96.15 12.9 10 1.10 1248.45 100.59 12.4 11 1.20 1238.96 101.91 12.2 F = 1248.45 mJy E = 100.59 mJy S = 12.41 D = 2.20 arcsec ''' # ======================== C18O line ========================================== # don't bother with selfcal as it doesn't help... # os.system('rm -rf f'+str(field)+'_c18o32*') clean(vis = 'f'+str(field)+'.vis.contsub', imagename = 'f'+str(field)+'_c18o32', mode = 'velocity', start = start, width = width, nchan = nchan, outframe = outframe, veltype = veltype, restfreq = '329.33055GHz', niter = 500, threshold = 0, interactive = False, mask = 'f48_cont_mask.crtf', imsize = imsize, cell = cell, weighting ='briggs', robust = robust, imagermode = imagermode) # couldn't see line, just cleaned in automated mode # in continuum region for all channels # use viewer to check channel maps and spectrum # make sure that velocity range is adequate and continuum subtraction ok imview(raster = [{'file':'f'+str(field)+'_c18o32.image'}], contour = [{'file':'f'+str(field)+'_cont.fits'}]) # export cube to fits file fbase = 'f'+str(field)+'_c18o32' os.system('rm -rf '+fbase+'.cube.fits') exportfits(imagename=fbase+'.image',fitsimage=fbase+'.cube.fits') # redo moment0 maps (now excluding noisy first channel!!) # os.system('rm -rf f'+str(field)+'_c18o32.mom0*') immoments(imagename = 'f'+str(field)+'_c18o32_b4sc.image', # CHANGEME (based on above analysis!!) outfile = 'f'+str(field)+'_c18o32.mom0', moments = [0], includepix = [-10.0,100.0], chans = ('range=[2km/s,5km/s]')) # can't see emission, using same velocity range as 13CO # export to fits file fbase = 'f'+str(field)+'_c18o32.mom0' os.system('rm -rf '+fbase+'.fits') exportfits(imagename=fbase,fitsimage=fbase+'.fits') # measure flux im_rms = imstat(imagename = 'f'+str(field)+'_c18o32.mom0', region='annulus[['+str(xc)+'pix,'+str(yc)+'pix],['+str(in_a)+'pix,'+str(out_a)+'pix]]')['rms'][0] im_flux = imstat(imagename = 'f'+str(field)+'_c18o32.mom0', region='circle[['+str(xc)+'pix,'+str(yc)+'pix],'+str(aper)+'arcsec]')['flux'][0] print 'Flux = {0:.2f} mJy, rms = {1:.2f} mJy, S/N = {2:.1f}'.format(1000*im_flux, 1000*im_rms, im_flux/im_rms) # Flux = 328.08 mJy, rms = 33.24 mJy, S/N = 9.9 # view continuum and gas imview(raster=[{'file':'f'+str(field)+'_c18o32.mom0'}], contour = [{'file':'f'+str(field)+'_cont.fits'}]) # re-center image on source and use measure.py to get COG flux # os.system('rm -rf f'+str(field)+'_c18o32.mom0_cropped*') ia.fromimage(outfile = 'f'+str(field)+'_c18o32.mom0_cropped.image', infile = 'f'+str(field)+'_c18o32.mom0.fits', region = box ) ia.close() exportfits(imagename = 'f'+str(field)+'_c18o32.mom0_cropped.image', fitsimage = 'f'+str(field)+'_c18o32.mom0_cropped.fits') ''' Measuring COG for M/f48_c18o32.mom0_cropped.fits Assuming object center (300.0,300.0) Background: -0.48 mJy/beam km/s RMS in annulus 4.0-6.0 arcsec = 33.79 mJy/beam km/s i radius flux err snr (asec) (mJy) (mJy) 0 0.10 36.26 8.18 4.4 1 0.20 109.28 21.35 5.1 2 0.30 175.22 37.40 4.7 3 0.40 203.27 48.34 4.2 4 0.50 229.70 69.57 3.3 5 0.60 269.48 63.85 4.2 6 0.70 294.99 69.35 3.1 <---where levels off 7 0.80 304.52 92.53 3.3 8 0.90 314.45 90.52 3.5 9 1.00 347.84 99.62 3.5 10 1.10 388.06 134.15 2.9 11 1.20 362.46 135.78 2.7 F = 388.06 mJy E = 134.15 mJy S = 2.89 D = 2.20 arcsec ''' # ======================== Image CN ================== # don't bother with selfcal as it doesn't help... # os.system('rm -rf f'+str(field)+'_cn32_b4sc*') clean(vis = 'f'+str(field)+'.vis.contsub', imagename = 'f'+str(field)+'_cn32_b4sc', mode = 'velocity', start = start, width = width, nchan = nchan, outframe = outframe, veltype = veltype, restfreq = '340.24777GHz', niter = 2000, threshold = 0, interactive = True, imsize = imsize, cell = cell, weighting ='briggs', robust = robust, imagermode = imagermode) # use viewer to check channel maps and spectrum # make sure that velocity range is adequate and continuum subtraction ok imview(raster = [{'file':'f'+str(field)+'_cn32_b4sc.image'}], contour = [{'file':'f'+str(field)+'_cont.fits'}]) # redo moment0 maps (now excluding noisy first channel!!) # os.system('rm -rf f'+str(field)+'_cn32.mom0*') immoments(imagename = 'f'+str(field)+'_cn32_b4sc.image', # CHANGEME (based on above analysis!!) outfile = 'f'+str(field)+'_cn32.mom0', moments = [0], includepix = [-10.0,100.0], chans = ('range=[2km/s,4km/s]')) # CHANGEME (based on above analysis!!) # can see emission; highest snr # export to fits file fbase = 'f'+str(field)+'_cn32.mom0' os.system('rm -rf '+fbase+'.fits') exportfits(imagename=fbase,fitsimage=fbase+'.fits') # measure flux im_rms = imstat(imagename = 'f'+str(field)+'_cn32.mom0', region='annulus[['+str(xc)+'pix,'+str(yc)+'pix],['+str(in_a)+'pix,'+str(out_a)+'pix]]')['rms'][0] im_flux = imstat(imagename = 'f'+str(field)+'_cn32.mom0', region='circle[['+str(xc)+'pix,'+str(yc)+'pix],'+str(aper)+'arcsec]')['flux'][0] print 'Flux = {0:.2f} mJy, rms = {1:.2f} mJy, S/N = {2:.1f}'.format(1000*im_flux, 1000*im_rms, im_flux/im_rms) # Flux = 564.58 mJy, rms = 19.91 mJy, S/N = 28.4 # view continuum and gas imview(raster=[{'file':'f'+str(field)+'_cn32.mom0'}], contour = [{'file':'f'+str(field)+'_cont.fits'}]) # re-center image on source and use measure.py to get COG flux # os.system('rm -rf f'+str(field)+'_cn32.mom0_cropped*') ia.fromimage(outfile = 'f'+str(field)+'_cn32.mom0_cropped.image', infile = 'f'+str(field)+'_cn32.mom0.fits', region = box ) ia.close() exportfits(imagename = 'f'+str(field)+'_cn32.mom0_cropped.image', fitsimage = 'f'+str(field)+'_cn32.mom0_cropped.fits') ''' diameter total rms snr 2.0 608.13 84.069 7.2 diameter total rms snr 0.8 274.70 29.267 9.4 diameter peak rms snr 0.6 102.59 20.135 5.1 '''