#!/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_141 16:13:44.11 -37:36:46.4 field = 24 # 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 = 334 # CHANGEME yc = 313 # 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 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 = 'f24_cont_mask.crtf', imsize = imsize, cell = cell, weighting ='briggs', robust = robust, imagermode = imagermode) # could not see line, so just cleaned in automated mode around continuum # 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.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]')) # can't see emission but this gave 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 = 381.48 mJy, rms = 65.29 mJy, S/N = 5.8 # 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/f24_13co32.mom0_cropped.fits Assuming object center (300.0,300.0) Background: 0.33 mJy/beam km/s RMS in annulus 4.0-9.0 arcsec = 64.69 mJy/beam km/s i radius flux err snr (asec) (mJy) (mJy) 0 0.1 47.9 20.5 2.3 1 0.2 141.5 36.6 3.9 2 0.3 221.1 72.7 3.0 3 0.4 281.2 85.9 3.3 4 0.5 381.5 106.4 3.6 5 0.6 468.4 108.0 4.3 6 0.7 459.1 150.5 3.1 7 0.8 411.8 150.4 2.7 8 0.9 366.0 173.6 2.1 9 1.0 333.8 168.4 2.0 F = 468.4 mJy E = 108.0 mJy S = 4.3 D = 1.2 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 = 'f24_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 # 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') # 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'}]) # 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=[1km/s,7km/s]')) # cannot see emission, 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 = -351.47 mJy, rms = 84.30 mJy, S/N = -4.2 # 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/f24_c18o32.mom0_cropped.fits Assuming object center (300.0,300.0) Background: -0.62 mJy/beam km/s RMS in annulus 4.0-6.0 arcsec = 79.58 mJy/beam km/s i radius flux err snr (asec) (mJy) (mJy) 0 0.10 -27.50 19.26 -1.4 1 0.20 -103.87 55.03 -1.9 2 0.30 -206.38 85.78 -2.4 <-- ND 3 0.40 -291.18 99.62 -2.9 4 0.50 -351.47 127.71 -2.8 5 0.60 -323.70 145.56 -2.2 6 0.70 -238.70 201.06 -1.2 7 0.80 -266.39 168.29 -1.6 F = -27.50 mJy E = 19.26 mJy S = -1.43 D = 0.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 = 500, threshold = 0, interactive = False, mask = 'f24_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=[0km/s,6km/s]')) # CHANGEME (based on above analysis!!) # difficult to see emission, this gives best 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 = 235.00 mJy, rms = 45.72 mJy, S/N = 5.1 # 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/f24_cn32.mom0_cropped.fits Assuming object center (300.0,300.0) Background: 0.28 mJy/beam km/s RMS in annulus 4.0-9.0 arcsec = 46.45 mJy/beam km/s i radius flux err snr (asec) (mJy) (mJy) 0 0.1 19.1 10.2 1.9 1 0.2 53.5 30.1 1.8 2 0.3 97.7 57.0 1.7 3 0.4 163.0 61.0 2.7 4 0.5 235.0 78.8 3.0 5 0.6 267.2 112.9 2.4 6 0.7 269.1 125.4 2.1 7 0.8 267.7 136.7 2.0 8 0.9 243.7 141.8 1.7 9 1.0 178.9 148.6 1.2 F = 269.1 mJy E = 125.4 mJy S = 2.1 D = 1.4 arcsec '''