""" 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 17 September 2017 and 10 October 2017 (2 execution blocks) reducer: J. Huang """ """ Starting matter """ import os execfile('reduction_utils.py') skip_plots = True # if True, can run script non-interactively """ Input for loading data """ prefix = 'DoAr33' 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': 'DoAr_33', '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' : 'DoAr_33', 'line_spws': np.array([3,7]), # CO SPWs 'line_freqs': np.array([2.30538e11, 2.30538e11]), } } """ Check data (various options here; an example) """ if not skip_plots: for i in data_params.keys(): plotms(vis=data_params[i]['vis'], xaxis='channel', yaxis='amplitude', field=data_params[i]['field'], ydatacolumn='data', avgtime='1e8', avgscan=True, avgbaseline=True, iteraxis='spw') """ Identify 50 km/s-wide region containing CO emission; then flag that and do a spectral average to a pseudo-continuum MS """ for i in data_params.keys(): flagchannels_string = get_flagchannels(data_params[i], prefix, velocity_range=np.array([-15, 25])) avg_cont(data_params[i], prefix, flagchannels=flagchannels_string) """ Define simple masks and clean scales for imaging """ mask_pa = 79 # position angle of mask in degrees mask_maj = 0.7 # semimajor axis of mask in arcsec mask_min = 0.5 # semiminor axis of mask in arcsec mask_ra = '16h27m39.00s' mask_dec = '-23.58.19.17' SB1_mask = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ (mask_ra, mask_dec, mask_maj, mask_min, mask_pa) LB1_mask = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ (mask_ra, mask_dec, mask_maj, mask_min, mask_pa) SB_scales = [0, 5, 10] LB_scales = [0, 5, 20, 60, 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.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 16h27m39.004011s -23d58m19.18062s fit_gaussian(prefix+'_SB1_initcont_exec1.image', region=SB1_mask) #Peak : ICRS 16h27m39.002326s -23d58m19.16918s fit_gaussian(prefix+'_SB1_initcont_exec2.image', region=SB1_mask) #Peak : ICRS 16h27m39.001420s -23d58m19.17448s fit_gaussian(prefix+'_LB1_initcont_exec0.image', region=LB1_mask) #Peak : ICRS 16h27m39.003192s -23d58m19.18784s fit_gaussian(prefix+'_LB1_initcont_exec1.image', region=LB1_mask) #Peak : ICRS 16h27m39.003034s -23d58m19.18898s #PA of Gaussian component: 83.14 deg #Inclination of Gaussian component: 42.62 deg """ Note that there are two point sources found in the SB maps, located about 5 arcsec almost due E (about 0.5 mJy), and a little more than 7" to the SW (about 0.3 mJy). Both of these fall outside the default LB image size. """ """ The emission centers are slightly misaligned. So we split out the individual executions, shift the peaks to the phase center, and reassign the phase centers to a common direction. (Not in script: check that these shifts do what you think they should by re-imaging! They do for us.) """ """ 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') """ Define a common direction (peak of 2nd LB EB) """ common_dir = 'J2000 16h27m39.00363s -023.58.19.176191' """ Shift each MS so emission center is at same phase center """ 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 16h27m39.004011s -23d58m19.18062s') fixplanets(vis=SB0_shift, field=data_params['SB1']['field'], direction=common_dir) 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 16h27m39.002326s -23d58m19.16918s') fixplanets(vis=SB1_shift, field=data_params['SB1']['field'], direction=common_dir) SB2_shift = prefix+'_SB1_initcont_exec2_shift.ms' os.system('rm -rf '+SB2_shift+'*') fixvis(vis=prefix+'_SB1_initcont_exec2.ms', outputvis=SB2_shift, field=data_params['SB1']['field'], phasecenter='ICRS 16h27m39.001420s -23d58m19.17448s') fixplanets(vis=SB2_shift, field=data_params['SB1']['field'], direction=common_dir) 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 16h27m39.003192s -23d58m19.18784s') fixplanets(vis=LB0_shift, field=data_params['LB1']['field'], direction=common_dir) LB1_shift = prefix+'_LB1_initcont_exec1_shift.ms' os.system('rm -rf '+LB1_shift+'*') fixvis(vis=prefix+'_LB1_initcont_exec1.ms', outputvis=LB1_shift, field=data_params['LB1']['field'], phasecenter='ICRS 16h27m39.003034s -23d58m19.18898') fixplanets(vis=LB1_shift, field=data_params['LB1']['field'], direction=common_dir) """ Inspect the flux calibration. """ if not skip_plots: """ Assign rough emission geometry parameters. """ PA, incl = 83, 42 """ Export MS contents into Numpy save files """ for msfile in [prefix+'_SB1_initcont_exec0_shift.ms', prefix+'_SB1_initcont_exec1_shift.ms', prefix+'_SB1_initcont_exec2_shift.ms', prefix+'_LB1_initcont_exec0_shift.ms', prefix+'_LB1_initcont_exec1_shift.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+'_SB1_initcont_exec2_shift.vis.npz', prefix+'_LB1_initcont_exec0_shift.vis.npz', prefix+'_LB1_initcont_exec1_shift.vis.npz'], fluxscale=[1.0, 1.0, 1.0, 1.0, 1.0], PA=PA, incl=incl, show_err=False) # small, but non-negligible discrepancies; use SB1xEB1 as reference """ Now inspect offsets by comparing against a reference """ estimate_flux_scale(reference=prefix+'_SB1_initcont_exec1_shift.vis.npz', comparison=prefix+'_SB1_initcont_exec0_shift.vis.npz', offx=offx, offy=offy, incl=incl, PA=PA) #The ratio of comparison : reference is 0.94202 #The scaling factor for gencal is 0.971 for your comparison measurement estimate_flux_scale(reference=prefix+'_SB1_initcont_exec1_shift.vis.npz', comparison=prefix+'_SB1_initcont_exec2_shift.vis.npz', offx=offx, offy=offy, incl=incl, PA=PA) #The ratio of comparison : reference is 0.98665 #The scaling factor for gencal is 0.993 for your comparison measurement estimate_flux_scale(reference=prefix+'_SB1_initcont_exec1_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.02093 #The scaling factor for gencal is 1.010 for your comparison measurement estimate_flux_scale(reference=prefix+'_SB1_initcont_exec1_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.92562 #The scaling factor for gencal is 0.962 for your comparison measurement """ Correct the flux scales where appropriate. """ rescale_flux(prefix+'_SB1_initcont_exec0_shift.ms', [0.971]) rescale_flux(prefix+'_LB1_initcont_exec1_shift.ms', [0.962]) """ 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_rescaled.ms', prefix+'_SB1_initcont_exec1_shift.ms', prefix+'_SB1_initcont_exec2_shift.ms'], concatvis=SB_cont_p0+'.ms', dirtol='0.1arcsec', copypointing=False) """ Set up a clean mask """ mask_maj = 0.70 mask_min = 0.52 mask_pa = 83 mask_ra = '16h27m39.003s' mask_dec = '-23.58.19.189' common_mask = ['ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ (mask_ra, mask_dec, mask_maj, mask_min, mask_pa), 'circle[[16h27m39.363s, -23d58m19.017s], 0.4arcsec]', 'circle[[16h27m38.845s, -23d58m26.047s], 0.4arcsec]'] disk_mask = common_mask[0] # note the common mask includes the two outlier point sources """ 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=disk_mask, noise_mask=noise_annulus) #DoAr33_SB_contp0.image #Beam 0.268 arcsec x 0.225 arcsec (87.12 deg) #Flux inside disk mask: 34.37 mJy #Peak intensity of source: 20.99 mJy/beam #rms: 5.14e-02 mJy/beam #Peak SNR: 408.62 """ Self-calibration parameters """ SB_contspws = '0~11' SB_refant = 'DA46@A034' SB1_obs0_timerange = '2017/05/13/00~2017/05/15/00' SB1_obs1_timerange = '2017/05/15/00~2017/05/18/00' SB1_obs2_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='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]) plotcal(caltable=SB_p1, xaxis='time', yaxis='phase',subplot=441, iteration='antenna', timerange=SB1_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=common_mask, scales=SB_scales, threshold='0.1mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p1+'.image', disk_mask=disk_mask, noise_mask=noise_annulus) #DoAr33_SB_contp1.image #Beam 0.269 arcsec x 0.225 arcsec (87.08 deg) #Flux inside disk mask: 34.74 mJy #Peak intensity of source: 22.48 mJy/beam #rms: 3.12e-02 mJy/beam #Peak SNR: 719.39 """ 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]) plotcal(caltable=SB_p2, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=SB1_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=common_mask, scales=SB_scales, threshold='0.07mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p2+'.image', disk_mask=disk_mask, noise_mask=noise_annulus) #DoAr33_SB_contp2.image #Beam 0.269 arcsec x 0.225 arcsec (86.74 deg) #Flux inside disk mask: 34.81 mJy #Peak intensity of source: 22.66 mJy/beam #rms: 3.13e-02 mJy/beam #Peak SNR: 724.49 """ 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]) plotcal(caltable=SB_ap, xaxis='time', yaxis='amp', subplot=441, iteration='antenna', timerange=SB1_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=common_mask, scales=SB_scales, threshold='0.07mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_ap+'.image', disk_mask=disk_mask, noise_mask=noise_annulus) #DoAr33_SB_contap.image #Beam 0.270 arcsec x 0.224 arcsec (86.54 deg) #Flux inside disk mask: 34.86 mJy #Peak intensity of source: 22.63 mJy/beam #rms: 2.86e-02 mJy/beam #Peak SNR: 790.06 """ 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_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=common_mask, imsize=5600, scales=LB_scales, threshold='0.05mJy', gain=0.1, savemodel='modelcolumn') # increased imsize to make sure we image those outlier point sources estimate_SNR(combined_cont_p0+'.image', disk_mask=disk_mask, noise_mask=noise_annulus) #DoAr33_combined_contp0.image #Beam 0.043 arcsec x 0.024 arcsec (79.86 deg) #Flux inside disk mask: 35.08 mJy #Peak intensity of source: 1.78 mJy/beam #rms: 1.42e-02 mJy/beam #Peak SNR: 125.61 """ Self-calibration parameters """ combined_contspws = '0~19' combined_refant = 'DA61@A015,DV25@A087,DA46@A034,DA51@A124' 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/16/00~2017/09/19/00' LB2_obs1_timerange = '2017/10/09/00~2017/10/11/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]) 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, imsize=5600, scales=LB_scales, threshold='0.05mJy', gain=0.1, savemodel='modelcolumn') estimate_SNR(combined_cont_p1+'.image', disk_mask=disk_mask, noise_mask=noise_annulus) #DoAr33_combined_contp1.image #Beam 0.043 arcsec x 0.024 arcsec (79.97 deg) #Flux inside disk mask: 34.30 mJy #Peak intensity of source: 1.82 mJy/beam #rms: 1.42e-02 mJy/beam #Peak SNR: 128.49 # ** map SNR improvement negligible, but quality improves ** """ 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]) 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, imsize=5600, scales=LB_scales, threshold='0.05mJy', gain=0.1, savemodel='modelcolumn') estimate_SNR(combined_cont_p2+'.image', disk_mask=disk_mask, noise_mask=noise_annulus) #DoAr33_combined_contp2.image #Beam 0.043 arcsec x 0.024 arcsec (79.97 deg) #Flux inside disk mask: 34.05 mJy #Peak intensity of source: 1.86 mJy/beam #rms: 1.42e-02 mJy/beam #Peak SNR: 131.32 # ** map quality still improving ** """ 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]) 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, imsize=5600, scales=LB_scales, threshold='0.05mJy', gain=0.1, savemodel='modelcolumn') estimate_SNR(combined_cont_p3+'.image', disk_mask=disk_mask, noise_mask=noise_annulus) #DoAr33_combined_contp3.image #Beam 0.043 arcsec x 0.024 arcsec (79.97 deg) #Flux inside disk mask: 33.89 mJy #Peak intensity of source: 2.01 mJy/beam #rms: 1.41e-02 mJy/beam #Peak SNR: 142.88 # Additional phase-only cal on shorter intervals or amp self-cal attempts are # not helpful. """ Final outputs """ """ Save the final MS """ os.system('cp -r '+combined_cont_p3+'.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_p3+'.ms', imagename=prefix+'_continuum', mask=disk_mask, scales=LB_scales, threshold='0.06mJy', robust=0., uvtaper=['0.02arcsec','0.01arcsec','-13deg']) exportfits(prefix+'_continuum.image', prefix+'_continuum.fits')