""" 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, 17 May 2017, and 19 May 2017 (3 execution blocks) LB1: 2016.1.00484.L Observed 29 September 2017 (1 execution block) 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 = 'AS205' 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' : 'V866 Sco', 'line_spws': np.array([0,2,3]), # CO SPWs 'line_freqs': np.array([2.30538e11, 2.2039868420e11, 2.1956035410e11]), }, 'SB2': {'vis' : SB2_path, 'name' : 'SB2', 'field' : 'V866 Sco', 'line_spws': np.array([0,2,3]), # CO SPWs 'line_freqs': np.array([2.30538e11, 2.2039868420e11, 2.1956035410e11]), }, 'SB3': {'vis' : SB3_path, 'name' : 'SB3', 'field': 'AS_205', 'line_spws': np.array([0,4,8]), # CO SPWs 'line_freqs': np.array([2.30538e11, 2.30538e11, 2.30538e11]), }, 'LB1': {'vis' : LB1_path, 'name' : 'LB1', 'field' : 'AS_205', 'line_spws': np.array([3]), # CO SPWs 'line_freqs': np.array([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') """ Since SB1 and SB2 are from Cycle 0, need to re-calculate their visibility weights properly before continuum averaging. """ os.system('rm -rf '+prefix+'_SB1_in.ms*') split(vis=data_params['SB1']['vis'], field=data_params['SB1']['field'], spw='0~3', outputvis=prefix+'_SB1_in.ms', datacolumn='data') cchans = '0:0~2200;3500~3839, 1:0~3839, 2:0~650;1860~3839, 3:0~140;1350~3839' fchans = '0:2232~3492, 2:651~1855, 3:148~1348' statwt(vis=prefix+'_SB1_in.ms', fitspw=cchans, datacolumn='data') flagmanager(vis=prefix+'_SB1_in.ms', mode='save', versionname='before_cont_flags', comment='Flag states before spectral lines are flagged') flagdata(vis=prefix+'_SB1_in.ms', mode='manual', spw=fchans, flagbackup=False, field=data_params['SB1']['field']) os.system('rm -rf '+prefix+'_SB1_initcont.ms*') tb.open(prefix+'_SB1_in.ms/SPECTRAL_WINDOW') total_bw = tb.getcol('TOTAL_BANDWIDTH') num_chan = tb.getcol('NUM_CHAN') tb.close() width_array = (num_chan/np.ceil(total_bw/(125e6))).astype('int').tolist() split(vis=prefix+'_SB1_in.ms', field=data_params['SB1']['field'], spw='0~3', outputvis=prefix+'_SB1_initcont.ms', width=width_array, timebin='0s', datacolumn='data', intent='OBSERVE_TARGET#ON_SOURCE', keepflags=False) flagmanager(vis=prefix+'_SB1_in.ms', mode='restore', versionname='before_cont_flags') os.system('rm -rf '+prefix+'_SB2_in.ms*') split(vis=data_params['SB2']['vis'], field=data_params['SB2']['field'], spw='0~3', outputvis=prefix+'_SB2_in.ms', datacolumn='data') cchans = '0:0~1860;3130~3839, 1:0~3839, 2:0~990;2210~3839, 3:0~490;1700~3839' fchans = '0:1868~3128, 2:999~2204, 3:495~1695' statwt(vis=prefix+'_SB2_in.ms', fitspw=cchans, datacolumn='data') flagmanager(vis=prefix+'_SB2_in.ms', mode='save', versionname='before_cont_flags', comment='Flag states before spectral lines are flagged') flagdata(vis=prefix+'_SB2_in.ms', mode='manual', spw=fchans, flagbackup=False, field=data_params['SB2']['field']) os.system('rm -rf '+prefix+'_SB2_initcont.ms*') tb.open(prefix+'_SB2_in.ms/SPECTRAL_WINDOW') total_bw = tb.getcol('TOTAL_BANDWIDTH') num_chan = tb.getcol('NUM_CHAN') tb.close() width_array = (num_chan/np.ceil(total_bw/(125e6))).astype('int').tolist() split(vis=prefix+'_SB2_in.ms', field=data_params['SB2']['field'], spw='0~3', outputvis=prefix+'_SB2_initcont.ms', width=width_array, timebin='0s', datacolumn='data', intent='OBSERVE_TARGET#ON_SOURCE', keepflags=False) flagmanager(vis=prefix+'_SB2_in.ms', mode='restore', versionname='before_cont_flags') """ Identify region containing CO emission; then flag that and do a spectral average to a pseudo-continuum MS """ for i in ['SB3', 'LB1']: 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 ra_p = '16h11m31.351s' dec_p = '-18d38m26.247s' SBmaj_p = 0.82 SBmin_p = 0.55 SBpa_p = 90 LBmaj_p = 0.60 LBmin_p = 0.50 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) # Secondary ra_s = '16h11m31.295s' dec_s = '-18d38m27.276s' SBmaj_s = 0.52 SBmin_s = 0.30 SBpa_s = 108.5 LBmaj_s = 0.48 LBmin_s = 0.26 LBpa_s = 108.5 mask_SB_s = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ (ra_s, dec_s, SBmaj_s, SBmin_s, SBpa_s) mask_LB_s = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ (ra_s, dec_s, LBmaj_s, LBmin_s, LBpa_s) # Combined mask_SB = [mask_SB_p, mask_SB_s] mask_LB = [mask_LB_p, mask_LB_s] noise_annulus_SB = "annulus[[%s, %s], ['3arcsec', '5arcsec']]" % \ (ra_p, dec_p) noise_annulus_LB = "annulus[[%s, %s], ['1.5arcsec', '2.5arcsec']]" % \ (ra_p, dec_p) SB_scales = [0] LB_scales = [0, 11, 33, 55] if not skip_plots: """ Image each dataset individually """ tclean_wrapper(vis=prefix+'_SB1_initcont.ms', imagename=prefix+'_SB1_initcont', scales=SB_scales, mask=mask_SB, savemodel='modelcolumn', threshold='1mJy') tclean_wrapper(vis=prefix+'_SB2_initcont.ms', imagename=prefix+'_SB2_initcont', scales=SB_scales, mask=mask_SB, savemodel='modelcolumn', threshold='1mJy') tclean_wrapper(vis=prefix+'_SB3_initcont.ms', imagename=prefix+'_SB3_initcont', scales=SB_scales, mask=mask_SB, savemodel='modelcolumn', threshold='1mJy') image_each_obs(data_params['SB3'], prefix, mask=mask_SB, scales=SB_scales, threshold='0.1mJy', interactive=False) tclean_wrapper(vis=prefix+'_LB1_initcont.ms', imagename=prefix+'_LB1_initcont', scales=LB_scales, mask=mask_LB, savemodel='modelcolumn', threshold='0.1mJy', uvtaper=['0mas','40mas','90deg']) """ Fit Gaussians to roughly estimate centers, inclinations, PAs """ fit_gaussian(prefix+'_SB1_initcont.image', region=mask_SB_p) #Peak : J2000 16h11m31.354744s -18d38m26.09769s fit_gaussian(prefix+'_SB2_initcont.image', region=mask_SB_p) #Peak : J2000 16h11m31.356286s -18d38m26.13514s fit_gaussian(prefix+'_SB3_initcont.image', region=mask_SB_p) #Peak : ICRS 16h11m31.351009s -18d38m26.25311s fit_gaussian(prefix+'_SB3_initcont_exec0.image', region=mask_SB_p) #Peak : ICRS 16h11m31.351910s -18d38m26.26593s fit_gaussian(prefix+'_SB3_initcont_exec1.image', region=mask_SB_p) #Peak : ICRS 16h11m31.352010s -18d38m26.22882s fit_gaussian(prefix+'_SB3_initcont_exec2.image', region=mask_SB_p) #Peak : ICRS 16h11m31.349837s -18d38m26.25763s fit_gaussian(prefix+'_LB1_initcont.image', region=mask_LB_p) #Peak : ICRS 16h11m31.351310s -18d38m26.24660s fit_gaussian(prefix+'_LB1_initcont.image', region=mask_LB_s) #Peak : ICRS 16h11m31.294774s -18d38m27.27574s """ There are significant misalignments that we need to fix. """ """ Define a common direction (peak of 2nd LB EB) """ common_dir = 'J2000 16h11m31.35202s -18d38m26.232956s' """ 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 16h11m31.35473s -18d38m26.09852s') fixplanets(vis=SB_shift+'.ms', field=data_params['SB1']['field'], direction=common_dir) SB_shift = prefix+'_SB2_initcont_shift' os.system('rm -rf '+SB_shift+'.ms*') fixvis(vis=prefix+'_SB2_initcont.ms', outputvis=SB_shift+'.ms', field=data_params['SB2']['field'], phasecenter='J2000 16h11m31.35629s -18d38m26.13575s') fixplanets(vis=SB_shift+'.ms', field=data_params['SB2']['field'], direction=common_dir) split_all_obs(prefix+'_SB3_initcont.ms', prefix+'_SB3_initcont_exec') SB_shift = prefix+'_SB3_initcont_exec0_shift' os.system('rm -rf '+SB_shift+'.ms*') fixvis(vis=prefix+'_SB3_initcont_exec0.ms', outputvis=SB_shift+'.ms', field=data_params['SB3']['field'], phasecenter='J2000 16h11m31.35262s -18d38m26.252276s') fixplanets(vis=SB_shift+'.ms', field=data_params['SB3']['field'], direction=common_dir) SB_shift = prefix+'_SB3_initcont_exec1_shift' os.system('rm -rf '+SB_shift+'.ms*') fixvis(vis=prefix+'_SB3_initcont_exec1.ms', outputvis=SB_shift+'.ms', field=data_params['SB3']['field'], phasecenter='J2000 16h11m31.35272s -18d38m26.215186s') fixplanets(vis=SB_shift+'.ms', field=data_params['SB3']['field'], direction=common_dir) SB_shift = prefix+'_SB3_initcont_exec2_shift' os.system('rm -rf '+SB_shift+'.ms*') fixvis(vis=prefix+'_SB3_initcont_exec2.ms', outputvis=SB_shift+'.ms', field=data_params['SB3']['field'], phasecenter='J2000 16h11m31.35055s -18d38m26.243996s') fixplanets(vis=SB_shift+'.ms', field=data_params['SB3']['field'], direction=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) """ Inspect the flux calibration. """ if not skip_plots: """ Assign rough emission geometry parameters. """ PA, incl = 111, 24 """ Export MS contents into Numpy save files """ for msfile in [prefix+'_SB1_initcont_shift.ms', prefix+'_SB2_initcont_shift.ms', prefix+'_SB3_initcont_exec0_shift.ms', prefix+'_SB3_initcont_exec1_shift.ms', prefix+'_SB3_initcont_exec2_shift.ms', prefix+'_LB1_initcont_shift.ms']: export_MS(msfile) """ Plot deprojected visibility profiles for all data together """ plot_deprojected([prefix+'_SB1_initcont_shift.vis.npz', prefix+'_SB2_initcont_shift.vis.npz', prefix+'_SB3_initcont_exec0_shift.vis.npz', prefix+'_SB3_initcont_exec1_shift.vis.npz', prefix+'_SB3_initcont_exec2_shift.vis.npz', prefix+'_LB1_initcont_shift.vis.npz'], fluxscale=[1.,1.,1.,1.,1.,1.], PA=PA, incl=incl, show_err=False) """ Correct the flux scales where appropriate. """ rescale_flux(prefix+'_SB1_initcont_shift.ms', [1.080]) rescale_flux(prefix+'_SB2_initcont_shift.ms', [1.033]) rescale_flux(prefix+'_SB3_initcont_exec1_shift.ms', [1.044]) rescale_flux(prefix+'_SB3_initcont_exec2_shift.ms', [1.040]) rescale_flux(prefix+'_LB1_initcont_shift.ms', [1.009]) """ SELF-CAL for short-baseline data """ """ Merge the SB executions back into a single MS """ SB_cont_p0 = prefix+'_SB_contp0' os.system('rm -rf %s*' % SB_cont_p0) concat(vis=[prefix+'_SB1_initcont_shift_rescaled.ms', prefix+'_SB2_initcont_shift_rescaled.ms', prefix+'_SB3_initcont_exec0_shift.ms', prefix+'_SB3_initcont_exec1_shift_rescaled.ms', prefix+'_SB3_initcont_exec2_shift_rescaled.ms'], concatvis=SB_cont_p0+'.ms', dirtol='0.1arcsec', copypointing=False) """ Initial clean """ SB_scales = [0, 9] 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) #AS205_SB_contp0.image #Beam 0.268 arcsec x 0.228 arcsec (86.44 deg) #Flux inside disk mask: 346.66 mJy #Peak intensity of source: 97.39 mJy/beam #rms: 1.69e-01 mJy/beam #Peak SNR: 576.63 """ Self-calibration parameters """ SB_contspws = '0~15' SB_refant = 'DV09@A046,DA46@A034' SB1_timerange = '2012/03/27/00:00:01~2012/03/27/23:59:59' SB2_timerange = '2012/05/04/00:00:01~2012/05/04/23:59:59' SB3_obs0_timerange = '2017/05/14/00:00:01~2017/05/14/23:59:59' SB3_obs1_timerange = '2017/05/17/00:00:01~2017/05/17/23:59:59' SB3_obs2_timerange = '2017/05/19/00:00:01~2017/05/19/23:59:59' """ 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_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p1, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB2_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p1, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB3_obs0_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p1, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB3_obs1_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p1, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB3_obs2_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) #AS205_SB_contp1.image #Beam 0.270 arcsec x 0.228 arcsec (84.95 deg) #Flux inside disk mask: 352.65 mJy #Peak intensity of source: 103.14 mJy/beam #rms: 5.90e-02 mJy/beam #Peak SNR: 1749.40 """ Second round of phase-only self-cal (short baselines only) """ SB_p2 = prefix+'_SB.p2' os.system('rm -rf '+SB_p2) gaincal(vis=SB_cont_p1+'.ms', caltable=SB_p2, gaintype='T', spw=SB_contspws, refant=SB_refant, calmode='p', solint='60s', minsnr=1.5, minblperant=4) if not skip_plots: """ Inspect gain tables """ plotcal(caltable=SB_p2, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB1_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p2, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB2_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p2, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB3_obs0_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p2, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB3_obs1_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=SB_p2, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB3_obs2_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.09mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p2+'.image', disk_mask=mask_SB_p, noise_mask=noise_annulus_SB) #AS205_SB_contp2.image #Beam 0.270 arcsec x 0.228 arcsec (84.79 deg) #Flux inside disk mask: 352.93 mJy #Peak intensity of source: 103.18 mJy/beam #rms: 5.83e-02 mJy/beam #Peak SNR: 1770.71 """ Amplitude self-cal (short baselines only) """ SB_ap = prefix+'_SB.ap' os.system('rm -rf '+SB_ap) gaincal(vis=SB_cont_p2+'.ms', caltable=SB_ap, gaintype='T', spw=SB_contspws, refant=SB_refant, calmode='ap', solint='inf', minsnr=3.0, minblperant=4, solnorm=False) if not skip_plots: """ Inspect gain tables """ plotcal(caltable=SB_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=SB1_timerange, plotrange=[0,0,0,2]) plotcal(caltable=SB_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=SB2_timerange, plotrange=[0,0,0,2]) plotcal(caltable=SB_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=SB3_obs0_timerange, plotrange=[0,0,0,2]) plotcal(caltable=SB_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=SB3_obs1_timerange, plotrange=[0,0,0,2]) plotcal(caltable=SB_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=SB3_obs2_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=mask_SB, scales=SB_scales, threshold='0.06mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_ap+'.image', disk_mask=mask_SB_p, noise_mask=noise_annulus_SB) #AS205_SB_contap.image #Beam 0.270 arcsec x 0.227 arcsec (84.89 deg) #Flux inside disk mask: 353.81 mJy #Peak intensity of source: 102.50 mJy/beam #rms: 3.10e-02 mJy/beam #Peak SNR: 3310.36 """ 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_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=mask_LB, scales=LB_scales, threshold='0.1mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p0+'.image', disk_mask=mask_LB_p, noise_mask=noise_annulus_LB) #AS205_combined_contp0.image #Beam 0.039 arcsec x 0.025 arcsec (-81.58 deg) #Flux inside disk mask: 374.69 mJy #Peak intensity of source: 4.89 mJy/beam #rms: 2.90e-02 mJy/beam #Peak SNR: 168.75 """ Self-calibration parameters """ combined_contspws = '0~19' combined_refant = 'DA61@A015,DV09@A046,DA46@A034' combined_spwmap = [0,0,0,0,4,4,4,4,8,8,8,8,12,12,12,12,16,16,16,16] LB1_obs0_timerange = '2017/09/29/00:00:01~2017/09/29/23:59:59' """ 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.06mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p1+'.image', disk_mask=mask_LB_p, noise_mask=noise_annulus_LB) #AS205_combined_contp1.image #Beam 0.039 arcsec x 0.025 arcsec (-81.58 deg) #Flux inside disk mask: 363.44 mJy #Peak intensity of source: 5.71 mJy/beam #rms: 1.96e-02 mJy/beam #Peak SNR: 291.45 """ 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.04mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p2+'.image', disk_mask=mask_LB_p, noise_mask=noise_annulus_LB) #AS205_combined_contp2.image #Beam 0.039 arcsec x 0.025 arcsec (-81.58 deg) #Flux inside disk mask: 358.54 mJy #Peak intensity of source: 6.13 mJy/beam #rms: 1.80e-02 mJy/beam #Peak SNR: 341.17 """ 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.03mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p3+'.image', disk_mask=mask_LB_p, noise_mask=noise_annulus_LB) #AS205_combined_contp3.image #Beam 0.039 arcsec x 0.025 arcsec (-81.58 deg) #Flux inside disk mask: 356.45 mJy #Peak intensity of source: 6.45 mJy/beam #rms: 1.78e-02 mJy/beam #Peak SNR: 362.7 """ 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.03mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p4+'.image', disk_mask=mask_LB_p, noise_mask=noise_annulus_LB) #AS205_combined_contp4.image #Beam 0.039 arcsec x 0.025 arcsec (-81.58 deg) #Flux inside disk mask: 356.56 mJy #Peak intensity of source: 6.61 mJy/beam #rms: 1.78e-02 mJy/beam #Peak SNR: 371.72 """ 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.03mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_ap+'.image', disk_mask=mask_LB_p, noise_mask=noise_annulus_LB) #AS205_combined_contap.image #Beam 0.039 arcsec x 0.025 arcsec (-83.55 deg) #Flux inside disk mask: 357.24 mJy #Peak intensity of source: 6.32 mJy/beam #rms: 1.64e-02 mJy/beam #Peak SNR: 384.05 """ Amplitude self-cal (all data) """ combined_ap1 = prefix+'_combined.ap1' os.system('rm -rf '+combined_ap1) gaincal(vis=combined_cont_ap+'.ms', caltable=combined_ap1, gaintype='T', combine='spw,scan', spw=combined_contspws, refant=combined_refant, calmode='ap', solint='180s', minsnr=3.0, minblperant=4, solnorm=False) if not skip_plots: """ Inspect gain tables """ plotcal(caltable=combined_ap1, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=LB1_obs0_timerange, plotrange=[0,0,0,2]) """ Apply the solutions """ applycal(vis=combined_cont_ap+'.ms', spw=combined_contspws, spwmap=combined_spwmap, gaintable=[combined_ap1], interp='linearPD', calwt=True, applymode='calonly') """ Split off a corrected MS """ combined_cont_ap1 = prefix+'_combined_contap1' os.system('rm -rf %s*' % combined_cont_ap1) split(vis=combined_cont_ap+'.ms', outputvis=combined_cont_ap1+'.ms', datacolumn='corrected') """ Image the results; check the resulting map """ tclean_wrapper(vis=combined_cont_ap1+'.ms', imagename=combined_cont_ap1, mask=mask_LB, scales=LB_scales, threshold='0.03mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_ap1+'.image', disk_mask=mask_LB_p, noise_mask=noise_annulus_LB) #AS205_combined_contap1.image #Beam 0.038 arcsec x 0.025 arcsec (-84.63 deg) #Flux inside disk mask: 358.01 mJy #Peak intensity of source: 6.15 mJy/beam #rms: 1.61e-02 mJy/beam #Peak SNR: 381.28 """ Final outputs """ """ Save the final MS """ os.system('cp -r '+combined_cont_ap1+'.ms '+prefix+'_continuum.ms') os.system('tar cvzf '+prefix+'_continuum.ms.tgz '+prefix+'_continuum.ms') """ Make a fiducial continuum image (based on experimentation) """ os.system('cp -r '+combined_cont_ap1+'.image '+prefix+'_continuum.image') exportfits(prefix+'_continuum.image', prefix+'_continuum.fits')