######################################################### # This script is for CASA Version 3.3 and WILL NOT WORK # # in CASA Version 3.4 # ######################################################### ############################### # Concat and fix a few issues # ############################### split_data = ['X1ed.decimated.ms.I16293', 'X39.decimated.ms.I16293', 'X3bf.decimated.ms.I16293', 'X4fb.decimated.ms.I16293'] # Zero pointing table (clean does not currently use these correctly) for vis in split_data: myMsName = vis tb.open(myMsName+'/POINTING', nomodify = False) a = tb.rownumbers() tb.removerows(a) tb.close() # This fixes a small problem with the X4fb source table. This # would not be necessary for future datasets tb.open('X39.decimated.ms.I16293/SOURCE') s=tb.getcol('DIRECTION') tb.close() tb.open('X4fb.decimated.ms.I16293/SOURCE',nomodify=False) tb.putcol('DIRECTION',s) tb.close() os.system('rm -rf I16293_corrected.ms') concat(vis=split_data,concatvis='I16293_corrected.ms') os.system('rm -rf I16293_corrected.listobs.txt') listobs(vis='I16293_corrected.ms', listfile='I16293_corrected.listobs.txt') ################################## # Use clean to make a dirty cube # ################################## os.system('rm -rf dirty*') clean(vis='I16293_corrected.ms', imagename='dirty', field='0,1', mode='channel', niter=0, imagermode='mosaic', phasecenter='0', cell=['0.2arcsec'], imsize=[220,220], outframe='LSRK') # Use the Viewer to find line-free channels # Then do a continuum subtraction os.system('rm -rf I16293_corrected.ms.contsub') uvcontsub(vis='I16293_corrected.ms', fitorder=0, fitspw='0:410~430') # If you want, check quality of continuum subtraction os.system('rm -rf dirty_contsub*') clean(vis='I16293_corrected.ms.contsub', imagename='dirty_contsub', field='0,1', mode='channel', niter=0, imagermode='mosaic', phasecenter='0', cell=['0.2arcsec'], imsize=[220,220], outframe='LSRK') # Looks pretty good, but it's very hard to find line-free channels, so # buyer beware ########################################################################## # Clean the continuum only channels # # Note that the rms per channel in the dirty cube is ~ 5mJy/beam # # Since we'll be averaging together 21 channels, the rms in our # # final image should be about 1 mJy/beam # # Set the clean threshold to 3 mJy/beam so that we're not cleaning noise # ########################################################################## # Select the continuum channels using spw. Use interactive and stop at # a conservative spot, so the model is conservative for # self-calibration os.system('rm -rf I16293_clean_continuum*') clean(vis='I16293_corrected.ms', imagename='I16293_clean_continuum', field='0,1', phasecenter='0', mode='mfs',spw='0:410~430', imagermode='mosaic', imsize=[220,220],cell='0.2arcsec', interactive=T, pbcor=False, weighting='briggs',robust=0.5, niter=10000, threshold='3.0mJy') ########################################## # Self-calibration on the continuum data # ########################################## # For the first iteration, get one solution per scan using # 'inf' with combine=''. # Only calibrate with the continuum channels, since they're the # only channels with useful model column data # Gaincal os.system('rm -rf I16293_corrected.ms.selfcal.pcal1') vis = 'I16293_corrected.ms' gaincal(vis='I16293_corrected.ms', caltable='I16293_corrected.ms.selfcal.pcal1', field='0,1',refant='DV03', gaintype='T', spw='0:410~430', calmode='p',solint='inf',combine='', minsnr=2.0,minblperant=4) # There are no failed solutions which is a good sign. # Plot solutions to check that they seem acceptable - # they should show smooth variations with time, and consistent # from antenna to antenna plotcal(caltable='I16293_corrected.ms.selfcal.pcal1', xaxis='time',yaxis='phase', iteration='antenna',subplot=421, plotrange=[0,0,-180,180]) # Apply the solutions applycal(vis='I16293_corrected.ms', gaintable=['I16293_corrected.ms.selfcal.pcal1'], field='', calwt=F) # Re-clean to see the improvement os.system('rm -rf I16293_clean_continuum_pcal1*') clean(vis='I16293_corrected.ms', imagename='I16293_clean_continuum_pcal1', field='0,1', phasecenter='0', mode='mfs',spw='0:410~430', imagermode='mosaic', imsize=[220,220],cell='0.2arcsec', interactive=T, pbcor=False, weighting='briggs',robust=0.5, niter=10000, threshold='3.0mJy') # Selfcal seems to have lowered the rms from 12 to 7 mJy/beam # The theorectical rms is this map is: 0.5 mJy # Try selfcal with shorter solint os.system('rm -rf I16293_corrected.ms.selfcal.pcal2') gaincal(vis='I16293_corrected.ms', caltable='I16293_corrected.ms.selfcal.pcal2', field='0,1',refant='DV03', gaintype='T', spw='0:410~430', calmode='p',solint='30s',minsnr=3.0,minblperant=4) # Plot solutions again plotcal(caltable='I16293_corrected.ms.selfcal.pcal2', xaxis='time',yaxis='phase', iteration='antenna',subplot=421, plotrange=[0,0,-180,180]) # Apply solutions applycal(vis='I16293_corrected.ms', gaintable=['I16293_corrected.ms.selfcal.pcal2'], field='', calwt=F) # Clean again, stopping clean when the image looks noise-like os.system('rm -rf I16293_clean_continuum_pcal2*') clean(vis='I16293_corrected.ms', imagename='I16293_clean_continuum_pcal2', field='0,1', phasecenter='0', mode='mfs',spw='0:410~430', imagermode='mosaic', imsize=[220,220],cell='0.2arcsec', interactive=T, pbcor=False, weighting='briggs',robust=0.5, niter=10000, threshold='1.5mJy') # Now try an amplitude solution. Use a long solint, we mostly # want to take out any variations between the different datasets. # You need to apply the best phase only self-cal os.system('rm -rf I16293_corrected.ms.selfcal.ampcal') gaincal(vis='I16293_corrected.ms', caltable='I16293_corrected.ms.selfcal.ampcal', field='0,1',refant='DV03', gaintype='T', calmode='ap',solint='1200s',combine='scan', spw='0:410~430', minsnr=3.0,minblperant=4,gaintable=['I16293_corrected.ms.selfcal.pcal2']) # Plot solutions again plotcal(caltable='I16293_corrected.ms.selfcal.ampcal', xaxis='time',yaxis='amp', iteration='antenna',subplot=421) # Apply solutions applycal(vis='I16293_corrected.ms', gaintable=['I16293_corrected.ms.selfcal.pcal2','I16293_corrected.ms.selfcal.ampcal'], field='', calwt=F) # Clean again, stopping clean when the image looks noise-like os.system('rm -rf I16293_clean_continuum_ampcal*') clean(vis='I16293_corrected.ms', imagename='I16293_clean_continuum_ampcal', field='0,1', phasecenter='0', mode='mfs',spw='0:410~430', imagermode='mosaic', imsize=[220,220],cell='0.2arcsec', interactive=T, pbcor=False, weighting='briggs',robust=0.5, niter=10000, threshold='0.5mJy') # RMS noise immediately to south-west of sources has gone down from 18 to 4 mJy/beam ##################################### # Now apply self-cal to line-only data ##################################### # Apply solutions applycal(vis='I16293_corrected.ms.contsub', gaintable=['I16293_corrected.ms.selfcal.pcal2','I16293_corrected.ms.selfcal.ampcal'], field='', calwt=F) # Make this image in frequency space, to get best regridding os.system('rm -rf I16293_corrected.ms.contsub.selfcube*') clean(vis='I16293_corrected.ms.contsub', imagename='I16293_corrected.ms.selfcube', field='0,1',mode='frequency', niter=10000, imagermode='mosaic', phasecenter='0',cell=['0.2arcsec'], imsize=[220,220], nchan=-1, start='220313MHz',width='-488.311kHz', outframe='LSRK',interactive=T, threshold='1.5mJy') ######################################## # Convert to fits format exportfits(imagename='I16293_clean_continuum_ampcal.image', fitsimage='I16293_1.3mm_continuum_ampcal.image.fits') exportfits(imagename='I16293_corrected.ms.selfcube.image', fitsimage='I16293_220GHz_corrected.ms.selfcube.image.fits')