#!/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_103 # Class II, K6 field = 18 # CHANGEME fitspw = '2,3,4' # line-free channels for fitting continuum linespw = '0,1,3' # 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 = '-3km/s' nchan = 16 xc = 320 # CHANGEME yc = 310 # CHANGEME in_a = 80 out_a = 120 aper = 1.25 # CHANGEME 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 f'+str(field)+'_13co32*') clean(vis = 'f'+str(field)+'.vis.contsub', imagename = 'f'+str(field)+'_13co32', mode = 'velocity', start = start, width = width, nchan = nchan, outframe = outframe, veltype = veltype, restfreq = '330.58797GHz', niter = 500, threshold = 0, interactive = False, mask = 'f18_cont_mask.crtf', imsize = imsize, cell = cell, weighting ='briggs', robust = robust, imagermode = imagermode) # could not see line, so just cleaned in automated mode # 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'}]) # export cube to fits file fbase = '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_b4sc.image', # CHANGEME (based on above analysis!!) outfile = 'f'+str(field)+'_13co32.mom0', moments = [0], includepix = [-10.0,100.0], chans = ('range=[-1km/s,7km/s]')) # broad wings but this gives best SNR # export to fits file fbase = 'f'+str(field)+'_13co32.mom0' os.system('rm -rf '+fbase+'.fits') exportfits(imagename=fbase,fitsimage=fbase+'.fits') # measure flux # imview(raster=[{'file':'f'+str(field)+'_13co32.mom0'}]) 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 = 614.69 mJy, rms = 46.23 mJy, S/N = 13.3 # 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 G/f18_13co32.mom0_cropped.fits Assuming object center (300.0,300.0) Background: 0.18 mJy/beam km/s RMS in annulus 4.0-9.0 arcsec = 45.26 mJy/beam km/s i radius flux err snr (asec) (mJy) (mJy) 0 0.1 36.7 11.5 3.2 1 0.2 124.9 31.0 4.0 2 0.3 261.7 63.8 4.1 3 0.4 407.1 62.6 6.5 4 0.5 537.0 89.3 6.0 5 0.6 607.6 108.1 5.6 6 0.7 648.8 99.5 6.5 7 0.8 668.9 119.6 5.6 8 0.9 681.2 164.4 4.1 9 1.0 694.6 133.0 5.2 10 1.1 682.6 164.7 4.1 11 1.2 636.7 194.0 3.3 12 1.3 607.9 204.0 3.0 13 1.4 622.9 242.1 2.6 14 1.5 619.3 233.6 2.7 F = 694.6 mJy E = 133.0 mJy S = 5.2 D = 2.0 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 = 'f18_cont_mask.crtf', imsize = imsize, cell = cell, weighting ='briggs', robust = robust, imagermode = imagermode) # could not see line, so just cleaned lightly # within continuum region for all channels at once # 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't see line so will use same velocity range as 13co32 # 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 = 494.79 mJy, rms = 75.97 mJy, S/N = 6.5 # 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 G/f18_c18o32.mom0_cropped.fits Assuming object center (300.0,300.0) Background: 0.00 mJy/beam km/s RMS in annulus 4.0-9.0 arcsec = 78.44 mJy/beam km/s i radius flux err snr (asec) (mJy) (mJy) 0 0.1 26.2 20.9 1.3 1 0.2 92.7 54.6 1.7 2 0.3 188.1 84.1 2.2 <-- ND 3 0.4 286.3 94.9 3.0 4 0.5 387.0 127.6 3.0 5 0.6 454.4 165.3 2.7 <-- flux levels off 6 0.7 457.0 165.5 2.8 7 0.8 447.8 179.4 2.5 8 0.9 475.9 212.6 2.2 9 1.0 524.8 217.3 2.4 10 1.1 541.1 205.2 2.6 11 1.2 514.4 300.0 1.7 12 1.3 472.6 248.9 1.9 13 1.4 406.9 245.0 1.7 14 1.5 315.6 269.7 1.2 F = 541.1 mJy E = 205.2 mJy S = 2.6 D = 2.2 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 = 500, threshold = 0, interactive = False, mask = 'f18_cont_mask.crtf', imsize = imsize, cell = cell, weighting ='briggs', robust = robust, imagermode = imagermode) # could not see line clearly, so just cleaned lightly # within continuum region for all channels at once # 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'}]) # can see peaks at 2 and 6 km/s # 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,6km/s]')) # CHANGEME (based on above analysis!!) # can see line emission # 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 = 641.22 mJy, rms = 34.28 mJy, S/N = 18.7 # 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') ''' Measuring COG for G/f18_cn32.mom0_cropped.fits Assuming object center (300.0,300.0) Background: 0.14 mJy/beam km/s RMS in annulus 4.0-9.0 arcsec = 34.78 mJy/beam km/s i radius flux err snr (asec) (mJy) (mJy) 0 0.1 13.9 8.9 1.6 1 0.2 62.3 20.1 3.1 2 0.3 166.9 41.4 4.0 3 0.4 302.5 43.4 7.0 4 0.5 448.8 63.1 7.1 5 0.6 570.5 69.9 8.2 6 0.7 661.7 87.7 7.5 7 0.8 695.7 99.5 7.0 8 0.9 679.9 113.7 6.0 9 1.0 664.3 113.5 5.9 10 1.1 663.6 170.8 3.9 11 1.2 655.4 131.0 5.0 12 1.3 621.5 137.4 4.5 13 1.4 568.1 155.1 3.7 14 1.5 519.2 142.0 3.7 F = 695.7 mJy E = 99.5 mJy S = 7.0 D = 1.6 arcsec '''