""" This script was written for CASA 5.1.1 Datasets calibrated (in order of date observed): SB1: 2013.1.00226.S Observed 02 July 2014 (1 execution block) SB2: 2013.1.00226.S Observed 17 July 2014 (1 execution block) SB3: 2015.1.00486.S Observed 22 September 2016 (1 execution block) SB4: 2015.1.00486.S Observed 26 September 2016 (1 execution block) SB5: 2016.1.00484.L Observed 09 May 2017 (1 execution block) LB1: 2016.1.00484.L Observed 07 September 2017 and 19 September 2017 (2 execution blocks) reducer: V. Guzman """ """ Starting matter """ import os execfile('reduction_utils.py') skip_plots = True # if True, can run script non-interactively """ Input for loading data """ prefix = 'AS209' SB1_path = '/full_path/to_calibrated/msfile.ms' # Note we only use the SB1 upper sideband, with frequencies near the LP data SB2_path = '/full_path/to_calibrated/msfile.ms' # Note we only use the SB2 lower sideband, with frequencies near the LP data SB3_path = '/full_path/to_calibrated/msfile.ms' # Note we only use the SB3 upper sideband, with frequencies near the LP data SB4_path = '/full_path/to_calibrated/msfile.ms' # Note we only use the SB4 upper sideband, with frequencies near the LP data SB5_path = '/full_path/to_calibrated/msfile.ms' LB1_path = '/full_path/to_calibrated/msfile.ms' # Note that if you are downloading data from the archive, your SPW numbering # may differ from this script, depending on how you split your data out! data_params = {'SB1': {'vis' : SB1_path, 'name' : 'SB1', 'field': 'as_209', 'line_spws': np.array([15]), # CO SPWs 'line_freqs': np.array([2.30538e11]), 'cont_spws': '14,15,16', 'width_array': [960,960,960], }, 'SB2': {'vis' : SB2_path, 'name' : 'SB2', 'field': 'as_209', 'line_spws': np.array([]), # CO SPWs 'line_freqs': np.array([]), 'cont_spws': '12,13,14,15,16,17', 'width_array': [1920,1920,960,960,960,960], }, 'SB3': {'vis' : SB3_path, 'name' : 'SB3', 'field': 'AS_209', 'line_spws': np.array([1]), # CO SPWs 'line_freqs': np.array([2.30538e11]), 'cont_spws': '0,1', 'width_array': [8,480], }, 'SB4': {'vis' : SB4_path, 'name' : 'SB4', 'field': 'AS_209', 'line_spws': np.array([1]), # CO SPWs 'line_freqs': np.array([2.30538e11]), 'cont_spws': '0,1', 'width_array': [8,480], }, 'SB5': {'vis' : SB5_path, 'name' : 'SB5', 'field': 'AS_209', 'line_spws': np.array([0]), # CO SPWs 'line_freqs': np.array([2.30538e11]), 'cont_spws': None, 'width_array': None, }, 'LB1': {'vis' : LB1_path, 'name' : 'LB1', 'field' : 'AS_209', 'line_spws': np.array([3,7]), # CO SPWs 'line_freqs': np.array([2.30538e11, 2.30538e11]), 'cont_spws': None, 'width_array': None, } } """ Check data (various options here; an example) """ if not skip_plots: for i in data_params.keys(): plotms(vis=data_params[i]['vis'], xaxis='channel', yaxis='amplitude', field=data_params[i]['field'], ydatacolumn='data', avgtime='1e8', avgscan=True, avgbaseline=True, iteraxis='spw') """ Identify region containing CO emission; then flag that and do a spectral average to a pseudo-continuum MS """ for i in data_params.keys(): flagchannels_string = get_flagchannels(data_params[i], prefix, velocity_range=np.array([-10, 20])) avg_cont(data_params[i], prefix, flagchannels=flagchannels_string, contspws=data_params[i]['cont_spws'], width_array=data_params[i]['width_array']) # Flagchannels input string for SB1: '15:461~839' # No CO in SB2 # Flagchannels input string for SB3: '1:136~325' # Flagchannels input string for SB4: '1:136~325' # Flagchannels input string for SB5: '0:1874~1968' # Flagchannels input string for LB1: '3:1874~1968, 7:1874~1968' """ Define simple masks and clean scales for imaging """ mask_pa = 86 # position angle of mask in degrees mask_maj = 1.3 # semimajor axis of mask in arcsec mask_min = 1.1 # semiminor axis of mask in arcsec mask_ra = '16h49m15.29s' mask_dec = '-14.22.09.04' SB1_mask = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ (mask_ra, mask_dec, mask_maj, mask_min, mask_pa) LB1_mask = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ (mask_ra, mask_dec, mask_maj, mask_min, mask_pa) SB_scales = [0, 5, 10, 20] LB_scales = [0, 5, 30, 100, 200] if not skip_plots: """ Image each dataset individually """ # images are saved in the format prefix+'_name_initcont_exec#.ms' image_each_obs(data_params['SB1'], prefix, mask=SB1_mask, scales=SB_scales, threshold='1.0mJy', interactive=False) image_each_obs(data_params['SB2'], prefix, mask=SB1_mask, scales=SB_scales, threshold='0.5mJy', interactive=False) image_each_obs(data_params['SB3'], prefix, mask=SB1_mask, scales=SB_scales, threshold='1.0mJy', interactive=False) image_each_obs(data_params['SB4'], prefix, mask=SB1_mask, scales=SB_scales, threshold='1.0mJy', interactive=False) image_each_obs(data_params['SB5'], prefix, mask=SB1_mask, scales=SB_scales, threshold='0.25mJy', interactive=False) image_each_obs(data_params['LB1'], prefix, mask=LB1_mask, scales=LB_scales, threshold='0.06mJy', interactive=False) """ Fit Gaussians to roughly estimate centers, inclinations, PAs """ fit_gaussian(prefix+'_SB1_initcont_exec0.image', region=SB1_mask) #Peak : J2000 16h49m15.299722s -14d22m08.85131s fit_gaussian(prefix+'_SB2_initcont_exec0.image', region=SB1_mask) #Peak : J2000 16h49m15.301179s -14d22m08.95623s fit_gaussian(prefix+'_SB3_initcont_exec0.image', region=SB1_mask) #Peak : ICRS 16h49m15.292292s -14d22m09.02254s fit_gaussian(prefix+'_SB4_initcont_exec0.image', region=SB1_mask) #Peak : J2000 16:49:15.29443, -014:22:09.070855 fit_gaussian(prefix+'_SB5_initcont_exec0.image', region=SB1_mask) #Peak : J2000 16:49:15.29399, -014:22:09.050565 fit_gaussian(prefix+'_LB1_initcont_exec0.image', region=LB1_mask) #Peak : J2000 16:49:15.29425, -014:22:09.037705 fit_gaussian(prefix+'_LB1_initcont_exec1.image', region=LB1_mask) #Peak : J2000 16:49:15.29463, -014:22:09.048165 #PA of Gaussian component: 85.41 deg #Inclination of Gaussian component: 34.93 deg """ The emission centers are slightly misaligned. So we split out the individual executions, shift the peaks to the phase center, and reassign the phase centers to a common direction. (Not in script: check that these shifts do what you think they should by re-imaging! They do for us.) """ """ Split out individual MSs for each execution """ split_all_obs(prefix+'_LB1_initcont.ms', prefix+'_LB1_initcont_exec') """ Define a common direction (peak of 2nd LB EB) """ common_dir = 'J2000 16h49m15.29463s -014.22.09.048165' """ 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=prefix+'_SB1_initcont.ms', outputvis=SB1_shift, field=data_params['SB1']['field'], phasecenter='J2000 16h49m15.299722s -14d22m08.85131s') fixplanets(vis=SB1_shift, field=data_params['SB1']['field'], direction=common_dir) SB2_shift = prefix+'_SB2_initcont_shift.ms' os.system('rm -rf '+SB2_shift+'*') fixvis(vis=prefix+'_SB2_initcont.ms', outputvis=SB2_shift, field=data_params['SB2']['field'], phasecenter='J2000 16h49m15.301179s -14d22m08.95623s') fixplanets(vis=SB2_shift, field=data_params['SB2']['field'], direction=common_dir) SB3_shift = prefix+'_SB3_initcont_shift.ms' os.system('rm -rf '+SB3_shift+'*') fixvis(vis=prefix+'_SB3_initcont.ms', outputvis=SB3_shift, field=data_params['SB3']['field'], phasecenter='ICRS 16h49m15.292292s -14d22m09.02254s') fixplanets(vis=SB3_shift, field=data_params['SB3']['field'], direction=common_dir) SB4_shift = prefix+'_SB4_initcont_shift.ms' os.system('rm -rf '+SB4_shift+'*') fixvis(vis=prefix+'_SB4_initcont.ms', outputvis=SB4_shift, field=data_params['SB4']['field'], phasecenter='ICRS 16h49m15.293687s -14d22m09.08240s') fixplanets(vis=SB4_shift, field=data_params['SB4']['field'], direction=common_dir) SB5_shift = prefix+'_SB5_initcont_shift.ms' os.system('rm -rf '+SB5_shift+'*') fixvis(vis=prefix+'_SB5_initcont.ms', outputvis=SB5_shift, field=data_params['SB5']['field'], phasecenter='ICRS 16h49m15.293254s -14d22m09.06211s') fixplanets(vis=SB5_shift, field=data_params['SB5']['field'], direction=common_dir) LB0_shift = prefix+'_LB1_initcont_exec0_shift.ms' os.system('rm -rf '+LB0_shift+'*') fixvis(vis=prefix+'_LB1_initcont_exec0.ms', outputvis=LB0_shift, field=data_params['LB1']['field'], phasecenter='ICRS 16h49m15.293513s -14d22m09.04925s') fixplanets(vis=LB0_shift, field=data_params['LB1']['field'], direction=common_dir) LB1_shift = prefix+'_LB1_initcont_exec1_shift.ms' os.system('rm -rf '+LB1_shift+'*') fixvis(vis=prefix+'_LB1_initcont_exec1.ms', outputvis=LB1_shift, field=data_params['LB1']['field'], phasecenter='ICRS 16h49m15.293891s -14d22m09.05974s') fixplanets(vis=LB1_shift, field=data_params['LB1']['field'], direction=common_dir) """ Now that everything is aligned, we inspect the flux calibration. """ if not skip_plots: """ Assign rough emission geometry parameters. """ PA, incl = 86, 35 """ Export MS contents into Numpy save files """ for msfile in [prefix+'_SB1_initcont_shift.ms', prefix+'_SB2_initcont_shift.ms', prefix+'_SB3_initcont_shift.ms', prefix+'_SB4_initcont_shift.ms', prefix+'_SB5_initcont_shift.ms', prefix+'_LB1_initcont_exec0_shift.ms', prefix+'_LB1_initcont_exec1_shift.ms']: export_MS(msfile) """ Plot deprojected visibility profiles for all data together """ plot_deprojected([prefix+'_SB1_initcont_shift.vis.npz', prefix+'_SB2_initcont_shift.vis.npz', prefix+'_SB3_initcont_shift.vis.npz', prefix+'_SB4_initcont_shift.vis.npz', prefix+'_SB5_initcont_shift.vis.npz', prefix+'_LB1_initcont_exec0_shift.vis.npz', prefix+'_LB1_initcont_exec1_shift.vis.npz'], fluxscale=[1.,1.,1.,1.,1.,1.,1.], PA=PA, incl=incl, show_err=False) """ Now inspect offsets by comparing against a reference """ estimate_flux_scale(reference=prefix+'_LB1_initcont_exec0_shift.vis.npz', comparison=prefix+'_SB3_initcont_shift.vis.npz', incl=incl, PA=PA) estimate_flux_scale(reference=prefix+'_LB1_initcont_exec0_shift.vis.npz', comparison=prefix+'_SB4_initcont_shift.vis.npz', incl=incl, PA=PA) """ Correct the flux scales where appropriate. """ """ Here we did a single iteration of phase-only self-cal to make sure that phase noise was not the issue """ rescale_flux(prefix+'_SB1_initcont_shift.ms', [0.919]) rescale_flux(prefix+'_SB2_initcont_shift.ms', [1.072]) rescale_flux(prefix+'_SB4_initcont_shift.ms', [0.965]) """ SELF-CAL for short-baseline data """ """ Merge the SB executions back into a single MS """ SB_cont_p0 = prefix+'_SB_contp0' os.system('rm -rf %s*' % SB_cont_p0) concat(vis=[prefix+'_SB1_initcont_shift_rescaled.ms', prefix+'_SB2_initcont_shift_rescaled.ms', prefix+'_SB3_initcont_shift.ms', prefix+'_SB4_initcont_shift_rescaled.ms', prefix+'_SB5_initcont_shift.ms'], concatvis=SB_cont_p0+'.ms', dirtol='0.1arcsec', copypointing=False) """ Set up a clean mask """ mask_ra = '16h49m15.29463s' mask_dec = '-14.22.09.048165' common_mask = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ (mask_ra, mask_dec, mask_maj, mask_min, mask_pa) """ Initial clean """ tclean_wrapper(vis=SB_cont_p0+'.ms', imagename=SB_cont_p0, mask=common_mask, scales=SB_scales, threshold='0.2mJy', savemodel='modelcolumn') """ Define a noise annulus, measure the peak SNR in map """ noise_annulus = "annulus[[%s, %s],['%.2farcsec', '4.25arcsec']]" % \ (mask_ra, mask_dec, 1.1*mask_maj) estimate_SNR(SB_cont_p0+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #AS209_SB_contp0.image #Beam 0.216 arcsec x 0.190 arcsec (-86.70 deg) #Flux inside disk mask: 246.34 mJy #Peak intensity of source: 18.33 mJy/beam #rms: 1.20e-01 mJy/beam #Peak SNR: 153.13 """ Self-calibration parameters """ SB_contspws = '0~16' SB_refant = 'DA48@A046,DV22@A011,DA41@A004' SB1_timerange = '2014/07/02/00~2014/07/03/00' SB2_timerange = '2014/07/17/00~2014/07/18/00' SB3_timerange = '2016/09/22/00~2016/09/24/00' SB4_timerange = '2016/09/26/00~2016/09/28/00' SB5_timerange = '2017/05/09/00~2017/05/10/00' """ First round of phase-only self-cal (short baselines only) """ # Use scan-length interval to align SPWs as much as possible SB_p1 = prefix+'_SB.p1' os.system('rm -rf '+SB_p1) gaincal(vis=SB_cont_p0+'.ms', caltable=SB_p1, gaintype='T', spw=SB_contspws, refant=SB_refant, calmode='p', solint='inf', minsnr=1.5, minblperant=4) if not skip_plots: """ Inspect gain tables """ plotcal(caltable=SB_p1, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB1_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p1, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB2_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p1, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB3_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p1, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB4_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p1, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB5_timerange, plotrange=[0,0,-180,180]) """ Apply the solutions """ applycal(vis=SB_cont_p0+'.ms', spw=SB_contspws, gaintable=[SB_p1], interp='linearPD', calwt=True) """ Split off a corrected MS """ SB_cont_p1 = prefix+'_SB_contp1' os.system('rm -rf %s*' % SB_cont_p1) split(vis=SB_cont_p0+'.ms', outputvis=SB_cont_p1+'.ms', datacolumn='corrected') """ Image the results; check the resulting map """ tclean_wrapper(vis=SB_cont_p1+'.ms', imagename=SB_cont_p1, mask=common_mask, scales=SB_scales, threshold='0.1mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p1+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #AS209_SB_contp1.image #Beam 0.219 arcsec x 0.191 arcsec (-87.60 deg) #Flux inside disk mask: 258.83 mJy #Peak intensity of source: 19.74 mJy/beam #rms: 3.48e-02 mJy/beam #Peak SNR: 566.72 """ Second round of phase-only self-cal (short baselines only) """ 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='60s', minsnr=1.5, minblperant=4) if not skip_plots: """ Inspect gain tables """ plotcal(caltable=SB_p2, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB1_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p2, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB2_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p2, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB3_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p2, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB4_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p2, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB5_timerange, plotrange=[0,0,-180,180]) """ Apply the solutions """ applycal(vis=SB_cont_p1+'.ms', spw=SB_contspws, gaintable=[SB_p2], interp='linearPD', calwt=True) """ Split off a corrected MS """ 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') """ Image the results; check the resulting map """ tclean_wrapper(vis=SB_cont_p2+'.ms', imagename=SB_cont_p2, mask=common_mask, scales=SB_scales, threshold='0.1mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p2+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #AS209_SB_contp2.image #Beam 0.224 arcsec x 0.193 arcsec (-89.21 deg) #Flux inside disk mask: 260.45 mJy #Peak intensity of source: 20.88 mJy/beam #rms: 3.06e-02 mJy/beam #Peak SNR: 682.66 # Additional phase-only self-cal with shorter intervals ends up flagging a lot # of the data due to noisy individual SPWs. This will get cleaned up in the # combined self-cal when the SPWs are averaged. """ Amplitude self-cal (short baselines only) """ SB_ap = prefix+'_SB.ap' os.system('rm -rf '+SB_ap) gaincal(vis=SB_cont_p2+'.ms', caltable=SB_ap, gaintype='T', spw=SB_contspws, refant=SB_refant, calmode='ap', solint='inf', minsnr=3.0, minblperant=4, solnorm=False) if not skip_plots: """ Inspect gain tables """ plotcal(caltable=SB_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=SB1_timerange, plotrange=[0,0,0,2]) plotcal(caltable=SB_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=SB2_timerange, plotrange=[0,0,0,2]) plotcal(caltable=SB_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=SB3_timerange, plotrange=[0,0,0,2]) plotcal(caltable=SB_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=SB4_timerange, plotrange=[0,0,0,2]) plotcal(caltable=SB_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=SB5_timerange, plotrange=[0,0,0,2]) """ Apply the solutions """ applycal(vis=SB_cont_p2+'.ms', spw=SB_contspws, gaintable=[SB_ap], interp='linearPD', calwt=True) """ Split off a corrected MS """ 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') """ Image the results; check the resulting map """ tclean_wrapper(vis=SB_cont_ap+'.ms', imagename=SB_cont_ap, mask=common_mask, scales=SB_scales, threshold='0.1mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_ap+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #AS209_SB_contap.image #Beam 0.219 arcsec x 0.189 arcsec (-89.80 deg) #Flux inside disk mask: 259.63 mJy #Peak intensity of source: 20.14 mJy/beam #rms: 2.57e-02 mJy/beam #Peak SNR: 784.05 """ SELF-CAL for the combined (short-baseline + long-baseline) data """ """ Merge the SB+LB executions into a single MS """ combined_cont_p0 = prefix+'_combined_contp0' os.system('rm -rf %s*' % combined_cont_p0) concat(vis=[SB_cont_ap+'.ms', prefix+'_LB1_initcont_exec0_shift.ms', prefix+'_LB1_initcont_exec1_shift.ms'], concatvis=combined_cont_p0+'.ms', dirtol='0.1arcsec', copypointing=False) """ Initial clean """ tclean_wrapper(vis=combined_cont_p0+'.ms', imagename=combined_cont_p0, mask=common_mask, scales=LB_scales, threshold='0.06mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p0+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #AS209_combined_contp0.image #Beam 0.056 arcsec x 0.039 arcsec (-86.32 deg) #Flux inside disk mask: 256.32 mJy #Peak intensity of source: 2.30 mJy/beam #rms: 1.26e-02 mJy/beam #Peak SNR: 181.62 """ Self-calibration parameters """ combined_contspws = '0~24' combined_refant = 'DA61@A015,DA48@A046,DV22@A017,DA41@A004' combined_spwmap = [0,0,0,3,3,3,3,3,3,9,9,11,11,13,13,13,13,17,17,17,17,21,21,21,21] LB1_obs0_timerange = '2017/09/06/00~2017/09/09/00' LB2_obs1_timerange = '2017/09/18/00~2017/09/22/00' """ First round of phase-only self-cal (all data) """ combined_p1 = prefix+'_combined.p1' os.system('rm -rf '+combined_p1) gaincal(vis=combined_cont_p0+'.ms', caltable=combined_p1, gaintype='T', combine='spw,scan', spw=combined_contspws, refant=combined_refant, calmode='p', solint='900s', minsnr=1.5, minblperant=4) if not skip_plots: """ Inspect gain tables """ plotcal(caltable=combined_p1, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=LB1_obs0_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=combined_p1, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=LB2_obs1_timerange, plotrange=[0,0,-180,180]) """ Apply the solutions """ applycal(vis=combined_cont_p0+'.ms', spw=combined_contspws, spwmap=combined_spwmap, gaintable=[combined_p1], interp='linearPD', calwt=True, applymode='calonly') """ Split off a corrected MS """ combined_cont_p1 = prefix+'_combined_contp1' os.system('rm -rf %s*' % combined_cont_p1) split(vis=combined_cont_p0+'.ms', outputvis=combined_cont_p1+'.ms', datacolumn='corrected') """ Image the results; check the resulting map """ tclean_wrapper(vis=combined_cont_p1+'.ms', imagename=combined_cont_p1, mask=common_mask, scales=LB_scales, threshold='0.06mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p1+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #AS209_combined_contp1.image #Beam 0.056 arcsec x 0.039 arcsec (-86.32 deg) #Flux inside disk mask: 256.33 mJy #Peak intensity of source: 2.38 mJy/beam #rms: 1.18e-02 mJy/beam #Peak SNR: 201.03 """ Second round of phase-only self-cal (all data) """ 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=combined_contspws, refant=combined_refant, calmode='p', solint='360s', minsnr=1.5, minblperant=4) if not skip_plots: """ Inspect gain tables """ plotcal(caltable=combined_p2, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=LB1_obs0_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=combined_p2, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=LB2_obs1_timerange, plotrange=[0,0,-180,180]) """ Apply the solutions """ applycal(vis=combined_cont_p1+'.ms', spw=combined_contspws, spwmap=combined_spwmap, gaintable=[combined_p2], interp='linearPD', calwt=True, applymode='calonly') """ Split off a corrected MS """ combined_cont_p2 = prefix+'_combined_contp2' os.system('rm -rf %s*' % combined_cont_p2) split(vis=combined_cont_p1+'.ms', outputvis=combined_cont_p2+'.ms', datacolumn='corrected') """ Image the results; check the resulting map """ tclean_wrapper(vis=combined_cont_p2+'.ms', imagename=combined_cont_p2, mask=common_mask, scales=LB_scales, threshold='0.04mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p2+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #AS209_combined_contp2.image #Beam 0.056 arcsec x 0.039 arcsec (-86.32 deg) #Flux inside disk mask: 257.83 mJy #Peak intensity of source: 2.45 mJy/beam #rms: 1.12e-02 mJy/beam #Peak SNR: 219.57 """ Third round of phase-only self-cal (all data) """ 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=combined_contspws, refant=combined_refant, calmode='p', solint='180s', minsnr=1.5, minblperant=4) if not skip_plots: """ Inspect gain tables """ plotcal(caltable=combined_p3, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=LB1_obs0_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=combined_p3, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=LB2_obs1_timerange, plotrange=[0,0,-180,180]) """ Apply the solutions """ applycal(vis=combined_cont_p2+'.ms', spw=combined_contspws, spwmap=combined_spwmap, gaintable=[combined_p3], interp='linearPD', calwt=True, applymode='calonly') """ Split off a corrected MS """ combined_cont_p3 = prefix+'_combined_contp3' os.system('rm -rf %s*' % combined_cont_p3) split(vis=combined_cont_p2+'.ms', outputvis=combined_cont_p3+'.ms', datacolumn='corrected') """ Image the results; check the resulting map """ tclean_wrapper(vis=combined_cont_p3+'.ms', imagename=combined_cont_p3, mask=common_mask, scales=LB_scales, threshold='0.04mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p3+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #AS209_combined_contp3.image #Beam 0.056 arcsec x 0.039 arcsec (-86.32 deg) #Flux inside disk mask: 257.62 mJy #Peak intensity of source: 2.48 mJy/beam #rms: 1.11e-02 mJy/beam #Peak SNR: 224.30 # *** map quality still improving *** """ Fourth round of phase-only self-cal (all data) """ combined_p4 = prefix+'_combined.p4' os.system('rm -rf '+combined_p4) gaincal(vis=combined_cont_p3+'.ms', caltable=combined_p4, gaintype='T', combine='spw,scan', spw=combined_contspws, refant=combined_refant, calmode='p', solint='60s', minsnr=1.5, minblperant=4) if not skip_plots: """ Inspect gain tables """ plotcal(caltable=combined_p4, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=LB1_obs0_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=combined_p4, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=LB2_obs1_timerange, plotrange=[0,0,-180,180]) """ Apply the solutions """ applycal(vis=combined_cont_p3+'.ms', spw=combined_contspws, spwmap=combined_spwmap, gaintable=[combined_p4], interp='linearPD', calwt=True, applymode='calonly') """ Split off a corrected MS """ combined_cont_p4 = prefix+'_combined_contp4' os.system('rm -rf %s*' % combined_cont_p4) split(vis=combined_cont_p3+'.ms', outputvis=combined_cont_p4+'.ms', datacolumn='corrected') """ Image the results; check the resulting map """ tclean_wrapper(vis=combined_cont_p4+'.ms', imagename=combined_cont_p4, mask=common_mask, scales=LB_scales, threshold='0.04mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p4+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #AS209_combined_contp4.image #Beam 0.056 arcsec x 0.039 arcsec (-86.32 deg) #Flux inside disk mask: 258.02 mJy #Peak intensity of source: 2.55 mJy/beam #rms: 1.10e-02 mJy/beam #Peak SNR: 232.09 # Phase-only cal on shorter intervals or amp self-cal are not helpful. """ Final outputs """ """ Save the final MS """ os.system('rm -r '+prefix+'_continuum.ms') split(vis=combined_cont_p4+'.ms', outputvis=prefix+'_continuum.ms', datacolumn='data') os.system('tar cvzf '+prefix+'_continuum.ms.tgz '+prefix+'_continuum.ms') """ Make a fiducial continuum image (based on experimentation) """ tclean_wrapper(vis=combined_cont_p4+'.ms', imagename=prefix+'_continuum', mask=common_mask, scales=LB_scales, robust=-0.5, gain=0.2, threshold='.08mJy', uvtaper=['.037arcsec','.01arcsec','162deg']) exportfits(prefix+'_continuum.image', prefix+'_continuum.fits')