#!/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 ===================================== # Lupus I_14 15:45:08.88 -34:17:33.7 # Class II, K0 field = 6 # 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 = '0.5km/s' start = '0km/s' nchan = 17 xc = 331 # CHANGEME yc = 316 # CHANGEME in_a = 80 out_a = 120 aper = 0.5 # 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 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 = 500, interactive = True, threshold = 0, 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':'test_f'+str(field)+'_13co32.image'}], contour = [{'file':'f'+str(field)+'_cont.fits'}]) # 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', # CHANGEME (based on above analysis!!) outfile = 'test_f'+str(field)+'_13co32.mom0', moments = [0], includepix = [-10.0,100.0], chans = ('range=[3km/s,6km/s]')) # can see emission, gives highest snr # first moment map sigma = 30e-3 # Jy/beam in peak velocity channel os.system('rm -rf test_f6_13co32_mom1.image') immoments(imagename = 'test_f6_13co32.cube.fits', outfile = 'test_f6_13co32_mom1.image', moments = [1], includepix = [3.0*sigma,100.0], chans = ('range=[3km/s,6km/s]')) os.system('rm -rf test_f6_13co32_mom1.fits') exportfits(imagename='test_f6_13co32_mom1.image',fitsimage='test_f6_13co32_mom1.fits') os.system('rm -rf test_f6_13co32_mom1.image') # 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_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 = 739.23 mJy, rms = 45.14 mJy, S/N = 16.4 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/f6_13co32.mom0_cropped.fits Assuming object center (300.0,300.0) Background: -0.26 mJy/beam km/s RMS in annulus 4.0-6.0 arcsec = 49.67 mJy/beam km/s i radius flux err snr (asec) (mJy) (mJy) 0 0.10 51.56 11.20 4.6 1 0.20 177.60 31.99 5.6 2 0.30 381.18 45.61 8.4 3 0.40 591.31 63.59 9.3 4 0.50 739.23 84.64 8.7 5 0.60 758.95 87.34 8.7 6 0.70 694.70 86.16 8.1 7 0.80 608.69 124.18 4.9 8 0.90 576.47 129.47 4.5 9 1.00 622.20 114.85 5.4 F = 758.95 mJy E = 87.34 mJy S = 8.69 D = 1.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 = 'f6_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.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.image', # CHANGEME (based on above analysis!!) outfile = 'f'+str(field)+'_c18o32.mom0', moments = [0], includepix = [-10.0,100.0], chans = ('range=[3km/s,6km/s]')) # emission not clear, used 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 = 358.91 mJy, rms = 59.33 mJy, S/N = 6.0 # 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/f6_c18o32.mom0_cropped.fits Assuming object center (300.0,300.0) Background: 0.20 mJy/beam km/s RMS in annulus 4.0-6.0 arcsec = 60.72 mJy/beam km/s i radius flux err snr (asec) (mJy) (mJy) 0 0.10 29.32 15.72 1.9 1 0.20 112.83 35.83 3.1 2 0.30 231.04 61.38 3.8 3 0.40 308.17 81.73 3.8 4 0.50 358.91 127.24 2.8 5 0.60 421.55 129.92 3.2 6 0.70 509.85 127.93 4.0 7 0.80 573.33 145.52 3.9 <--where levels off 8 0.90 580.29 165.97 3.5 9 1.00 518.74 150.03 3.5 10 1.10 418.38 164.57 2.5 11 1.20 345.35 184.34 1.9 12 1.30 344.91 224.90 1.5 13 1.40 338.75 195.50 1.7 14 1.50 255.58 181.82 1.4 F = 580.29 mJy E = 165.97 mJy S = 3.50 D = 1.80 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 = 'f6_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)+'_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,6km/s]')) # CHANGEME (based on above analysis!!) # can't see line, max 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 = 778.89 mJy, rms = 37.84 mJy, S/N = 20.6 # 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/f6_cn32.mom0_cropped.fits Assuming object center (300.0,300.0) Background: 0.07 mJy/beam km/s RMS in annulus 4.0-9.0 arcsec = 39.86 mJy/beam km/s i radius flux err snr (asec) (mJy) (mJy) 0 0.1 39.1 9.1 4.3 1 0.2 135.0 24.6 5.5 2 0.3 269.7 40.1 6.7 3 0.4 391.1 50.7 7.7 4 0.5 480.2 75.8 6.3 5 0.6 515.9 102.2 5.0 <-- used this (flux levels off; looks best in CASA) 6 0.7 544.9 96.2 5.7 7 0.8 589.4 119.3 4.9 8 0.9 624.8 120.0 5.2 9 1.0 655.5 127.3 5.2 10 1.1 716.0 152.5 4.7 11 1.2 776.4 193.5 4.0 F = 776.4 mJy E = 193.5 mJy S = 4.0 D = 2.4 arcsec '''