#!/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_80 M4.8 II # 16:09:01.85 -39:05:12.5 field = 42 # 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 = '-3km/s' nchan = 16 xc = 322 # CHANGEME yc = 311 # CHANGEME in_a = 80 out_a = 120 aper = 0.5 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 = 'f42_cont_mask.crtf', imsize = imsize, cell = cell, weighting ='briggs', robust = robust, imagermode = imagermode) # used continuum region in automated mode # can see emission but faint and extended so not clear where to clean # 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'}]) # maybe see gas # 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.image', # CHANGEME (based on above analysis!!) outfile = 'f'+str(field)+'_13co32.mom0', moments = [0], includepix = [-10.0,100.0], chans = ('range=[2km/s,6km/s]')) # CHANGEME (based on above analysis!!) # export to fits file fbase = 'f'+str(field)+'_13co32.mom0' os.system('rm -rf '+fbase+'.fits') exportfits(imagename=fbase,fitsimage=fbase+'.fits') # 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 = 343.68 mJy, rms = 29.33 mJy, S/N = 11.7 # 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/f42_13co32.mom0_cropped.fits Assuming object center (300.0,300.0) Background: 0.40 mJy/beam km/s RMS in annulus 4.0-6.0 arcsec = 30.01 mJy/beam km/s i radius flux err snr (asec) (mJy) (mJy) 0 0.10 29.84 9.33 3.2 1 0.20 97.09 21.69 4.5 2 0.30 181.28 32.84 5.5 3 0.40 251.09 41.46 6.1 4 0.50 343.68 53.84 6.4 5 0.60 465.53 59.04 7.9 6 0.70 562.24 79.62 7.1 7 0.80 584.89 71.52 8.2 8 0.90 581.85 78.53 7.4 9 1.00 618.47 88.17 7.0 10 1.10 680.36 93.12 7.3 11 1.20 716.21 105.43 6.8 12 1.30 697.79 102.05 6.8 13 1.40 666.94 106.43 6.3 14 1.50 632.96 143.34 4.4 15 1.60 580.58 115.15 5.0 16 1.70 539.61 107.86 5.0 17 1.80 531.85 122.76 4.3 18 1.90 543.96 95.09 5.7 19 2.00 574.46 156.56 3.7 F = 716.21 mJy E = 105.43 mJy S = 6.79 D = 2.40 arcsec # LARGE APERTURE MAY BE REAL # GAS LOOKS EXTENDED ''' # ======================== C18O line ========================================== # don't bother with selfcal as it doesn't help... # os.system('rm -rf f'+str(field)+'_c18o32_b4sc*') clean(vis = 'f'+str(field)+'.vis.contsub', imagename = 'f'+str(field)+'_c18o32_b4sc', mode = 'velocity', start = start, width = width, nchan = nchan, outframe = outframe, veltype = veltype, restfreq = '329.33055GHz', niter = 500, threshold = 0, interactive = False, mask = 'f42_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_b4sc.image'}], contour = [{'file':'f'+str(field)+'_cont.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,6km/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 = -32.58 mJy, rms = 34.75 mJy, S/N = -0.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') # COULD NOT SEE GAS, SO USED 1" APER TO MEASURE FLUX # 5 0.500 877 -28.98 59.410 -0.5 49.45 33.612 1.5 0.999 ''' Measuring COG for M/f42_c18o32.mom0_cropped.fits Assuming object center (300.0,300.0) Background: 0.27 mJy/beam km/s RMS in annulus 4.0-6.0 arcsec = 34.39 mJy/beam km/s i radius flux err snr (asec) (mJy) (mJy) 0 0.10 10.97 9.19 1.2 1 0.20 27.01 21.57 1.3 2 0.30 20.37 36.58 0.6 <--- ND 3 0.40 -14.23 45.58 -0.3 4 0.50 -32.58 56.42 -0.6 5 0.60 -11.35 68.51 -0.2 6 0.70 23.21 67.71 0.3 7 0.80 75.33 74.08 1.0 8 0.90 125.44 91.45 1.4 9 1.00 177.49 98.57 1.8 10 1.10 222.89 102.21 2.2 11 1.20 216.82 119.03 1.8 12 1.30 141.30 153.40 0.9 13 1.40 42.35 142.53 0.3 14 1.50 -50.02 164.14 -0.3 15 1.60 -104.75 175.21 -0.6 16 1.70 -142.66 161.23 -0.9 17 1.80 -174.34 171.97 -1.0 18 1.90 -180.72 147.28 -1.2 19 2.00 -189.41 157.17 -1.2 F = 222.89 mJy E = 102.21 mJy S = 2.18 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=[3km/s,6km/s]')) # CHANGEME (based on above analysis!!) # can see 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 = 858.25 mJy, rms = 20.34 mJy, S/N = 42.2 # 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.0 1648.64 144.098 11.4 diameter total rms snr 0.6 472.58 24.050 19.6 diameter peak rms snr 0.2 231.89 20.068 11.6 '''