""" This script was written for CASA 5.1.1 Datasets calibrated (in order of date observed): SB1: 2015.1.00964.S Observed 14 June 2016 (2 execution blocks) SB2: 2013.1.00964.S Observed 02 July 2016 (1 execution block) SB3: 2015.1.00484.L Observed 14 May 2017, 17 May 2017, 19 May 2017 (3 execution blocks) LB1: 2016.1.00484.L Observed 26 September 2017 and 26 November 2017 (2 execution blocks) reducer: L. Perez """ """ Starting matter """ import os execfile('reduction_utils.py') skip_plots = True # if True, can run script non-interactively """ Input for loading data """ prefix = 'HD143006' SB1_path = '/full_path/to_calibrated/msfile.ms' SB2_path = '/full_path/to_calibrated/msfile.ms' SB3_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': 'HD143006', 'line_spws': np.array([0,9]), # CO SPWs 'line_freqs': np.array([2.30538e11,2.30538e11]), 'cont_spws': '0,1,2,9,10,11', 'width_array': [960,960,256,960,960,256], }, 'SB2': {'vis' : SB2_path, 'name' : 'SB2', 'field': 'HD143006', 'line_spws': np.array([0]), # CO SPWs 'line_freqs': np.array([2.30538e11]), 'cont_spws': '0,1,2', 'width_array': [960,960,256], }, 'SB3': {'vis' : SB3_path, 'name' : 'SB3', 'field': 'HD_143006', 'line_spws': np.array([0,4,8]), # CO SPWs 'line_freqs': np.array([2.30538e11,2.30538e11,2.30538e11]), 'cont_spws': None, 'width_array': None, }, 'LB1': {'vis' : LB1_path, 'name' : 'LB1', 'field' : 'HD_143006', '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') """ Additional flagging """ flagmanager(vis=data_params['SB2']['vis'], mode='save', versionname='original_flags', comment='Original flag states') flagdata(vis=data_params['SB2']['vis'], mode='manual', spw='0:0~200', flagbackup=False, field=data_params['SB2']['field']) """ 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, contspws=data_params[i]['cont_spws'], width_array=data_params[i]['width_array']) """ Define simple masks and clean scales for imaging """ mask_rad = 0.8 # radius of mask in arcsec mask_ra = '15h58m36.9s' mask_dec = '-22.57.15.60' SB1_mask = 'circle[[%s, %s], %.1farcsec]' % (mask_ra, mask_dec, mask_rad) LB1_mask = 'circle[[%s, %s], %.1farcsec]' % (mask_ra, mask_dec, mask_rad) SB_scales = [0, 10, 30] LB_scales = [0, 25, 50, 75, 100] 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.8mJy', robust=-2, interactive=False) image_each_obs(data_params['SB2'], prefix, mask=SB1_mask, scales=SB_scales, threshold='0.5mJy', robust=-2, interactive=False) image_each_obs(data_params['SB3'], prefix, mask=SB1_mask, scales=SB_scales, threshold='1.0mJy', robust=-2, interactive=False) image_each_obs(data_params['LB1'], prefix, mask=LB1_mask, scales=LB_scales, threshold='0.06mJy', interactive=False) """ The emission appears to be aligned. But, we want to force everything to the same phase center to avoid confusion in the visibilities. """ for i in data_params.keys(): shift = prefix+'_'+i+'_initcont_shift.ms' os.system('rm -rf '+shift+'*') split(vis=prefix+'_'+i+'_initcont.ms', outputvis=shift, datacolumn='data') fixvis(vis=shift, outputvis=shift, field=data_params[i]['field'], phasecenter='J2000 15h58m36.89967s -22d57m15.602730s') """ Now that everything is aligned, we inspect the flux calibration. """ split_all_obs(prefix+'_SB1_initcont_shift.ms', prefix+'_SB1_initcont_exec') split_all_obs(prefix+'_SB2_initcont_shift.ms', prefix+'_SB2_initcont_exec') split_all_obs(prefix+'_SB3_initcont_shift.ms', prefix+'_SB3_initcont_exec') split_all_obs(prefix+'_LB1_initcont.ms', prefix+'_LB1_initcont_exec') if not skip_plots: """ Assign rough emission geometry parameters and centers. """ PA, incl = 80, 10 phasecenter = au.radec2deg('15:58:36.899002, -22.57.15.61701') peakpos = au.radec2deg('15:58:36.899, -22.57.15.596') offsets = au.angularSeparation(peakpos[0], peakpos[1], phasecenter[0], phasecenter[1], True) offx, offy = 3600.*offsets[3], 3600.*offsets[2] # since the SB3 and LB1 data only overlap in frequency in the LSB with the # SB1, SB2 data, we split out to make the comparison more fair split(vis=prefix+'_SB3_initcont_exec0.ms', outputvis=prefix+'_SB3_initcont_exec0_LSB.ms', spw='0,1', datacolumn='data') split(vis=prefix+'_SB3_initcont_exec1.ms', outputvis=prefix+'_SB3_initcont_exec1_LSB.ms', spw='0,1', datacolumn='data') split(vis=prefix+'_SB3_initcont_exec2.ms', outputvis=prefix+'_SB3_initcont_exec2_LSB.ms', spw='0,1', datacolumn='data') """ Export MS contents into Numpy save files """ for msfile in [prefix+'_SB1_initcont_exec0.ms', prefix+'_SB1_initcont_exec1.ms', prefix+'_SB2_initcont_exec0.ms', prefix+'_SB3_initcont_exec0.ms', prefix+'_SB3_initcont_exec0_LSB.ms', prefix+'_SB3_initcont_exec1_LSB.ms', prefix+'_SB3_initcont_exec2_LSB.ms', prefix+'_LB1_initcont_exec0.ms', prefix+'_LB1_initcont_exec1.ms']: export_MS(msfile) """ Plot deprojected visibility profiles for all data together """ plot_deprojected([prefix+'_SB1_initcont_exec0.vis.npz', prefix+'_SB1_initcont_exec1.vis.npz', prefix+'_SB2_initcont_exec0.vis.npz', prefix+'_SB3_initcont_exec0_LSB.vis.npz', prefix+'_SB3_initcont_exec1_LSB.vis.npz', prefix+'_SB3_initcont_exec2_LSB.vis.npz', prefix+'_LB1_initcont_exec0.vis.npz', prefix+'_LB1_initcont_exec1.vis.npz'], fluxscale=[1.,1.,1.,1.,1.,1.,1.,1.], PA=PA, incl=incl, offx=offx, offy=offy, show_err=False) """ Now inspect offsets by comparing against a reference """ estimate_flux_scale(reference=prefix+'_LB1_initcont_exec0.vis.npz', comparison=prefix+'_LB1_initcont_exec1.vis.npz', incl=incl, PA=PA, offx=offx, offy=offy) #The ratio of comparison : reference is 0.76910 #The scaling factor for gencal is 0.877 for your comparison measurement estimate_flux_scale(reference=prefix+'_SB3_initcont_exec0_LSB.vis.npz', comparison=prefix+'_SB1_initcont_exec0.vis.npz', incl=incl, PA=PA, offx=offx, offy=offy) #The ratio of comparison : reference is 1.12843 #The scaling factor for gencal is 1.062 for your comparison measurement estimate_flux_scale(reference=prefix+'_SB3_initcont_exec0_LSB.vis.npz', comparison=prefix+'_SB1_initcont_exec1.vis.npz', incl=incl, PA=PA, offx=offx, offy=offy) #The ratio of comparison : reference is 1.04710 #The scaling factor for gencal is 1.023 for your comparison measurement estimate_flux_scale(reference=prefix+'_SB3_initcont_exec0_LSB.vis.npz', comparison=prefix+'_SB2_initcont_exec0.vis.npz', incl=incl, PA=PA, offx=offx, offy=offy) #The ratio of comparison : reference is 1.00313 #The scaling factor for gencal is 1.002 for your comparison measurement estimate_flux_scale(reference=prefix+'_SB3_initcont_exec0_LSB.vis.npz', comparison=prefix+'_SB3_initcont_exec1_LSB.vis.npz', incl=incl, PA=PA, offx=offx, offy=offy) #The ratio of comparison : reference is 1.07657 #The scaling factor for gencal is 1.038 for your comparison measurement estimate_flux_scale(reference=prefix+'_SB3_initcont_exec0.vis.npz', comparison=prefix+'_SB3_initcont_exec2.vis.npz', incl=incl, PA=PA, offx=offx, offy=offy) #The ratio of comparison : reference is 1.15169 #The scaling factor for gencal is 1.073 for your comparison measurement estimate_flux_scale(reference=prefix+'_SB3_initcont_exec0.vis.npz', comparison=prefix+'_LB1_initcont_exec0.vis.npz', incl=incl, PA=PA, offx=offx, offy=offy) #The ratio of comparison : reference is 1.16634 #The scaling factor for gencal is 1.080 for your comparison measurement estimate_flux_scale(reference=prefix+'_SB3_initcont_exec0.vis.npz', comparison=prefix+'_LB1_initcont_exec1.vis.npz', incl=incl, PA=PA, offx=offx, offy=offy) #The ratio of comparison : reference is 0.92359 #The scaling factor for gencal is 0.961 for your comparison measurement """ Correct the flux scales where appropriate. """ rescale_flux(prefix+'_SB1_initcont_exec0.ms', [1.062]) rescale_flux(prefix+'_SB1_initcont_exec1.ms', [1.023]) rescale_flux(prefix+'_SB2_initcont_exec0.ms', [1.002]) rescale_flux(prefix+'_SB3_initcont_exec1.ms', [1.038]) rescale_flux(prefix+'_SB3_initcont_exec2.ms', [1.073]) rescale_flux(prefix+'_LB1_initcont_exec0.ms', [1.080]) rescale_flux(prefix+'_LB1_initcont_exec1.ms', [0.961]) """ 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_rescaled.ms', prefix+'_SB1_initcont_exec1_rescaled.ms', prefix+'_SB2_initcont_exec0_rescaled.ms', prefix+'_SB3_initcont_exec0.ms', prefix+'_SB3_initcont_exec1_rescaled.ms', prefix+'_SB3_initcont_exec2_rescaled.ms'], concatvis=SB_cont_p0+'.ms', dirtol='0.1arcsec', copypointing=False) """ Set up a clean mask """ mask_ra = '15h58m36.9s' mask_dec = '-22.57.15.60' common_mask = 'circle[[%s, %s], %.1farcsec]' % (mask_ra, mask_dec, mask_rad) """ Initial clean """ tclean_wrapper(vis=SB_cont_p0+'.ms', imagename=SB_cont_p0, mask=common_mask, scales=SB_scales, threshold='0.04mJy', 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_rad) estimate_SNR(SB_cont_p0+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD143006_SB_contp0.image #Beam 0.287 arcsec x 0.244 arcsec (-89.21 deg) #Flux inside disk mask: 56.02 mJy #Peak intensity of source: 7.19 mJy/beam #rms: 3.98e-02 mJy/beam #Peak SNR: 180.66 """ Self-calibration parameters """ SB_contspws = '0~20' SB_refant = 'DA41@A004,DA49@A002,DV16@A036,DA46@A034,DA51@A023' SB1_obs0_timerange = '2016/06/14/00~2016/06/15/00' SB1_obs2_timerange = '2016/07/02/00~2016/07/03/00' SB1_obs3_timerange = '2017/05/14/00~2017/05/15/00' SB1_obs4_timerange = '2017/05/17/00~2017/05/18/00' SB1_obs5_timerange = '2017/05/18/00~2017/05/20/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_obs1_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p1, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB1_obs2_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p1, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB1_obs3_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p1, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB1_obs4_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p1, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB1_obs5_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.035mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p1+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD143006_SB_contp1.image #Beam 0.288 arcsec x 0.244 arcsec (-89.16 deg) #Flux inside disk mask: 57.74 mJy #Peak intensity of source: 7.94 mJy/beam #rms: 2.34e-02 mJy/beam #Peak SNR: 339.68 """ Second round of phase-only self-cal (short baselines only) """ # will combine SPWs for increased SNR 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, combine='spw', 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_obs1_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p2, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB1_obs2_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p2, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB1_obs3_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p2, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB1_obs4_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p2, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB1_obs5_timerange, plotrange=[0,0,-180,180]) """ Apply the solutions """ SB_spwmap=[0,0,0,3,3,3,6,6,6,9,9,9,9,13,13,13,13,17,17,17,17] applycal(vis=SB_cont_p1+'.ms', spw=SB_contspws, gaintable=[SB_p2], spwmap=SB_spwmap, 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.035mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p2+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD143006_SB_contp2.image #Beam 0.288 arcsec x 0.244 arcsec (-89.93 deg) #Flux inside disk mask: 57.84 mJy #Peak intensity of source: 8.10 mJy/beam #rms: 2.33e-02 mJy/beam #Peak SNR: 347.61 """ 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, combine='scan', 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_obs1_timerange, plotrange=[0,0,0,2]) plotcal(caltable=SB_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=SB1_obs2_timerange, plotrange=[0,0,0,2]) plotcal(caltable=SB_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=SB1_obs3_timerange, plotrange=[0,0,0,2]) plotcal(caltable=SB_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=SB1_obs4_timerange, plotrange=[0,0,0,2]) plotcal(caltable=SB_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=SB1_obs5_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.023mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_ap+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD143006_SB_contap.image #Beam 0.284 arcsec x 0.243 arcsec (-89.87 deg) #Flux inside disk mask: 58.00 mJy #Peak intensity of source: 7.98 mJy/beam #rms: 2.24e-02 mJy/beam #Peak SNR: 355.71 # Additional rounds or variants of amplitude self-cal do not help. """ 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_rescaled.ms', prefix+'_LB1_initcont_exec1_rescaled.ms'], concatvis=combined_cont_p0+'.ms', dirtol='0.1arcsec', copypointing=False) """ Modify mask and scales """ mask_ra = '15h58m36.899s' mask_dec = '-22.57.15.583' mask_rad = 0.7 common_mask = 'circle[[%s, %s], %.1farcsec]' % (mask_ra, mask_dec, mask_rad) noise_annulus = "annulus[[%s, %s],['%.2farcsec', '4.25arcsec']]" % \ (mask_ra, mask_dec, 1.1*mask_rad) LB_scales = [0,17,52,87,170] """ Initial clean """ tclean_wrapper(vis=combined_cont_p0+'.ms', imagename=combined_cont_p0, mask=common_mask, scales=LB_scales, threshold='0.02mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p0+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD143006_combined_contp0.image #Beam 0.052 arcsec x 0.037 arcsec (86.45 deg) #Flux inside disk mask: 57.85 mJy #Peak intensity of source: 0.60 mJy/beam #rms: 1.29e-02 mJy/beam #Peak SNR: 46.44 """ Self-calibration parameters """ combined_contspws = '0~28' combined_refant = 'DV16@A036,DA41@A004,DA49@A002,DA46@A034,DA51@A023,DA61@A015,DA47@A074,DV20@A072,DV08@A042' combined_spwmap = [0,0,0,3,3,3,6,6,6,9,9,9,9,13,13,13,13,17,17,17,17,21,21,21,21,25,25,25,25] LB1_obs0_timerange = '2017/09/26/00~2017/09/27/00' LB2_obs1_timerange = '2017/11/26/00~2017/11/27/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.019mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p1+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD143006_combined_contp1.image #Beam 0.052 arcsec x 0.037 arcsec (86.45 deg) #Flux inside disk mask: 57.95 mJy #Peak intensity of source: 0.61 mJy/beam #rms: 1.27e-02 mJy/beam #Peak SNR: 48.02 # ** map quality improving ** """ 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.019mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p2+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD143006_combined_contp2.image #Beam 0.052 arcsec x 0.037 arcsec (86.45 deg) #Flux inside disk mask: 57.84 mJy #Peak intensity of source: 0.62 mJy/beam #rms: 1.27e-02 mJy/beam #Peak SNR: 49.16 """ 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.019mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p3+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD143006_combined_contp3.image #Beam 0.052 arcsec x 0.037 arcsec (86.45 deg) #Flux inside disk mask: 57.87 mJy #Peak intensity of source: 0.65 mJy/beam #rms: 1.26e-02 mJy/beam #Peak SNR: 51.40 """ 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.019mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p4+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD143006_combined_contp4.image #Beam 0.052 arcsec x 0.037 arcsec (86.45 deg) #Flux inside disk mask: 58.04 mJy #Peak intensity of source: 0.69 mJy/beam #rms: 1.25e-02 mJy/beam #Peak SNR: 55.10 """ 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.013mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p5+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD143006_combined_contp5.image #Beam 0.052 arcsec x 0.037 arcsec (86.45 deg) #Flux inside disk mask: 58.16 mJy #Peak intensity of source: 0.71 mJy/beam #rms: 1.25e-02 mJy/beam #Peak SNR: 56.86 """ 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.012mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_ap+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #HD143006_combined_contap.image #Beam 0.046 arcsec x 0.031 arcsec (78.50 deg) #Flux inside disk mask: 58.25 mJy #Peak intensity of source: 0.56 mJy/beam #rms: 1.14e-02 mJy/beam #Peak SNR: 48.91 """ 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) """ scales = [0, 5, 30, 75] tclean_wrapper(vis=combined_cont_ap+'.ms', imagename=prefix+'_continuum', mask=common_mask, scales=scales, threshold='0.05mJy', robust=0., uvtaper=['.042arcsec', '.020arcsec', '172.1deg']) exportfits(prefix+'_continuum.image', prefix+'_continuum.fits')