""" This script was written for CASA 5.1.1 Datasets calibrated (in order of date observed): SB1: 2013.1.00498.S Observed 21 July 2015 (1 execution block) LB1: 2016.1.00484.L Observed 25 September 2017 and 09 November 2017 (2 execution blocks) reducer: S. Andrews """ """ Starting matter """ import os execfile('reduction_utils.py') skip_plots = True # if True, can run script non-interactively """ Input for loading data """ prefix = 'HD142666' SB1_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': 'HD_142666', 'line_spws': np.array([1,5,6]), # CO 13CO, C18O SPWs 'line_freqs': np.array([2.30538e11, 2.2039868420e11, 2.1956035410e11]), }, 'LB1': {'vis' : LB1_path, 'name' : 'LB1', 'field' : 'HD_142666', 'line_spws': np.array([3,7]), # CO SPWs 'line_freqs': np.array([2.30538e11, 2.30538e11]), } } """ 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 50 km/s-wide 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([-15, 25])) avg_cont(data_params[i], prefix, flagchannels=flagchannels_string) """ Define simple masks and clean scales for imaging """ mask_pa = 164 # position angle of mask in degrees mask_maj = 0.6 # semimajor axis of mask in arcsec mask_min = 0.4 # semiminor axis of mask in arcsec SB1_mask = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ ('15h56m40.00s', '-22.01.40.40', mask_maj, mask_min, mask_pa) LB1_mask = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ ('15h56m40.00s', '-22.01.40.40', mask_maj, mask_min, mask_pa) SB_scales = [0, 5, 10, 20] LB_scales = [0, 5, 20, 60, 180] 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='0.2mJy', interactive=False) image_each_obs(data_params['LB1'], prefix, mask=LB1_mask, scales=LB_scales, threshold='0.1mJy', interactive=False) """ Fit Gaussians to roughly estimate centers, inclinations, PAs """ fit_gaussian(prefix+'_SB1_initcont_exec0.image', region=SB1_mask) #Peak : J2000 15h56m40.005498s -22d01m40.38913s #PA of Gaussian component: 161.38 deg #Inclination of Gaussian component: 58.49 deg fit_gaussian(prefix+'_LB1_initcont_exec0.image', region=LB1_mask) #Peak : ICRS 15h56m40.004870s -22d01m40.39670s #PA of Gaussian component: 162.35 deg #Inclination of Gaussian component: 61.70 deg fit_gaussian(prefix+'_LB1_initcont_exec1.image', region=LB1_mask) #Peak : ICRS 15h56m40.005249s -22d01m40.39771s #PA of Gaussian component: 163.32 deg #Inclination of Gaussian component: 61.95 deg """ The emission appears to be aligned. But, we want to force everything to the same phase center to avoid confusion in the visibilities. """ """ Split out individual MSs for each execution """ split_all_obs(prefix+'_SB1_initcont.ms', prefix+'_SB1_initcont_exec') split_all_obs(prefix+'_LB1_initcont.ms', prefix+'_LB1_initcont_exec') """ Shift each MS to same phase center (LB1 EB1 as reference) """ SB0_shift = prefix+'_SB1_initcont_exec0_shift.ms' os.system('rm -rf '+SB0_shift+'*') fixvis(vis=prefix+'_SB1_initcont_exec0.ms', outputvis=SB0_shift, field=data_params['SB1']['field'], phasecenter='ICRS 15h56m40.008076s -22d01m40.42728s') 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 15h56m40.008076s -22d01m40.42728s') """ Now that everything is aligned, we inspect the flux calibration. """ if not skip_plots: """ Assign rough emission geometry parameters and centers. """ PA, incl = 163., 62. phasecenter = au.radec2deg('15:56:40.008076, -22.01.40.42728') peakpos = au.radec2deg('15:56:40.005, -22.01.40.398') offsets = au.angularSeparation(peakpos[0], peakpos[1], phasecenter[0], phasecenter[1], True) offx, offy = 3600.*offsets[3], 3600.*offsets[2] """ Export MS contents into Numpy save files """ for msfile in [prefix+'_SB1_initcont_exec0_shift.ms', prefix+'_LB1_initcont_exec0_shift.ms', prefix+'_LB1_initcont_exec1.ms']: export_MS(msfile) """ Plot deprojected visibility profiles for all data together """ plot_deprojected([prefix+'_SB1_initcont_exec0_shift.vis.npz', prefix+'_LB1_initcont_exec0_shift.vis.npz', prefix+'_LB1_initcont_exec1.vis.npz'], fluxscale=[1.0, 1.0, 1.0], offx=offx, offy=offy, PA=PA, incl=incl, show_err=False) """ Now inspect offsets by comparing against a reference """ estimate_flux_scale(reference=prefix+'_SB1_initcont_exec0_shift.vis.npz', comparison=prefix+'_LB1_initcont_exec0_shift.vis.npz', offx=offx, offy=offy, incl=incl, PA=PA) #There is discrepancy, but this turns out to be purely phase noise estimate_flux_scale(reference=prefix+'_SB1_initcont_exec0_shift.vis.npz', comparison=prefix+'_LB1_initcont_exec1.vis.npz', offx=offx, offy=offy, incl=incl, PA=PA) #More phase noise, but there is a genuine offset. """ Correct the flux scales where appropriate. """ rescale_flux(prefix+'_LB1_initcont_exec1.ms', [0.884]) """ SELF-CAL for short-baseline data """ """ Copy the single SB execution """ SB_cont_p0 = prefix+'_SB_contp0' os.system('rm -rf %s*' % SB_cont_p0) os.system('cp -r '+prefix+'_SB1_initcont_exec0_shift.ms '+SB_cont_p0+'.ms') """ Set up a clean mask """ mask_ra = '15h56m40.005s' mask_dec = '-22.01.40.398' 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) #HD142666_SB_contp0.image #Beam 0.234 arcsec x 0.195 arcsec (60.13 deg) #Flux inside disk mask: 116.30 mJy #Peak intensity of source: 34.43 mJy/beam #rms: 2.19e-01 mJy/beam #Peak SNR: 156.86 """ Self-calibration parameters """ SB_contspws = '0~7' SB_refant = 'DV09, DV08, DV06' SB1_obs0_timerange = '2015/07/20/00~2015/07/22/00' """ First round of phase-only self-cal (short baselines only) """ 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='120s', 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_obs0_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) #HD142666_SB_contp1.image #Beam 0.236 arcsec x 0.199 arcsec (56.00 deg) #Flux inside disk mask: 124.35 mJy #Peak intensity of source: 39.88 mJy/beam #rms: 6.59e-02 mJy/beam #Peak SNR: 604.83 """ 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_obs0_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) #HD142666_SB_contp2.image #Beam 0.239 arcsec x 0.206 arcsec (48.39 deg) #Flux inside disk mask: 124.85 mJy #Peak intensity of source: 41.90 mJy/beam #rms: 6.51e-02 mJy/beam #Peak SNR: 643.11 """ Third round of phase-only self-cal (short baselines only) """ SB_p3 = prefix+'_SB.p3' os.system('rm -rf '+SB_p3) gaincal(vis=SB_cont_p2+'.ms', caltable=SB_p3, gaintype='T', spw=SB_contspws, refant=SB_refant, calmode='p', solint='30s', minsnr=1.5, minblperant=4) if not skip_plots: """ Inspect gain tables """ plotcal(caltable=SB_p3, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=SB1_obs0_timerange, plotrange=[0,0,-180,180]) """ Apply the solutions """ applycal(vis=SB_cont_p2+'.ms', spw=SB_contspws, gaintable=[SB_p3], interp='linearPD', calwt=True) """ Split off a corrected MS """ SB_cont_p3 = prefix+'_SB_contp3' os.system('rm -rf %s*' % SB_cont_p3) split(vis=SB_cont_p2+'.ms', outputvis=SB_cont_p3+'.ms', datacolumn='corrected') """ Image the results; check the resulting map """ tclean_wrapper(vis=SB_cont_p3+'.ms', imagename=SB_cont_p3, mask=common_mask, scales=SB_scales, threshold='0.07mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p3+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD142666_SB_contp3.image #Beam 0.241 arcsec x 0.208 arcsec (44.98 deg) #Flux inside disk mask: 125.29 mJy #Peak intensity of source: 43.11 mJy/beam #rms: 6.37e-02 mJy/beam #Peak SNR: 676.57 """ Fourth round of phase-only self-cal (short baselines only) """ SB_p4 = prefix+'_SB.p4' os.system('rm -rf '+SB_p4) gaincal(vis=SB_cont_p3+'.ms', caltable=SB_p4, gaintype='T', spw=SB_contspws, refant=SB_refant, calmode='p', solint='18s', minsnr=1.5, minblperant=4) if not skip_plots: """ Inspect gain tables """ plotcal(caltable=SB_p4, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=SB1_obs0_timerange, plotrange=[0,0,-180,180]) """ Apply the solutions """ applycal(vis=SB_cont_p3+'.ms', spw=SB_contspws, gaintable=[SB_p4], interp='linearPD', calwt=True) """ Split off a corrected MS """ SB_cont_p4 = prefix+'_SB_contp4' os.system('rm -rf %s*' % SB_cont_p4) split(vis=SB_cont_p3+'.ms', outputvis=SB_cont_p4+'.ms', datacolumn='corrected') """ Image the results; check the resulting map """ tclean_wrapper(vis=SB_cont_p4+'.ms', imagename=SB_cont_p4, mask=common_mask, scales=SB_scales, threshold='0.07mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p4+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD142666_SB_contp4.image #Beam 0.244 arcsec x 0.213 arcsec (37.91 deg) #Flux inside disk mask: 125.77 mJy #Peak intensity of source: 45.01 mJy/beam #rms: 6.40e-02 mJy/beam #Peak SNR: 703.39 """ Amplitude self-cal (short baselines only) """ SB_ap = prefix+'_SB.ap' os.system('rm -rf '+SB_ap) gaincal(vis=SB_cont_p4+'.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_obs0_timerange, plotrange=[0,0,0,2]) """ Apply the solutions """ applycal(vis=SB_cont_p4+'.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_p4+'.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.07mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_ap+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD142666_SB_contap.image #Beam 0.242 arcsec x 0.212 arcsec (37.16 deg) #Flux inside disk mask: 126.14 mJy #Peak intensity of source: 44.63 mJy/beam #rms: 5.01e-02 mJy/beam #Peak SNR: 890.02 """ 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_rescaled.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.1mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p0+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD142666_combined_contp0.image #Beam 0.033 arcsec x 0.023 arcsec (64.61 deg) #Flux inside disk mask: 127.79 mJy #Peak intensity of source: 1.42 mJy/beam #rms: 1.49e-02 mJy/beam #Peak SNR: 95.17 """ Self-calibration parameters """ combined_contspws = '0~15' combined_refant = 'DV24@A090,DA61@A015,DV09@A007,DV06@A027,DV08@A008,DV06@A014' combined_spwmap = [0,0,0,0,0,0,0,0,8,8,8,8,12,12,12,12] LB1_obs0_timerange = '2017/09/24/00~2017/09/26/00' LB2_obs1_timerange = '2017/11/08/00~2017/11/10/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.1mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p1+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD142666_combined_contp1.image #Beam 0.033 arcsec x 0.023 arcsec (64.61 deg) #Flux inside disk mask: 127.77 mJy #Peak intensity of source: 1.33 mJy/beam #rms: 1.55e-02 mJy/beam #Peak SNR: 86.14 """ 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.03mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p2+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD142666_combined_contp2.image #Beam 0.033 arcsec x 0.023 arcsec (64.61 deg) #Flux inside disk mask: 128.53 mJy #Peak intensity of source: 1.38 mJy/beam #rms: 1.46e-02 mJy/beam #Peak SNR: 94.80 """ 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.03mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p3+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD142666_combined_contp3.image #Beam 0.033 arcsec x 0.023 arcsec (64.61 deg) #Flux inside disk mask: 128.35 mJy #Peak intensity of source: 1.41 mJy/beam #rms: 1.43e-02 mJy/beam #Peak SNR: 98.89 """ 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.03mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p4+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD142666_combined_contp4.image #Beam 0.033 arcsec x 0.023 arcsec (64.61 deg) #Flux inside disk mask: 128.46 mJy #Peak intensity of source: 1.46 mJy/beam #rms: 1.42e-02 mJy/beam #Peak SNR: 103.25 """ Fifth round of phase-only self-cal (all data) """ combined_p5 = prefix+'_combined.p5' os.system('rm -rf '+combined_p5) gaincal(vis=combined_cont_p4+'.ms', caltable=combined_p5, gaintype='T', combine='spw,scan', spw=combined_contspws, refant=combined_refant, calmode='p', solint='30s', minsnr=1.5, minblperant=4) if not skip_plots: """ Inspect gain tables """ plotcal(caltable=combined_p5, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=LB1_obs0_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=combined_p5, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=LB2_obs1_timerange, plotrange=[0,0,-180,180]) """ Apply the solutions """ applycal(vis=combined_cont_p4+'.ms', spw=combined_contspws, spwmap=combined_spwmap, gaintable=[combined_p5], interp='linearPD', calwt=True, applymode='calonly') """ Split off a corrected MS """ combined_cont_p5 = prefix+'_combined_contp5' os.system('rm -rf %s*' % combined_cont_p5) split(vis=combined_cont_p4+'.ms', outputvis=combined_cont_p5+'.ms', datacolumn='corrected') """ Image the results; check the resulting map """ tclean_wrapper(vis=combined_cont_p5+'.ms', imagename=combined_cont_p5, mask=common_mask, scales=LB_scales, threshold='0.03mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p5+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD142666_combined_contp5.image #Beam 0.033 arcsec x 0.023 arcsec (64.61 deg) #Flux inside disk mask: 128.52 mJy #Peak intensity of source: 1.49 mJy/beam #rms: 1.42e-02 mJy/beam #Peak SNR: 104.61 """ Amplitude self-cal (all data) """ combined_ap = prefix+'_combined.ap' os.system('rm -rf '+combined_ap) gaincal(vis=combined_cont_p5+'.ms', caltable=combined_ap, gaintype='T', combine='spw,scan', spw=combined_contspws, refant=combined_refant, calmode='ap', solint='900s', minsnr=3.0, minblperant=4, solnorm=False) if not skip_plots: """ Inspect gain tables """ plotcal(caltable=combined_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=LB1_obs0_timerange, plotrange=[0,0,0,2]) plotcal(caltable=combined_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=LB2_obs1_timerange, plotrange=[0,0,0,2]) """ Apply the solutions """ applycal(vis=combined_cont_p5+'.ms', spw=combined_contspws, spwmap=combined_spwmap, gaintable=[combined_ap], interp='linearPD', calwt=True, applymode='calonly') """ Split off a corrected MS """ combined_cont_ap = prefix+'_combined_contap' os.system('rm -rf %s*' % combined_cont_ap) split(vis=combined_cont_p5+'.ms', outputvis=combined_cont_ap+'.ms', datacolumn='corrected') """ Image the results; check the resulting map """ tclean_wrapper(vis=combined_cont_ap+'.ms', imagename=combined_cont_ap, mask=common_mask, scales=LB_scales, threshold='0.02mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_ap+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD142666_combined_contap.image #Beam 0.031 arcsec x 0.022 arcsec (62.64 deg) #Flux inside disk mask: 130.60 mJy #Peak intensity of source: 1.28 mJy/beam #rms: 1.32e-02 mJy/beam #Peak SNR: 97.08 """ Final outputs """ """ Save the final MS """ os.system('cp -r '+combined_cont_ap+'.ms '+prefix+'_continuum.ms') os.system('tar cvzf '+prefix+'_continuum.ms.tgz '+prefix+'_continuum.ms') """ Make a fiducial continuum image (based on experimentation) """ tclean_wrapper(vis=combined_cont_ap+'.ms', imagename=prefix+'_continuum', mask=common_mask, scales=LB_scales, threshold='0.02mJy') exportfits(prefix+'_continuum.image', prefix+'_continuum.fits')