""" 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 29 September 2017 and 21 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 = 'RULup' 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': 'RU_Lup', 'line_spws': np.array([0, 4]), # CO SPWs 'line_freqs': np.array([2.30538e11, 2.30538e11]), }, 'LB1': {'vis' : LB1_path, 'name' : 'LB1', 'field' : 'RU_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) # Flagchannels input string for SB1: '0:1856~1982, 4:1856~1982' # Averaged continuum dataset saved to RULup_SB1_initcont.ms # Flagchannels input string for LB1: '3:1858~1984, 7:1858~1984' # Averaged continuum dataset saved to RULup_LB1_initcont.ms """ Define simple masks and clean scales for imaging """ mask_pa = 0 # position angle of mask in degrees mask_maj = 1.2 # semimajor axis of mask in arcsec mask_min = 1.2 # semiminor axis of mask in arcsec SB1_mask = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ ('15h56m42.29s', '-37.49.15.89', mask_maj, mask_min, mask_pa) LB1_mask = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ ('15h56m42.29s', '-37.49.15.89', mask_maj, mask_min, mask_pa) SB_scales = [0, 5, 10, 20] LB_scales = [0, 5, 20, 60, 120, 240] 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.15mJy', 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 : ICRS 15h56m42.294926s -37d49m15.88160s #PA of Gaussian component: 124.34 deg #Inclination of Gaussian component: 16.91 deg fit_gaussian(prefix+'_SB1_initcont_exec1.image', region=SB1_mask) #Peak : ICRS 15h56m42.294294s -37d49m15.89431s #PA of Gaussian component: 124.94 deg #Inclination of Gaussian component: 20.29 deg fit_gaussian(prefix+'_LB1_initcont_exec0.image', region=LB1_mask) #Peak : ICRS 15h56m42.293691s -37d49m15.88867s #PA of Gaussian component: 115.86 deg #Inclination of Gaussian component: 19.12 deg fit_gaussian(prefix+'_LB1_initcont_exec1.image', region=LB1_mask) #Peak : ICRS 15h56m42.293915s -37d49m15.88919s #PA of Gaussian component: 104.94 deg #Inclination of Gaussian component: 20.46 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 15h56m42.295491s -37d49m15.97202s') SB1_shift = prefix+'_SB1_initcont_exec1_shift.ms' os.system('rm -rf '+SB1_shift+'*') fixvis(vis=prefix+'_SB1_initcont_exec1.ms', outputvis=SB1_shift, field=data_params['SB1']['field'], phasecenter='ICRS 15h56m42.295491s -37d49m15.97202s') 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 15h56m42.295491s -37d49m15.97202s') """ Now that everything is aligned, we inspect the flux calibration. """ if not skip_plots: """ Assign rough emission geometry parameters and centers. """ PA, incl = 115., 20. phasecenter = au.radec2deg('15:56:42.295491, -37.49.15.97202') peakpos = au.radec2deg('15:56:42.294, -37.49.15.889') 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+'_SB1_initcont_exec1_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+'_SB1_initcont_exec1_shift.vis.npz', prefix+'_LB1_initcont_exec0_shift.vis.npz', prefix+'_LB1_initcont_exec1.vis.npz'], fluxscale=[1.0, 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+'_SB1_initcont_exec1_shift.vis.npz', offx=offx, offy=offy, incl=incl, PA=PA) #The ratio of comparison : reference is 1.02286 #The scaling factor for gencal is 1.011 for your comparison measurement 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) #The ratio of comparison : reference is 1.07773 #The scaling factor for gencal is 1.038 for your comparison measurement 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) #The ratio of comparison : reference is 0.83204 #The scaling factor for gencal is 0.912 for your comparison measurement """ The 2nd LB execution offset is a problem; this calibrator is found to be an issue in other DSHARP datasets. """ """ Correct the flux scales where appropriate. """ rescale_flux(prefix+'_LB1_initcont_exec1.ms', [0.900]) """ 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_exec0_shift.ms', prefix+'_SB1_initcont_exec1_shift.ms'], concatvis=SB_cont_p0+'.ms', dirtol='0.1arcsec', copypointing=False) """ Set up a clean mask """ mask_PA = 115 mask_maj = 1.2 mask_min = 1.1 mask_ra = '15h56m42.294s' mask_dec = '-37.49.15.889' 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.15mJy', 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) #RULup_SB_contp0.image #Beam 0.275 arcsec x 0.231 arcsec (-82.52 deg) #Flux inside disk mask: 197.37 mJy #Peak intensity of source: 61.48 mJy/beam #rms: 1.11e-01 mJy/beam #Peak SNR: 555.94 """ Self-calibration parameters """ SB_contspws = '0~7' SB_refant = 'DA49, DA59' 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='30s', 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=common_mask, scales=SB_scales, threshold='0.1mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p1+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #RULup_SB_contp1.image #Beam 0.275 arcsec x 0.231 arcsec (-82.52 deg) #Flux inside disk mask: 198.84 mJy #Peak intensity of source: 64.39 mJy/beam #rms: 4.07e-02 mJy/beam #Peak SNR: 1581.29 """ 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='18s', 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=common_mask, scales=SB_scales, threshold='0.07mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p2+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #RULup_SB_contp2.image #Beam 0.275 arcsec x 0.231 arcsec (-82.53 deg) #Flux inside disk mask: 198.95 mJy #Peak intensity of source: 64.64 mJy/beam #rms: 4.08e-02 mJy/beam #Peak SNR: 1585.75 """ 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_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_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.07mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_ap+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #RULup_SB_contap.image #Beam 0.274 arcsec x 0.230 arcsec (-82.50 deg) #Flux inside disk mask: 198.80 mJy #Peak intensity of source: 64.48 mJy/beam #rms: 3.04e-02 mJy/beam #Peak SNR: 2120.24 """ 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.05mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p0+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #RULup_combined_contp0.image #Beam 0.037 arcsec x 0.030 arcsec (88.65 deg) #Flux inside disk mask: 200.02 mJy #Peak intensity of source: 4.47 mJy/beam #rms: 1.77e-02 mJy/beam #Peak SNR: 252.14 """ Self-calibration parameters """ combined_contspws = '0~15' combined_refant = 'DV24@A090,DV08@A042,DV09@A007,DA52@A075,DA49@A002,DA59@A001' combined_spwmap = [0,0,0,0,4,4,4,4,8,8,8,8,12,12,12,12] LB1_obs0_timerange = '2017/09/28/00~2017/09/30/00' LB2_obs1_timerange = '2017/11/20/00~2017/11/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.04mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p1+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #RULup_combined_contp1.image #Beam 0.037 arcsec x 0.030 arcsec (88.65 deg) #Flux inside disk mask: 200.32 mJy #Peak intensity of source: 4.60 mJy/beam #rms: 1.50e-02 mJy/beam #Peak SNR: 305.88 """ 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) #RULup_combined_contp2.image #Beam 0.037 arcsec x 0.030 arcsec (88.65 deg) #Flux inside disk mask: 199.79 mJy #Peak intensity of source: 4.67 mJy/beam #rms: 1.52e-02 mJy/beam #Peak SNR: 306.51 """ 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) #RULup_combined_contp3.image #Beam 0.037 arcsec x 0.030 arcsec (88.65 deg) #Flux inside disk mask: 199.30 mJy #Peak intensity of source: 4.94 mJy/beam #rms: 1.42e-02 mJy/beam #Peak SNR: 347.53 """ 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) #RULup_combined_contp4.image #Beam 0.037 arcsec x 0.030 arcsec (88.65 deg) #Flux inside disk mask: 199.16 mJy #Peak intensity of source: 5.31 mJy/beam #rms: 1.38e-02 mJy/beam #Peak SNR: 383.98 """ 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) #RULup_combined_contp5.image #Beam 0.037 arcsec x 0.030 arcsec (88.65 deg) #Flux inside disk mask: 198.78 mJy #Peak intensity of source: 5.48 mJy/beam #rms: 1.42e-02 mJy/beam #Peak SNR: 386.89 """ 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) #RULup_combined_contap.image #Beam 0.036 arcsec x 0.029 arcsec (-88.53 deg) #Flux inside disk mask: 202.74 mJy #Peak intensity of source: 5.15 mJy/beam #rms: 1.26e-02 mJy/beam #Peak SNR: 406.75 """ 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) """ mask_PA, mask_maj, mask_min = 0, 1.2, 1.2 fid_mask = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ ('15h56m42.29s', '-37.49.15.89', mask_maj, mask_min, mask_PA) scales = [0, 5, 30, 75] tclean_wrapper(vis=combined_cont_ap+'.ms', imagename=prefix+'_continuum', mask=fid_mask, scales=scales, threshold='0.08mJy', robust=-0.5, uvtaper=['0.022arcsec', '0.01arcsec', '-6deg']) exportfits(prefix+'_continuum.image', prefix+'_continuum.fits')