#!/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 ===================================== # LupusIII_71 16:08:54.68 -39:37:43.9 # Class II, M1 field = 10 # 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 = '-1.0km/s' nchan = 11 xc = 322 # CHANGEME yc = 333 # 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) # 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') # 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)+'_13co32.image'}], contour = [{'file':'f'+str(field)+'_cont.fits'}]) # clear peaks # redo moment0 maps (now excluding noisy first channel!!) # os.system('rm -rf f'+str(field)+'_13co32.mom0*') immoments(imagename = 'test_f'+str(field)+'_13co32.image', outfile = 'test_f'+str(field)+'_13co32.mom0', moments = [0], includepix = [-10.0,100.0], chans = ('range=[1km/s,7km/s]')) # clear emission, gave best snr # export to fits file fbase = 'test_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_f10_13co32_mom1.image') immoments(imagename = 'test_f10_13co32.cube.fits', outfile = 'test_f10_13co32_mom1.image', moments = [1], includepix = [3.0*sigma,100.0], chans = ('range=[1km/s,7km/s]')) os.system('rm -rf test_f10_13co32_mom1.fits') exportfits(imagename='test_f10_13co32_mom1.image',fitsimage='test_f10_13co32_mom1.fits') os.system('rm -rf f10_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 = 2628.66 mJy, rms = 35.59 mJy, S/N = 73.9 # 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/f10_13co32.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 = 38.06 mJy/beam km/s i radius flux err snr (asec) (mJy) (mJy) 0 0.10 114.44 11.58 9.9 1 0.20 416.99 27.11 15.4 2 0.30 867.86 42.14 20.6 3 0.40 1268.69 50.60 25.1 4 0.50 1595.60 86.03 18.5 5 0.60 1831.60 112.38 16.3 6 0.70 2035.94 99.39 20.5 7 0.80 2194.36 113.84 19.3 8 0.90 2300.63 118.48 19.4 9 1.00 2380.78 146.14 16.3 <-- used thi 10 1.10 2468.25 161.28 15.3 11 1.20 2581.21 179.28 14.4 12 1.30 2682.09 222.64 12.0 13 1.40 2768.42 226.82 12.2 14 1.50 2847.59 198.59 14.3 F = 2847.59 mJy E = 198.59 mJy S = 14.34 D = 3.00 arcsec COG continues to grow with aperture size but change of slope around 2" diameter and in casa this looks like good aperture ''' # ======================== 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 = 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)+'_c18o32_b4sc.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=[1km/s,7km/s]')) # can see emission, same 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 = 794.41 mJy, rms = 44.92 mJy, S/N = 17.7 # 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/f10_c18o32.mom0_cropped.fits Assuming object center (300.0,300.0) Background: 0.25 mJy/beam km/s RMS in annulus 4.0-6.0 arcsec = 45.53 mJy/beam km/s i radius flux err snr (asec) (mJy) (mJy) 0 0.10 58.15 10.70 5.4 1 0.20 220.72 29.30 7.5 2 0.30 465.42 53.14 8.8 3 0.40 637.37 67.04 9.5 4 0.50 699.65 80.14 8.7 5 0.60 714.01 73.47 9.7 <--- used this 6 0.70 756.80 120.47 6.3 7 0.80 832.85 119.66 7.0 8 0.90 915.98 124.14 7.4 9 1.00 948.93 105.63 9.0 10 1.10 893.34 141.48 6.3 11 1.20 809.61 172.62 4.7 12 1.30 797.11 178.24 4.5 13 1.40 803.51 189.26 4.2 14 1.50 760.88 188.34 4.0 F = 948.93 mJy E = 105.63 mJy S = 8.98 D = 2.00 arcsec COG has bump after leveling off at 1.2" diameter which is casa looks like a good aperture size ''' # ======================== 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) # emission very extended, cleaned conservatively # 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=[1km/s,7km/s]')) # CHANGEME (based on above analysis!!) # clear emission, same as 13CO and C18O # 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 = 1834.72 mJy, rms = 31.98 mJy, S/N = 57.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 3.4 2158.85 354.942 6.1 diameter total rms snr 0.8 632.64 47.037 13.4 diameter peak rms snr 0.4 276.87 30.539 9.1 large aperture seems correct emission very extended '''