############################################################################## # Data reduction script for VY CMa Band 7 321-GHz 12-m data: Part 2: Imaging # VYCMa_Band7_321GHz_Imaging.py # Tested using CASA Version 4.2.1 (r29048) # # Datasets used: # # uid___A002_X6cdf83_X1118 - observed Aug 16 2013 with 16 antennas # Flux scale calibrator J0522-364 Bandpass calibrator J0522-3627 (same # as J0522-364) Phase reference source J0648-3044 Target Vy_CMa # # uid___A002_X6cdf83_X12ca - observed Aug 16 2013 with 16 antennas # Flux scale calibrator Pallas Bandpass calibrator J0522-3627 Phase # reference source J0648-3044 Target Vy_CMa # # uid___A002_X6d0f96_X2051 - observed Aug 17 2013 with 20 antennas # Flux scale calibrator Pallas # Bandpass calibrator J0522-3627 # Phase reference source J0648-3044 # Target Vy_CMa # ############################################################################## """ See accompanying README file for details of the necessary input files and comments on the data. Before running this script you should have run VYCMa_Band7_321GHz_Calibration.py or obtained the calibrated target dataset VYCMa_Band7_321.ms. This script assumes that VYCMa_321.ms has had all instrumental and reference source calibration applied, and has been adjusted to the LSRK frame, as at the end of VYCMa_Band7_321GHz_Calibration.py Prior to running this script you will need to obtain the Analysis Utilities package. See http://casaguides.nrao.edu/index.php?title=Analysis_Utilities for details and for how to download Analysis Utilities. """ #---------------------------------------------------------------------------------- #----- Optional Steps ------------------------------------------------------------- #---------------------------------------------------------------------------------- #----- Set this option to False if you do not want to do interactive clean. #----- Note that better results will be obtained if you do make masks interatively, #----- expanding the area as fainter emission emerges cleaninteractive=F #----- Set this option to true if you want to make diagnostic plots along the way. #----- Note that this will slow down the reduction. #----- Default is no plots (you can make them later anyway). calplots=T #---------------------------------------------------------------------------------- #----- Some setup steps ----------------------------------------------------------- #---------------------------------------------------------------------------------- version = casalog.version() print "You are using " + version if (version < '4.2.1'): print "\033[91m YOUR VERSION OF CASA IS TOO OLD FOR THIS GUIDE." print "\033[91m PLEASE UPDATE IT BEFORE PROCEEDING." else: print "Your version of CASA is appropriate for this guide." # Repeat this if you are using copy and paste, and you re-start CASA refantenna='DV15' thesteps=[] step_title = {0: 'List the data set and plot the visibility spectrum', 1: 'Make first image of maser peak', 2: 'Self-calibrate cycle 1, phase-only, re-image', 3: 'Self-calibrate cycle 2, phase-only, re-image', 4: 'Self-calibrate cycle 3, amp & phase, re-image', 5: 'Self-calibrate cycle 4, amp & phase, re-image', 6: 'Make low-res cubes to identify continuum \n(this step can be omitted if you trust our channel selection)', 7: 'Image continuum', 8: 'Subtract continuum', 9: 'Image a line in spw 0', 10: 'Image a line in spw 1', 11: 'Image a line in spw 2', 12: 'Image a line in spw 3' } # The Python variable 'mysteps' will control which steps # are executed when you start the script using # execfile('scriptForCalibration.py') # e.g. setting # mysteps = [2,3,4]# before starting the script will make the script execute # only steps 2, 3, and 4 # Setting mysteps = [] will make it execute all steps. try: print 'List of steps to be executed ...', mysteps thesteps = mysteps except: print 'global variable mysteps not set.' if (thesteps==[]): thesteps = range(0,len(step_title)) print 'Executing all steps: ', thesteps ###################################################### # Do listobs and write info to a text file; plot spectrum and brightest channel mystep = 0 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] os.system('rm VYCMa_Band7_321_listobs.txt') listobs(vis='VYCMa_Band7_321.ms', listfile='VYCMa_Band7_321_listobs.txt', verbose=True) print '<< Printed listobs output to VYCMa_Band7_321_listobs.txt' for s in ['0','1','2','3']: os.system('rm -rf VYCMa_Band7_321_spw'+s+'_1st-vis-spectrum.png') plotms(vis='VYCMa_Band7_321.ms', xaxis='frequency', yaxis='amp', selectdata=T, spw=s, avgtime='1e8',avgscan=T, coloraxis='baseline', showgui=F, plotfile='VYCMa_Band7_321_spw'+s+'_1st-vis-spectrum.png') # You will notice the very bright maser line in spw 0:1958 # You will also see that the end channels need flagging in spw 2, 3 # The appropriate selections are made below for self-calibration os.system('rm -rf VYCMa_Band7_321_0-1958_1st_amp.png') plotms(vis='VYCMa_Band7_321.ms', xaxis='time', yaxis='amp', selectdata=T, antenna='DV15&*', spw='0:1958', coloraxis='baseline', showgui=F, plotfile='VYCMa_Band7_321_0-1958_1st_amp.png') os.system('rm -rf VYCMa_Band7_321_0-1958_1st_phase.png') plotms(vis='VYCMa_Band7_321.ms', xaxis='time', yaxis='phase', selectdata=T, antenna='DV15&*', spw='0:1958', coloraxis='baseline', showgui=F, plotfile='VYCMa_Band7_321_0-1958_1st_phase.png') # From listobs, the scan numbers are 10,14,18,22,26,27,31,35,39,43,44,48,52 # You may want to plot baselines iteratively to check, but it turns out that # * The last few integrations of DA41 and DV24 scan 31 are irredeemably bad # * DV19 scan 10 has very fast phase # * The first scan has high amplitudes ###################################################### # Make first image of maser peak, omitting dubious data mystep = 1 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] os.system('rm -rf VYCMa_Band7_321_0_1958.clean*') clean(vis = 'VYCMa_Band7_321.ms', imagename='VYCMa_Band7_321_0_1958.clean', field='Vy_CMa', spw='0:1958', scan='14~52', mode='channel',start=1958,nchan=1, cell='0.01arcsec',npercycle=10,imsize=512, psfmode='hogbom', mask='ellipse[[256pix,275pix],[35pix,28pix],90deg]', interactive=cleaninteractive, niter=50) # Get image statistics imin='VYCMa_Band7_321_0_1958.clean.image' noise=imstat(imagename=imin, box='250,50,500,200')['rms'][0] peak=imstat(imagename=imin, box='200,200,330,330')['max'][0] print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise) #rms 6.734, peak 322.546, snr 47 ###################################################### # Self-calibrate, phase-only first, using all data mystep = 2 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] os.system('rm -rf cal-VYCMa_321_1958.ph1') gaincal(vis = 'VYCMa_Band7_321.ms', field='Vy_CMa', refant=refantenna, caltable='cal-VYCMa_321_1958.ph1', spw='0:1958', calmode='p', solint='60s') if (calplots): os.system('rm -rf cal-VYCMa_321_1958.ph1*.png') for ant in range(3): plotcal(caltable='cal-VYCMa_321_1958.ph1', xaxis = 'time', yaxis = 'phase', antenna=str(ant*8)+'~'+str(ant*8+7), subplot=421, fontsize=7.0, iteration='antenna', figfile='cal-VYCMa_321_1958.ph1_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png') applycal(vis = 'VYCMa_Band7_321.ms', field='Vy_CMa', spw='0,1,2,3', spwmap=['0,0,0,0'], gaintable='cal-VYCMa_321_1958.ph1', calwt = F, flagbackup = F) os.system('rm -rf VYCMa_Band7_321_0_1958_ph1.clean*') clean(vis = 'VYCMa_Band7_321.ms', imagename='VYCMa_Band7_321_0_1958_ph1.clean', field='Vy_CMa', spw='0:1958', scan='14~52', mode='channel',start=1958,nchan=1, cell='0.01arcsec',npercycle=50,imsize=512, psfmode='hogbom', mask='ellipse[[244pix,273pix],[44pix,38pix],90deg]', interactive=cleaninteractive, niter=500) # Get image statistics imin='VYCMa_Band7_321_0_1958_ph1.clean.image' noise=imstat(imagename=imin, box='250,50,500,200')['rms'][0] peak=imstat(imagename=imin, box='200,200,330,330')['max'][0] print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise) #rms 1.150, peak 423.971, snr 369 ###################################################### # Self-calibrate cycle 2, phase-only, shorter solint, using all data mystep = 3 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] os.system('rm -rf cal-VYCMa_321_1958.ph2') gaincal(vis = 'VYCMa_Band7_321.ms', field='Vy_CMa', refant=refantenna, caltable='cal-VYCMa_321_1958.ph2', spw='0:1958', calmode='p', solint='30s') if (calplots): os.system('rm -rf cal-VYCMa_321_1958.ph2*.png') for ant in range(3): plotcal(caltable='cal-VYCMa_321_1958.ph2', xaxis = 'time', yaxis = 'phase', antenna=str(ant*8)+'~'+str(ant*8+7), subplot=421, fontsize=7.0, iteration='antenna', figfile='cal-VYCMa_321_1958.ph2_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png') applycal(vis = 'VYCMa_Band7_321.ms', field='Vy_CMa', spw='0,1,2,3', spwmap=['0,0,0,0'], gaintable='cal-VYCMa_321_1958.ph2', calwt = F, flagbackup = F) os.system('rm -rf VYCMa_Band7_321_0_1958_ph2.clean*') clean(vis = 'VYCMa_Band7_321.ms', imagename='VYCMa_Band7_321_0_1958_ph2.clean', field='Vy_CMa', spw='0:1958', scan='14~52', mode='channel',start=1958,nchan=1, cell='0.01arcsec',npercycle=50,imsize=512, psfmode='hogbom', mask='ellipse[[248pix,275pix],[57pix,44pix],90deg]', interactive=cleaninteractive, niter=500) # Get image statistics imin='VYCMa_Band7_321_0_1958_ph2.clean.image' noise=imstat(imagename=imin, box='250,50,500,200')['rms'][0] peak=imstat(imagename=imin, box='200,200,330,330')['max'][0] print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise) #rms 0.883, peak 446.188, snr 505 ###################################################### # Self-calibrate cycle 3, a&p, using all data mystep = 4 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] os.system('rm -rf cal-VYCMa_321_1958.ap1') gaincal(vis = 'VYCMa_Band7_321.ms', field='Vy_CMa', refant=refantenna, caltable='cal-VYCMa_321_1958.ap1', spw='0:1958', gaintable='cal-VYCMa_321_1958.ph2', spwmap=['0,0,0,0'], calmode='ap', solint='120s') if (calplots): os.system('rm -rf cal-VYCMa_321_1958.ap1*.png') for ant in range(3): plotcal(caltable='cal-VYCMa_321_1958.ap1', xaxis = 'time', yaxis = 'amp', antenna=str(ant*8)+'~'+str(ant*8+7), subplot=421, fontsize=7.0, iteration='antenna', figfile='cal-VYCMa_321_1958.ap1_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png') applycal(vis = 'VYCMa_Band7_321.ms', field='Vy_CMa', spw='0,1,2,3', spwmap=[['0,0,0,0'],['0,0,0,0']], gaintable=['cal-VYCMa_321_1958.ph2','cal-VYCMa_321_1958.ap1'], calwt = F, flagbackup = F) # Now include all scans os.system('rm -rf VYCMa_Band7_321_0_1958_ap1.clean*') clean(vis = 'VYCMa_Band7_321.ms', imagename='VYCMa_Band7_321_0_1958_ap1.clean', field='Vy_CMa', spw='0:1958', mode='channel',start=1958,nchan=1, cell='0.01arcsec',npercycle=50,imsize=512, psfmode='hogbom', mask=['ellipse[[249pix,277pix],[66pix,49pix],90deg]','ellipse[[319pix,225pix],[44pix,27pix],0deg]'], interactive=cleaninteractive, niter=500) # Get image statistics imin='VYCMa_Band7_321_0_1958_ap1.clean.image' noise=imstat(imagename=imin, box='250,50,500,200')['rms'][0] peak=imstat(imagename=imin, box='200,200,330,330')['max'][0] print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise) #rms 0.697, peak 446.022, snr 640 ###################################################### # Self-calibrate cycle 4, ap, short solint, using all data mystep = 5 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] os.system('rm -rf cal-VYCMa_321_1958.ap2') gaincal(vis = 'VYCMa_Band7_321.ms', field='Vy_CMa', refant=refantenna, caltable='cal-VYCMa_321_1958.ap2', spw='0:1958', gaintable='cal-VYCMa_321_1958.ph2', spwmap=['0,0,0,0'], calmode='ap', solint='30s') if (calplots): os.system('rm -rf cal-VYCMa_321_1958.ap2*.png') for ant in range(3): plotcal(caltable='cal-VYCMa_321_1958.ap2', xaxis = 'time', yaxis = 'amp', antenna=str(ant*8)+'~'+str(ant*8+7), subplot=421, fontsize=7.0, iteration='antenna', figfile='cal-VYCMa_321_1958.ap2_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png') applycal(vis = 'VYCMa_Band7_321.ms', field='Vy_CMa', spw='0,1,2,3', spwmap=[['0,0,0,0'],['0,0,0,0']], gaintable=['cal-VYCMa_321_1958.ph2','cal-VYCMa_321_1958.ap2'], calwt = F, flagbackup = F) # os.system('rm -rf VYCMa_Band7_321_0_1958_ap2.clean*') clean(vis = 'VYCMa_Band7_321.ms', imagename='VYCMa_Band7_321_0_1958_ap2.clean', field='Vy_CMa', spw='0:1958', mode='channel',start=1958,nchan=1, cell='0.01arcsec',npercycle=50,imsize=512, psfmode='hogbom', mask=['ellipse[[247pix,305pix],[81pix,72pix],0deg]','ellipse[[317pix,249pix],[69pix,32pix],0deg]'], interactive=cleaninteractive, niter=1000) # Get image statistics imin='VYCMa_Band7_321_0_1958_ap2.clean.image' noise=imstat(imagename=imin, box='250,50,500,200')['rms'][0] peak=imstat(imagename=imin, box='200,200,330,330')['max'][0] print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise) #rms 0.533, peak 445.955, snr 836 ##################################### # Image whole cube at lower resolution to identify line-free channels mystep = 6 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] for s in ['0','1','2','3']: os.system('rm -rf VYCMa_Band7_321_spw'+s+'_5ap2.clean*') clean(vis = 'VYCMa_Band7_321.ms', imagename='VYCMa_Band7_321_spw'+s+'_5ap2.clean', field='Vy_CMa', spw=s, mode='channel',width=5, cell='0.02arcsec',imsize=512,usescratch=T, psfmode='hogbom', interactive=cleaninteractive, niter=1000) # The following is our best guess at a selection contchans='0:10~219;360~1799;3320~3819,1:190~1049;1540~1674;2180~2364;2710~2874;3600~3819,2:265~754;1315~1714;1810~2174;2905~3249;3440~3599,3:125~264;365~814;935~999;1500~1759;2250~2499;3095~3485;3540~3674' ############################################## # Image continuum mystep = 7 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] os.system('rm -rf VYCMa_Band7_321_contap2.clean*') clean(vis = 'VYCMa_Band7_321.ms', imagename='VYCMa_Band7_321_contap2.clean', field='Vy_CMa', spw=contchans, mode='mfs', multiscale=[0,4,8,12] , mask=['ellipse[[218pix,255pix],[64pix,52pix],90deg]','ellipse[[281pix,359pix],[140pix,70pix],20deg]'], cell='0.01arcsec',imsize=512, psfmode='hogbom',interactive=cleaninteractive) exportfits(imagename='VYCMa_Band7_321_contap2.clean.image', fitsimage='VYCMa_Band7_321_contap2.clean.fits') ##################################### # Subtract continuum mystep = 8 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] os.system('rm -rf VYCMa_Band7_321.ms.contsub') uvcontsub(vis = 'VYCMa_Band7_321.ms',field='V*',fitorder=1, combine='spw',solint='30s',spw='',fitspw=contchans) #---------------------------------------------------------------------------------- #----- Image maser and line emission and make moment maps ------------------------- #---------------------------------------------------------------------------------- # Some lines from Kaminski et al. 2013 ApJS 209 38 are imaged below. # However, note that it is possible that the named line may be blended with other lines, # especially for the thermal lines. Conversely, it is possible (especially # for the maser) that faint wings exceed the range imaged. # ##################################### # Image H2O 10(2,9)-9(3,6) line spw 0 mystep = 9 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] # Line frequencies taken from Kaminski et al. 2013 ApJS 209 38 # High resolution for maser os.system('rm -rf VYCMa_321.22_H2O.clean.*') clean(vis = 'VYCMa_Band7_321.ms.contsub', imagename='VYCMa_321.22_H2O.clean', field='Vy_CMa', spw='0', mode='channel',start=1825,nchan=255, restfreq='321.22564GHz', cell='0.01arcsec',imsize=512, psfmode='hogbom', weighting='briggs',robust=0.5, threshold='10mJy', mask=['ellipse[[247pix,305pix],[81pix,72pix],0deg]','ellipse[[317pix,249pix],[69pix,32pix],0deg]'], interactive=cleaninteractive, usescratch=T, niter=5000) # Make moment map thisimage='VYCMa_321.22_H2O.clean.image' chanstat=imstat(imagename=thisimage,box='250,50,500,100',chans='3') rms1= chanstat['rms'][0] chanstat=imstat(imagename=thisimage,box='250,50,500,100',chans='252') rms2= chanstat['rms'][0] rms=0.5*(rms1+rms2) print 'rms in a channel = '+str(rms) os.system('rm -rf '+thisimage+'.mom0') immoments(imagename = thisimage, moments = [0], axis = 'spectral', includepix = [rms*2,1000.], outfile = thisimage+'.mom0') # Export as fits exportfits(imagename='VYCMa_321.22_H2O.clean.image', fitsimage='VYCMa_321.22_H2O.clean.image.fits') exportfits(imagename=thisimage+'.mom0', fitsimage=thisimage+'.mom0.fits') # Note that since the dynamic range of the maser cube is 500~1000 around the peak, # implying a noise level more than an order of magnitude # higher than that for the fainter line wings, we don't make a first moment map. # However, it may be possible to obtain a meaningful first moment map with # additional data selection and/or blanking ##################################### # Image Si18O 8-7 line spw 1 mystep = 10 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] # Line frequencies taken from Kaminski et al. 2013 ApJS 209 38 # Lower resolution for weaker lines os.system('rm -rf VYCMa_322.77_Si18O.clean*') clean(vis = 'VYCMa_Band7_321.ms.contsub', imagename='VYCMa_322.77_Si18O.clean', field='Vy_CMa', spw='1', mode='channel',start=2850,nchan=100,width=5, restfreq='322.770135GHz', cell='0.01arcsec',imsize=512, psfmode='hogbom', threshold='4mJy', mask='ellipse[[284pix,263pix],[150pix,115pix],20deg]', interactive=cleaninteractive, usescratch=T, niter=1000) # Make moment maps thisimage='VYCMa_322.77_Si18O.clean.image' chanstat=imstat(imagename=thisimage,box='250,50,500,100',chans='3') rms1= chanstat['rms'][0] chanstat=imstat(imagename=thisimage,box='250,50,500,100',chans='97') rms2= chanstat['rms'][0] rms=0.5*(rms1+rms2) print 'rms in a channel = '+str(rms) os.system('rm -rf '+thisimage+'.mom0') immoments(imagename = thisimage, moments = [0], axis = 'spectral', includepix = [rms*2,100.], outfile = thisimage+'.mom0') os.system('rm -rf '+thisimage+'.mom1') immoments(imagename = thisimage, moments = [1], axis = 'spectral', includepix = [rms*5.5,100.], outfile = thisimage+'.mom1') # Export as fits exportfits(imagename='VYCMa_322.77_Si18O.clean.image', fitsimage='VYCMa_322.77_Si18O.clean.image.fits') exportfits(imagename=thisimage+'.mom0', fitsimage=thisimage+'.mom0.fits') exportfits(imagename=thisimage+'.mom1', fitsimage=thisimage+'.mom1.fits') ##################################### # Image TiO2 23(1,23)-22(0,22) line spw 2 mystep = 11 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] # Line frequencies taken from Kaminski et al. 2013 ApJS 209 38 # Lower resolution for weaker lines os.system('rm -rf VYCMa_310.78_TiO2.clean*') clean(vis = 'VYCMa_Band7_321.ms.contsub', imagename='VYCMa_310.78_TiO2.clean', field='Vy_CMa', spw='2', mode='channel',start=825,nchan=75,width=5, restfreq='310.7827126GHz', cell='0.01arcsec',imsize=512, threshold='4mJy', psfmode='hogbom', mask='ellipse[[284pix,263pix],[150pix,115pix],20deg]', interactive=cleaninteractive, usescratch=T, niter=1000) # Make moment maps thisimage='VYCMa_310.78_TiO2.clean.image' chanstat=imstat(imagename=thisimage,box='250,50,500,100',chans='3') rms1= chanstat['rms'][0] chanstat=imstat(imagename=thisimage,box='250,50,500,100',chans='72') rms2= chanstat['rms'][0] rms=0.5*(rms1+rms2) print 'rms in a channel = '+str(rms) os.system('rm -rf '+thisimage+'.mom0') immoments(imagename = thisimage, moments = [0], axis = 'spectral', includepix = [rms*2,100.], outfile = thisimage+'.mom0') os.system('rm -rf '+thisimage+'.mom1') immoments(imagename = thisimage, moments = [1], axis = 'spectral', includepix = [rms*4.5,100.], outfile = thisimage+'.mom1') # Export as fits exportfits(imagename='VYCMa_310.78_TiO2.clean.image', fitsimage='VYCMa_310.78_TiO2.clean.image.fits') exportfits(imagename=thisimage+'.mom0', fitsimage=thisimage+'.mom0.fits') exportfits(imagename=thisimage+'.mom1', fitsimage=thisimage+'.mom1.fits') ##################################### # Image NaCl v=1 24-23 line spw 3 mystep = 12 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] # Line frequencies taken from Kaminski et al. 2013 ApJS 209 38 # Lower resolution for weaker lines os.system('rm -rf VYCMa_309.79_NaCl.clean*') clean(vis = 'VYCMa_Band7_321.ms.contsub', imagename='VYCMa_309.79_NaCl.clean', field='Vy_CMa', spw='3', mode='channel',start=1000,nchan=35,width=5, restfreq='309.7878261GHz', threshold='4mJy', cell='0.01arcsec',imsize=512, psfmode='hogbom', mask='ellipse[[284pix,263pix],[150pix,115pix],20deg]', interactive=cleaninteractive, usescratch=T, niter=1000) # Make moment maps myimage='VYCMa_309.79_NaCl.clean.image' chanstat=imstat(imagename=myimage,box='250,50,500,100', chans='3') rms1= chanstat['rms'][0] chanstat=imstat(imagename=myimage,box='250,50,500,100', chans='32') rms2= chanstat['rms'][0] rms=0.5*(rms1+rms2) print 'rms in a channel = '+str(rms) os.system('rm -rf VYCMa_309.79_NaCl.clean.image.mom0') immoments(imagename = 'VYCMa_309.79_NaCl.clean.image', moments = [0], axis = 'spectral', includepix = [rms*2,100.], outfile = 'VYCMa_309.79_NaCl.clean.image.mom0') os.system('rm -rf VYCMa_309.79_NaCl.clean.image.mom1') immoments(imagename = 'VYCMa_309.79_NaCl.clean.image', moments = [1], axis = 'spectral', includepix = [rms*4,100.], outfile = 'VYCMa_309.79_NaCl.clean.image.mom1') # Export as fits exportfits(imagename='VYCMa_309.79_NaCl.clean.image', fitsimage='VYCMa_309.79_NaCl.clean.image.fits') exportfits(imagename='VYCMa_309.79_NaCl.clean.image.mom0', fitsimage='VYCMa_309.79_NaCl.clean.image.mom0.fits') exportfits(imagename='VYCMa_309.79_NaCl.clean.image.mom1', fitsimage='VYCMa_309.79_NaCl.clean.image.mom1.fits') #---------------------------------------------------------------------------------- #----- End of imaging script. #----------------------------------------------------------------------------------