""" This script was written for CASA 5.4.0 Datasets calibrated: SB1: 2017.1.01151.S Observed 04 October 2018 (1 execution block) LB1: 2017.1.01151.S Observed 27 October 2017 (2 execution blocks) """ execfile('reduction_utils.py') #change this path to whatever's appropriate SB1_path = '/path_to/SB1_calibrated_final.ms' LB1_path = '/path_to/LB1_calibrated_final.ms' field = 'GM_Aurigae' prefix = 'GMAurB4' contspws = '5,7,9,11' mask_pa = 53 # position angle of mask in degrees mask_maj = 2 # semimajor axis of mask in arcsec mask_min = 1.5 # semiminor axis of mask in arcsec ellipticalmask = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ ('04h55m10.989s', '30.21.58.95', mask_maj, mask_min, mask_pa) SB_scales = [0, 5, 10, 20,30] LB_scales = [0, 25, 50, 100,175] # optional for visual inspection of data plotms(vis=SB1_path, xaxis='channel', yaxis='amplitude', field=field, ydatacolumn='data', avgtime='1e8', avgscan=True, avgbaseline=True, iteraxis='spw') plotms(vis=LB1_path, xaxis='channel', yaxis='amplitude', field=field, ydatacolumn='data', avgtime='1e8', avgscan=True, avgbaseline=True, iteraxis='spw') # Split out channel-averaged datasets SB1_initcont = prefix+'_SB1_initcont.ms' os.system('rm -rf ' + SB1_initcont + '*') split(vis=SB1_path, field = field, spw=contspws, outputvis=SB1_initcont, timebin = '6s', width=[16,16,16,16], datacolumn='data', keepflags = False) LB1_initcont = prefix+'_LB1_initcont.ms' os.system('rm -rf ' + LB1_initcont + '*') split(vis=LB1_path, field = field, spw=contspws, outputvis=LB1_initcont, width=[16,16,16,16], timebin = '6s', datacolumn='data', keepflags = False) # check visibility amplitudes to see if they're reasonable plotms(vis=SB1_initcont,xaxis='uvdist',yaxis='amp',coloraxis='spw', avgtime = '30', avgchannel = '8') plotms(vis=LB1_initcont,xaxis='uvdist',yaxis='amp',coloraxis='spw', avgtime = '30', avgchannel = '8') #perform initial imaging for inspection SB_initcontimage = prefix+'SB1_initcontinuum' os.system('rm -rf '+SB_initcontimage+'.*') tclean(vis=SB1_initcont, imagename=SB_initcontimage, specmode='mfs', deconvolver = 'multiscale', scales = SB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=1000, cell='0.03arcsec', niter = 50000, interactive = True, threshold = '.05mJy', mask = ellipticalmask, nterms = 1) #Fit Gaussian in CASA to roughly estimate image center fit_gaussian(prefix+'SB1_initcontinuum.image', region=ellipticalmask) #Center from gaussian fit: ICRS 04:55:10.990 +30.21.58.952 LB_initcontimage = prefix+'LB1_initcontinuum' os.system('rm -rf '+LB_initcontimage+'.*') tclean(vis=LB1_initcont, imagename=LB_initcontimage, specmode='mfs', deconvolver = 'multiscale', scales = LB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=2400, cell='0.003arcsec', niter = 50000, interactive = True, threshold = '.03mJy', mask = ellipticalmask, nterms = 1) fit_gaussian(prefix+'LB1_initcontinuum.image', region=ellipticalmask) #Center from gaussian fit: ICRS 04:55:10.986+30.21.58.943 au.ICRSToJ2000('04:55:10.986 +30:21:58.943') #get the J2000 equivalent because fixplanets only works in J2000 coordinates in this version of CASA #define a common direction (center of LB observation) common_dir = 'J2000 04h55m10.986s +030.21.58.932' """ Shift each MS so emission center is at same phase center """ SB1_shift = prefix+'_SB1_initcont_shift.ms' os.system('rm -rf '+SB1_shift+'*') fixvis(vis=SB1_initcont, outputvis=SB1_shift, field=field, phasecenter='ICRS 04h55m10.990s +30d21m58.952s') fixplanets(vis=SB1_shift, field=field, direction=common_dir) LB1_shift = prefix+'_LB1_initcont_shift.ms' os.system('rm -rf '+LB1_shift+'*') fixvis(vis=LB1_initcont, outputvis=LB1_shift, field=field, phasecenter='ICRS 04h55m10.986s +30d21m58.943s') fixplanets(vis=LB1_shift, field=field,direction=common_dir) #Self-cal of SB data SB_cont_p0 = prefix+'SB_contp0' os.system('rm -rf '+SB_cont_p0+'.*') tclean(vis=SB1_shift, imagename=SB_cont_p0, specmode='mfs', deconvolver = 'multiscale', scales = SB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=1000, cell='0.03arcsec', niter = 50000, interactive = True, threshold = '.05mJy', mask = ellipticalmask, savemodel = 'modelcolumn', nterms = 1) #workaround for bug in CASA 5.4 that doesn't automatically fill out the model column under certain circumstances tclean(vis=SB1_shift, imagename=SB_cont_p0, specmode='mfs', deconvolver = 'multiscale', scales = SB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=1000, cell='0.03arcsec', niter = 0, interactive = False, threshold = '.05mJy', mask = '', #interestingly, this fails if the mask is "" rather than '' savemodel = 'modelcolumn', calcres = False, calcpsf = False, nterms = 1) #peak: 6.55 mJy/beam #rms: 0.024 mJy/beam #SNR: 273 SB_contspws = '0~3' SB_refant = 'DV06' #from the pipeline calibration logs SB_p1 = prefix+'_SB.p1' os.system('rm -rf '+SB_p1) gaincal(vis=SB1_shift, caltable=SB_p1, gaintype='T', spw=SB_contspws, refant=SB_refant, calmode='p', solint='30s', minsnr=1.5, minblperant=4, combine = 'spw') plotcal(caltable=SB_p1, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', plotrange=[0,0,-180,180]) applycal(vis=SB1_shift, spw=SB_contspws, gaintable=[SB_p1], spwmap = [0,0,0,0],interp='linearPD', calwt=True) SB_cont_p1 = prefix+'_SB_contp1' os.system('rm -rf %s*' % SB_cont_p1) split(vis=SB1_shift, outputvis=SB_cont_p1+'.ms', datacolumn='corrected') tclean(vis=SB_cont_p1+'.ms', imagename=SB_cont_p1, specmode='mfs', deconvolver = 'multiscale', scales = SB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=1000, cell='0.03arcsec', niter = 50000, interactive = True, threshold = '.05mJy', mask = ellipticalmask, savemodel = 'modelcolumn', nterms = 1) tclean(vis=SB_cont_p1+'.ms', imagename=SB_cont_p1, specmode='mfs', deconvolver = 'multiscale', scales = SB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=1000, cell='0.03arcsec', niter = 0, interactive = False, threshold = '.05mJy', mask = '', savemodel = 'modelcolumn', calcres = False, calcpsf = False, nterms = 1) #CASA throws some kind of error message about the mask after this, but the model plotted in plotms looks correct #peak: 6.88 mJy/beam #rms: 0.0188 mJy/beam #SNR: 366 SB_p2 = prefix+'_SB.p2' os.system('rm -rf '+SB_p2) gaincal(vis=SB_cont_p1+'.ms', caltable=SB_p2, gaintype='T', spw=SB_contspws, refant=SB_refant, calmode='p', solint='20s', minsnr=1.5, minblperant=4, combine = 'spw') plotcal(caltable=SB_p2, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', plotrange=[0,0,-180,180]) applycal(vis=SB_cont_p1+'.ms', spw=SB_contspws, gaintable=[SB_p2], spwmap = [0,0,0,0],interp='linearPD', calwt=True) SB_cont_p2 = prefix+'_SB_contp2' os.system('rm -rf %s*' % SB_cont_p2) split(vis=SB_cont_p1+'.ms', outputvis=SB_cont_p2+'.ms', datacolumn='corrected') tclean(vis=SB_cont_p2+'.ms', imagename=SB_cont_p2, specmode='mfs', deconvolver = 'multiscale', scales = SB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=1000, cell='0.03arcsec', niter = 50000, interactive = False, threshold = '.05mJy', mask = ellipticalmask, savemodel = 'modelcolumn', nterms = 1) tclean(vis=SB_cont_p2+'.ms', imagename=SB_cont_p2, specmode='mfs', deconvolver = 'multiscale', scales = SB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=1000, cell='0.03arcsec', niter = 0, interactive = False, threshold = '.05mJy', savemodel = 'modelcolumn', calcres = False, calcpsf = False, nterms = 1) #CASA throws some kind of error message about the mask after this, but the model plotted in plotms looks correct #peak: 6.89 mJy/beam #rms: 0.0187 mJy/beam #SNR: 368 SB_ap1 = prefix+'_SB.ap1' os.system('rm -rf '+SB_ap1) gaincal(vis=SB_cont_p2+'.ms', caltable=SB_ap1, gaintype='T', combine = 'spw', spw=SB_contspws, refant=SB_refant, calmode='ap', solint='inf', minsnr=3.0, minblperant=4, solnorm = True) plotcal(caltable=SB_ap1, xaxis='time', yaxis='amp',subplot=441, iteration='antenna', plotrange = [0,0,0.5, 1.5]) applycal(vis=SB_cont_p2+'.ms', spw=SB_contspws, gaintable=[SB_ap1], spwmap = [0,0,0,0],interp='linearPD', calwt=True) SB_cont_ap = prefix+'_SB_contap' os.system('rm -rf %s*' % SB_cont_ap) split(vis=SB_cont_p2+'.ms', outputvis=SB_cont_ap+'.ms', datacolumn='corrected') tclean(vis=SB_cont_ap+'.ms', imagename=SB_cont_ap, specmode='mfs', deconvolver = 'multiscale', scales = SB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=1000, cell='0.03arcsec', niter = 50000, interactive = True, threshold = '.05mJy', mask = ellipticalmask, nterms = 1) tclean(vis=SB_cont_ap+'.ms', imagename=SB_cont_ap, specmode='mfs', deconvolver = 'multiscale', scales = SB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=1000, cell='0.03arcsec', niter = 0, interactive = False, threshold = '.05mJy', savemodel = 'modelcolumn', calcres = False, calcpsf = False, nterms = 1) #marginal change #combining short and long baseline data concat(vis = [SB_cont_ap+'.ms', LB1_shift], concatvis = prefix+'_contcombined.ms', dirtol = '0.1arcsec', copypointing = False) #start of self-cal for combined data combined_cont_p0 = prefix+'combined_contp0' os.system('rm -rf '+combined_cont_p0+'.*') tclean(vis=prefix+'_contcombined.ms', imagename=combined_cont_p0, specmode='mfs', deconvolver = 'multiscale', scales = LB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=2400, cell='0.003arcsec', niter = 50000, interactive = True, threshold = '.03mJy', mask = ellipticalmask, savemodel = 'modelcolumn', nterms = 1) tclean(vis=prefix+'_contcombined.ms', imagename=combined_cont_p0, specmode='mfs', deconvolver = 'multiscale', scales = LB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=2400, cell='0.003arcsec', niter = 0, interactive = False, threshold = '.03mJy', savemodel = 'modelcolumn', calcres = False, calcpsf = False, nterms = 1) #peak: 0.61 mJy/beam #rms: 9.6 microJy/beam #SNR: 64 combined_spw = '0~3' combined_refant = 'DV25@A087, DA45@A135, DV06@A027' combined_p1 = prefix+'_combined.p1' os.system('rm -rf '+combined_p1) gaincal(vis=prefix+'_contcombined.ms', caltable=combined_p1, gaintype='T', combine = 'spw,scan', spw='0~3', refant=combined_refant, calmode='p', solint='300s', minsnr=1.5, minblperant=4) applycal(vis=prefix+'_contcombined.ms', spw='0~3', spwmap = [0]*4, gaintable=[combined_p1], calwt=True, applymode = 'calonly', flagbackup=False) combined_cont_p1 = prefix+'_combined_contp1' os.system('rm -rf '+combined_cont_p1+'.*') split(vis=prefix+'_contcombined.ms', outputvis=combined_cont_p1+'.ms', datacolumn='corrected', keepflags = False) tclean(vis=combined_cont_p1+'.ms', imagename=combined_cont_p1, specmode='mfs', deconvolver = 'multiscale', scales = LB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=2400, cell='0.003arcsec', niter = 50000, interactive = True, threshold = '.03mJy', mask = ellipticalmask, savemodel = 'modelcolumn', nterms = 1) tclean(vis=combined_cont_p1+'.ms', imagename=combined_cont_p1, specmode='mfs', deconvolver = 'multiscale', scales = LB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=2400, cell='0.003arcsec', niter = 0, interactive = False, threshold = '.03mJy', savemodel = 'modelcolumn', calcres = False, calcpsf = False, nterms = 1) #peak: 0.63 mJy/beam #rms: 9.2 microJy/beam #SNR: 68 combined_p2 = prefix+'_combined.p2' os.system('rm -rf '+combined_p2) gaincal(vis=combined_cont_p1+'.ms', caltable=combined_p2, gaintype='T', combine = 'spw,scan', spw='0~3', refant=combined_refant, calmode='p', solint='150s', minsnr=1.5, minblperant=4) applycal(vis=combined_cont_p1+'.ms', spw='0~3', spwmap = [0]*4, gaintable=[combined_p2], calwt=True, applymode = 'calonly', flagbackup=False) combined_cont_p2 = prefix+'_combined_contp2' os.system('rm -rf '+combined_cont_p2+'.*') split(vis=combined_cont_p1+'.ms', outputvis=combined_cont_p2+'.ms', datacolumn='corrected', keepflags = False) tclean(vis=combined_cont_p2+'.ms', imagename=combined_cont_p2, specmode='mfs', deconvolver = 'multiscale', scales = LB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=2400, cell='0.003arcsec', niter = 50000, interactive = True, threshold = '.03mJy', mask = ellipticalmask, savemodel = 'modelcolumn', nterms = 1) tclean(vis=combined_cont_p2+'.ms', imagename=combined_cont_p2, specmode='mfs', deconvolver = 'multiscale', scales = LB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=2400, cell='0.003arcsec', niter = 0, interactive = False, threshold = '.03mJy', savemodel = 'modelcolumn', calcres = False, calcpsf = False, nterms = 1) #peak: 0.66 mJy/beam #rms: 9.2 microJy/beam #SNR: 72 combined_p3 = prefix+'_combined.p3' os.system('rm -rf '+combined_p3) gaincal(vis=combined_cont_p2+'.ms', caltable=combined_p3, gaintype='T', combine = 'spw,scan', spw='0~3', refant=combined_refant, calmode='p', solint='60s', minsnr=1.5, minblperant=4) applycal(vis=combined_cont_p2+'.ms', spw='0~3', spwmap = [0]*4, gaintable=[combined_p3], calwt=True, applymode = 'calonly', flagbackup=False) combined_cont_p3 = prefix+'_combined_contp3' os.system('rm -rf '+combined_cont_p3+'.*') split(vis=combined_cont_p2+'.ms', outputvis=combined_cont_p3+'.ms', datacolumn='corrected', keepflags = False) tclean(vis=combined_cont_p3+'.ms', imagename=combined_cont_p3, specmode='mfs', deconvolver = 'multiscale', scales = LB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=2400, cell='0.003arcsec', niter = 50000, interactive = True, threshold = '.03mJy', mask = ellipticalmask, savemodel = 'modelcolumn', nterms = 1) tclean(vis=combined_cont_p3+'.ms', imagename=combined_cont_p3, specmode='mfs', deconvolver = 'multiscale', scales = LB_scales, weighting='briggs', robust=0.5, gain = 0.3, imsize=2400, cell='0.003arcsec', niter = 0, interactive = False, threshold = '.03mJy', savemodel = 'modelcolumn', calcres = False, calcpsf = False, nterms = 1) #not changed much from previous round #ampcal makes things worse, so we stop here #imaging with lower robust parameter to make synthesized beam more similar to B6 images tclean(vis=combined_cont_p3+'.ms', imagename=prefix+'_combinedcont_robust0', specmode='mfs', deconvolver = 'multiscale', scales = LB_scales, weighting='briggs', robust=0, gain = 0.1, imsize=2400, cell='0.003arcsec', niter = 50000, interactive = True, threshold = '.04mJy', mask = ellipticalmask, savemodel = 'modelcolumn', nterms = 1) #label final data products os.system('mv '+combined_cont_p3+'.ms GMAurB4continuumfinal.ms') exportfits(imagename=prefix+'_combinedcont_robust0.image', fitsimage = 'GMAurB4continuumfinal.image.fits')