""" This script was written for CASA 5.1.1 Datasets calibrated (in order of date observed): SB1: 2016.1.00484.L Observed 14 May 2017 and 17 May 2017 (2 execution blocks) LB1: 2016.1.00484.L Observed 24 September 2017 (2 execution blocks) reducer: N. Kurtovic """ """ Starting matter """ import os execfile('reduction_utils.py') skip_plots = True # if True, can run script non-interactively """ Input for loading data """ prefix = 'HTLup' 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': 'HT_Lup', 'line_spws': np.array([0, 4]), # CO SPWs 'line_freqs': np.array([2.30538e11, 2.30538e11]), }, 'LB1': {'vis' : LB1_path, 'name' : 'LB1', 'field' : 'HT_Lupi', '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 """ # Primary (+ Secondary) ra_p = '15h45m12.848s' dec_p = '-34d17m31.04s' SBmaj_p = 0.48 SBmin_p = 0.42 SBpa_p = 90 LBmaj_p = 0.26 LBmin_p = 0.23 LBpa_p = 90 mask_SB_p = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ (ra_p, dec_p, SBmaj_p, SBmin_p, SBpa_p) mask_LB_p = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ (ra_p, dec_p, LBmaj_p, LBmin_p, LBpa_p) # Tertiary ra_t = '15h45m12.645s' dec_t = '-34d17m29.74s' SBmaj_t = 0.45 SBmin_t = 0.30 SBpa_t = 90 LBmaj_t = 0.15 LBmin_t = 0.12 LBpa_t = 90 mask_SB_t = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ (ra_t, dec_t, SBmaj_t, SBmin_t, SBpa_t) mask_LB_t = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ (ra_t, dec_t, LBmaj_t, LBmin_t, LBpa_t) # Combined mask_SB = [mask_SB_p, mask_SB_t] noise_annulus_SB = "annulus[[%s, %s], ['4arcsec', '9arcsec']]" % (ra_p, dec_p) mask_LB = [mask_LB_p, mask_LB_t] noise_annulus_LB = "annulus[[%s, %s], ['4arcsec', '7arcsec']]" % (ra_p, dec_p) SB_scales = [0] LB_scales = [0, 10, 30] if not skip_plots: """ Image each dataset individually """ # images are saved in the format prefix+'_name_initcont_exec#.ms' tclean_wrapper(vis=prefix+'_SB_initcont.ms', imagename=prefix+'_SB_initcont', scales=SB_scales, mask=mask_SB, savemodel='modelcolumn', threshold='1.0mJy', interactive=False) image_each_obs(data_params['SB1'], prefix, mask=mask_SB, scales=SB_scales) tclean_wrapper(vis=prefix+'_LB_initcont.ms', imagename=prefix+'_LB_initcont', scales=LB_scales, mask=mask_SB, savemodel='modelcolumn', threshold='0.1mJy', interactive=False) image_each_obs(data_params['LB1'], prefix, mask=mask_LB, scales=LB_scales) """ Fit Gaussians to roughly estimate centers, inclinations, PAs """ fit_gaussian(prefix+'_SB1_initcont_exec0.image', region=mask_SB_p) #Peak : ICRS 15h45m12.847769s -34d17m31.03291s fit_gaussian(prefix+'_SB1_initcont_exec1.image', region=mask_SB_p) #Peak : ICRS 15h45m12.847681s -34d17m31.04498s fit_gaussian(prefix+'_SB1_initcont.image', region=mask_SB_p) #Peak : ICRS 15h45m12.847715s -34d17m31.04075s fit_gaussian(prefix+'_SB1_initcont.image', region=mask_SB_t) #Peak : ICRS 15h45m12.645265s -34d17m29.74878s # Peaks are essentially the same in individual executions, so we adopt the # centers from their combination # The LB data resolves the primary and secondary, so we need refined fit # regions to estimate centroids accurately. reg_maj_p = 0.18 reg_min_p = 0.12 reg_pa_p = 162 reg_p = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ (ra_p, dec_p, reg_maj_p, reg_min_p, reg_pa_p) reg_s = 'circle[[%s, %s], %.1farcsec]' % \ ('15h45m12.835s', '-34d17m31.09s', 0.04) fit_gaussian(prefix+'_LB1_initcont_exec0.image', region=reg_p) #Peak : ICRS 15h45m12.846759s -34d17m31.03121s fit_gaussian(prefix+'_LB1_initcont_exec1.image', region=reg_p) #Peak : ICRS 15h45m12.846793s -34d17m31.02979s fit_gaussian(prefix+'_LB1_initcont.image', region=reg_p) #Peak : ICRS 15h45m12.846796s -34d17m31.03108s fit_gaussian(prefix+'_LB1_initcont.image', region=reg_s) #Peak : ICRS 15h45m12.835149s -34d17m31.09333s fit_gaussian(prefix+'_LB1_initcont.image', region=mask_LB_t) #Peak : ICRS 15h45m12.644503s -34d17m29.73807s """ Because we cannot make a direct comparison of the SB and LB centroids for the primary disk (due to the secondary disk blend in the SB data), we use the tertiary centroid measurements to align the datasets. (Not in script: check that these shifts do what you think they should by re-imaging! They do for us.) """ """ Define a common direction (peak of 2nd LB EB) """ common_dir = 'J2000 15h45m12.84731s -34d17m31.016187s' """ Shift each MS so emission center is at same phase center """ SB_shift = prefix+'_SB1_initcont_shift' os.system('rm -rf '+SB_shift+'.ms*') fixvis(vis=prefix+'_SB1_initcont.ms', outputvis=SB_shift+'.ms', field=data_params['SB1']['field'], phasecenter='J2000 15h45m12.64578s -34d17m29.733887s') fixplanets(vis=SB_shift+'.ms', field=data_params['SB1']['field'], direction='J2000 15h45m12.64502s -34d17m29.723177s') fixvis(vis=SB_shift+'.ms', outputvis=SB_shift+'.ms', field=data_params['SB1']['field'], phasecenter=common_dir) LB_shift = prefix+'_LB1_initcont_shift' os.system('rm -rf '+LB_shift+'.ms*') fixvis(vis=prefix+'_LB1_initcont.ms', outputvis=LB_shift+'.ms', field=data_params['LB1']['field'], phasecenter=common_dir) """ Now that everything is aligned, we inspect the flux calibration. """ # Split up the shifted datasets split_all_obs(SB_shift+'.ms', SB_shift+'_exec') split_all_obs(LB_shift+'.ms', LB_shift+'_exec') if not skip_plots: """ Assign rough emission geometry parameters. """ PA, incl = 166, 47 """ Export MS contents into Numpy save files """ for msfile in [prefix+'_SB1_initcont_shift_exec0.ms', prefix+'_SB1_initcont_shift_exec1.ms', prefix+'_LB1_initcont_shift_exec0.ms', prefix+'_LB1_initcont_shift_exec1.ms']: export_MS(msfile) """ Plot deprojected visibility profiles for all data together """ plot_deprojected([prefix+'_SB1_initcont_shift_exec0.vis.npz', prefix+'_SB1_initcont_shift_exec1.vis.npz', prefix+'_LB1_initcont_shift_exec0.vis.npz', prefix+'_LB1_initcont_shift_exec1.vis.npz'], fluxscale=[1.0, 1.0, 1.0, 1.0], PA=PA, incl=incl, show_err=False) """ Correct the flux scales where appropriate. """ os.system('rm -rf *rescaled.*') rescale_flux(prefix+'_LB1_initcont_shift_exec0.ms', [0.975]) rescale_flux(prefix+'_LB1_initcont_shift_exec1.ms', [0.932]) """ SELF-CAL for short-baseline data """ """ Copy the shifted SB measurement set """ SB_cont_p0 = prefix+'_SB_contp0' os.system('rm -rf %s*' % SB_cont_p0) split(vis=SB_shift+'.ms', outputvis=SB_cont_p0+'.ms', datacolumn='data') """ Initial clean """ tclean_wrapper(vis=SB_cont_p0+'.ms', imagename=SB_cont_p0, mask=mask_SB, scales=SB_scales, threshold='0.2mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p0+'.image', disk_mask=mask_SB_p, noise_mask=noise_annulus_SB) #HTLup_SB_contp0.image #Beam 0.273 arcsec x 0.229 arcsec (-83.71 deg) #Flux inside disk mask: 73.14 mJy #Peak intensity of source: 50.57 mJy/beam #rms: 1.29e-01 mJy/beam #Peak SNR: 391.31 """ Self-calibration parameters """ SB_contspws = '0~7' SB_refant = 'DA61' SB1_obs0_timerange = '2017/05/13/00~2017/05/15/00' SB1_obs1_timerange = '2017/05/16/00~2017/05/18/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='60s', 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]) plotcal(caltable=SB_p1, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB1_obs1_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=mask_SB, scales=SB_scales, threshold='0.2mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p1+'.image', disk_mask=mask_SB_p, noise_mask=noise_annulus_SB) #HTLup_SB_contp1.image #Beam 0.273 arcsec x 0.229 arcsec (-83.72 deg) #Flux inside disk mask: 76.04 mJy #Peak intensity of source: 56.12 mJy/beam #rms: 3.86e-02 mJy/beam #Peak SNR: 1453.15 """ 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='30s', 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]) plotcal(caltable=SB_p2, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=SB1_obs1_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=mask_SB, scales=SB_scales, threshold='0.076mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p2+'.image', disk_mask=mask_SB_p, noise_mask=noise_annulus_SB) #HTLup_SB_contp2.image #Beam 0.273 arcsec x 0.229 arcsec (-83.72 deg) #Flux inside disk mask: 76.14 mJy #Peak intensity of source: 56.91 mJy/beam #rms: 3.72e-02 mJy/beam #Peak SNR: 1529.13 """ 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='18s', 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]) plotcal(caltable=SB_p3, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=SB1_obs1_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=mask_SB, scales=SB_scales, threshold='0.075mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p3+'.image', disk_mask=mask_SB_p, noise_mask=noise_annulus_SB) #HTLup_SB_contp3.image #Beam 0.273 arcsec x 0.229 arcsec (-83.72 deg) #Flux inside disk mask: 76.24 mJy #Peak intensity of source: 57.24 mJy/beam #rms: 3.74e-02 mJy/beam #Peak SNR: 1531.67 """ Amplitude self-cal (short baselines only) """ SB_ap = prefix+'_SB.ap' os.system('rm -rf '+SB_ap) gaincal(vis=SB_cont_p3+'.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]) plotcal(caltable=SB_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=SB1_obs1_timerange, plotrange=[0,0,0,2]) """ Apply the solutions """ applycal(vis=SB_cont_p3+'.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_p3+'.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=mask_SB, scales=SB_scales, threshold='0.055mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_ap+'.image', disk_mask=mask_SB_p, noise_mask=noise_annulus_SB) #HTLup_SB_contap.image #Beam 0.273 arcsec x 0.228 arcsec (-83.56 deg) #Flux inside disk mask: 76.25 mJy #Peak intensity of source: 57.17 mJy/beam #rms: 2.94e-02 mJy/beam #Peak SNR: 1944.31 """ 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_shift_exec0_rescaled.ms', prefix+'_LB1_initcont_shift_exec1_rescaled.ms'], concatvis=combined_cont_p0+'.ms', dirtol='0.1arcsec', copypointing=False) """ Initial clean """ LB_scales = [0, 11, 33, 55] tclean_wrapper(vis=combined_cont_p0+'.ms', imagename=combined_cont_p0, mask=mask_LB, scales=LB_scales, threshold='0.1mJy', imsize=4000, savemodel='modelcolumn') estimate_SNR(combined_cont_p0+'.image', disk_mask=mask_LB_p, noise_mask=noise_annulus_LB) #HTLup_combined_p0.image #Beam 0.038 arcsec x 0.033 arcsec (69.45 deg) #Flux inside disk mask: 77.26 mJy #Peak intensity of source: 7.01 mJy/beam #rms: 1.69e-02 mJy/beam #Peak SNR: 415.83 """ Self-calibration parameters """ combined_contspws = '0~15' combined_refant = 'DA61' combined_spwmap = [0,0,0,0,4,4,4,4,8,8,8,8,12,12,12,12] LB1_obs0_timerange = '2017/09/23/00~2017/09/26/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='360s', 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]) """ 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=mask_LB, scales=LB_scales, threshold='0.034mJy', imsize=4000, savemodel='modelcolumn') estimate_SNR(combined_cont_p1+'.image', disk_mask=mask_LB_p, noise_mask=noise_annulus_LB) #HTLup_combined_p1.image #Beam 0.038 arcsec x 0.033 arcsec (69.45 deg) #Flux inside disk mask: 76.59 mJy #Peak intensity of source: 7.24 mJy/beam #rms: 1.60e-02 mJy/beam #Peak SNR: 454.10 """ 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='180s', 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]) """ 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=mask_LB, scales=LB_scales, threshold='0.032mJy', imsize=4000, savemodel='modelcolumn') estimate_SNR(combined_cont_p2+'.image', disk_mask=mask_LB_p, noise_mask=noise_annulus_LB) #HTLup_combined_p2.image #Beam 0.038 arcsec x 0.033 arcsec (69.45 deg) #Flux inside disk mask: 76.62 mJy #Peak intensity of source: 7.49 mJy/beam #rms: 1.52e-02 mJy/beam #Peak SNR: 493.23 """ 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='60s', 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]) """ 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=mask_LB, scales=LB_scales, threshold='0.030mJy', imsize=4000, savemodel='modelcolumn') estimate_SNR(combined_cont_p3+'.image', disk_mask=mask_LB_p, noise_mask=noise_annulus_LB) #HTLup_combined_p3.image #Beam 0.038 arcsec x 0.033 arcsec (69.45 deg) #Flux inside disk mask: 76.31 mJy #Peak intensity of source: 8.00 mJy/beam #rms: 1.41e-02 mJy/beam #Peak SNR: 568.22 """ 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='30s', 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]) """ 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=mask_LB, scales=LB_scales, threshold='0.022mJy', imsize=4000, savemodel='modelcolumn') estimate_SNR(combined_cont_p4+'.image', disk_mask=mask_LB_p, noise_mask=noise_annulus_LB) #HTLup_combined_contp4.image #Beam 0.038 arcsec x 0.033 arcsec (69.45 deg) #Flux inside disk mask: 76.44 mJy #Peak intensity of source: 8.27 mJy/beam #rms: 1.41e-02 mJy/beam #Peak SNR: 587.77 """ Amplitude self-cal (all data) """ combined_ap = prefix+'_combined.ap' os.system('rm -rf '+combined_ap) gaincal(vis=combined_cont_p4+'.ms', caltable=combined_ap, gaintype='T', combine='spw,scan', spw=combined_contspws, refant=combined_refant, calmode='ap', solint='360s', 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]) """ Apply the solutions """ applycal(vis=combined_cont_p4+'.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_p4+'.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=mask_LB, scales=LB_scales, threshold='0.028mJy', imsize=4000, savemodel='modelcolumn') estimate_SNR(combined_cont_ap+'.image', disk_mask=mask_LB_p, noise_mask=noise_annulus_LB) #HTLup_combined_contap.image #Beam 0.038 arcsec x 0.033 arcsec (61.46 deg) #Flux inside disk mask: 76.71 mJy #Peak intensity of source: 8.28 mJy/beam #rms: 1.41e-02 mJy/beam #Peak SNR: 586.95 """ 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=mask_LB, scales=LB_scales, threshold='0.028mJy', imsize=4000, savemodel='modelcolumn') exportfits(prefix+'_continuum.image', prefix+'_continuum.fits')