""" This script was written for CASA 5.1.1 Datasets calibrated (in order of date observed): SB1: 2013.1.00498.S Observed 21 July 2015 (1 execution block) LB1: 2016.1.00484.L Observed 07 September 2017 and 03 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 = 'Elias27' SB1_path = '/full_path/to_calibrated/msfile.ms' LB1_path = '/full_path/to_calibrated/msfile.ms' """ Some preliminary flagging """ flagdata(vis=SB1_path, spw='1:167~215,5:49~72,6:806~829', field='Elia_2-27') # 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': 'Elia_2-27', 'line_spws': np.array([1,5,6]), # CO 13CO, C18O SPWs 'line_freqs': np.array([2.30538e11, 2.2039868420e11, 2.1956035410e11]), }, 'LB1': {'vis' : LB1_path, 'name' : 'LB1', 'field' : 'Elias_27', '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) """ Some additional flagging """ flagmanager(vis=prefix+'_SB1_initcont.ms', mode='save', versionname='init_cal_flags', comment='Flag states immediately after initial calibration') flagdata(vis=prefix+'_SB1_initcont.ms', mode='manual', spw='1,3,5', flagbackup=False, field=data_params['SB1']['field'], scan='32', antenna='DA46') flagdata(vis=prefix+'_SB1_initcont.ms', mode='manual', spw='3,4,6', flagbackup=False, field=data_params['SB1']['field'], scan='18,32,37', antenna='DA59') flagdata(vis=prefix+'_SB1_initcont.ms', mode='manual', spw='3,4', flagbackup=False, field=data_params['SB1']['field'], scan='13,18', antenna='DV06') flagdata(vis=prefix+'_SB1_initcont.ms', mode='manual', spw='3,4', flagbackup=False, field=data_params['SB1']['field'], scan='32', antenna='DV08') flagdata(vis=prefix+'_SB1_initcont.ms', mode='manual', spw='3', flagbackup=False, field=data_params['SB1']['field'], scan='32,37', antenna='DV18') """ Define simple masks and clean scales for imaging """ mask_PA = 118 # position angle of mask in degrees mask_maj = 2.5 # semimajor axis of mask in arcsec mask_min = 1.5 # semiminor axis of mask in arcsec mask_RA = '16h26m45.027s' mask_DEC = '-24.23.08.28' 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, 25, 50] LB_scales = [0, 10, 30, 100, 200] 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.2mJy', interactive=False) image_each_obs(data_params['LB1'], prefix, mask=LB1_mask, scales=LB_scales, threshold='0.08mJy', interactive=False) """ Fit Gaussians to roughly estimate centers, inclinations, PAs """ fit_gaussian(prefix+'_SB1_initcont_exec0.image', region=SB1_mask) #Peak : J2000 16h26m45.021955s -24d23m08.25057s #PA of Gaussian component: 121.00 deg #Inclination of Gaussian component: 49.72 deg fit_gaussian(prefix+'_LB1_initcont_exec0.image', region='circle[[%s, %s], %.1farcsec]' % \ (mask_RA, mask_DEC, 0.3), dooff = True) # shrink region size to avoid bias from large-scale spirals #Peak : 16:26:45.02183, -024:23:08.272621 #PA of Gaussian component: 120.58 deg #Inclination of Gaussian component: 56.86 deg fit_gaussian(prefix+'_LB1_initcont_exec1.image', region='circle[[%s, %s], %.1farcsec]' % \ (mask_RA, mask_DEC, 0.3), dooff = True) # shrink region size to avoid bias from large-scale spirals #Peak : 16:26:45.02187, -024:23:08.274981 #PA of Gaussian component: 121.02 deg #Inclination of Gaussian component: 59.04 deg """ The emission appears to be aligned. But, we want to force everything to the same phase center to avoid confusion in the visibilities. """ """ Shift each MS to same phase center (LB1 EB1 as reference) """ common_dir = 'J2000 16h26m45.022s -024.23.08.273' shiftname = prefix+'_SB1_initcont_shift' os.system('rm -rf %s.ms' % shiftname) fixvis(vis=prefix+'_SB1_initcont.ms', outputvis=shiftname+'.ms', field=data_params['SB1']['field'], phasecenter='J2000 16h26m45.021955s -24d23m08.25057s') fixplanets(vis=shiftname+'.ms', field=data_params['SB1']['field'], direction=common_dir) shiftname = prefix+'_LB1_initcont_shift' os.system('rm -rf %s.ms' % shiftname) fixvis(vis=prefix+'_LB1_initcont.ms', outputvis=shiftname+'.ms', field=data_params['LB1']['field'], phasecenter='ICRS 16h26m45.021309s -24d23m08.28623s') fixplanets(vis=shiftname+'.ms', field=data_params['LB1']['field'], direction=common_dir) """ Now that everything is aligned, we inspect the flux calibration. """ split_all_obs(prefix+'_LB1_initcont_shift.ms', prefix+'_LB1_initcont_shift_exec') if not skip_plots: """ Assign rough emission geometry parameters. """ PA, incl = 121., 57. """ Export MS contents into Numpy save files """ for msfile in [prefix+'_SB1_initcont_shift.ms', prefix+'_LB1_initcont_shift_exec0.ms', prefix+'_LB1_initcont_shift_exec1.ms']: export_MS(msfile) """ Plot deprojected visibility profiles for all data together """ plot_deprojected([prefix+'_SB1_initcont_shift.vis.npz', prefix+'_LB1_initcont_shift_exec0.vis.npz', prefix+'_LB1_initcont_shift_exec1.vis.npz'], fluxscale=[1., 1., 1.], PA=PA, incl=incl, show_err=False) """ Now inspect offsets by comparing against a reference """ estimate_flux_scale(reference=prefix+'_SB1_initcont_shift.vis.npz', comparison=prefix+'_LB1_initcont_shift_exec0.vis.npz', incl=incl, PA=PA) #There is discrepancy, but this turns out to be purely phase noise estimate_flux_scale(reference=prefix+'_SB1_initcont_shift.vis.npz', comparison=prefix+'_LB1_initcont_shift_exec1.vis.npz', incl=incl, PA=PA) #More phase noise. """ SELF-CAL for short-baseline data """ """ Copy the single SB execution """ SB_cont_p0 = prefix+'_SB_contp0' os.system('rm -rf %s*' % SB_cont_p0) os.system('cp -r '+prefix+'_SB1_initcont_shift.ms '+SB_cont_p0+'.ms') """ Set up a clean mask """ mask_ra = '16h26m45.022s' mask_dec = '-024.23.08.2733' 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.2mJy', 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) #Elias27_SB_contp0.image #Beam 0.233 arcsec x 0.194 arcsec (62.94 deg) #Flux inside disk mask: 299.04 mJy #Peak intensity of source: 23.99 mJy/beam #rms: 8.38e-02 mJy/beam #Peak SNR: 286.10 """ Self-calibration parameters """ SB_contspws = '0~7' SB_refant = 'DV09@A007' SB1_obs0_timerange = '2015/07/20/00~2015/07/22/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_obs0_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.15mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p1+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias27_SB_contp1.image #Beam 0.233 arcsec x 0.194 arcsec (62.59 deg) #Flux inside disk mask: 310.31 mJy #Peak intensity of source: 26.44 mJy/beam #rms: 5.61e-02 mJy/beam #Peak SNR: 471.18 """ 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_obs0_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.15mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p2+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias27_SB_contp2.image #Beam 0.234 arcsec x 0.196 arcsec (62.94 deg) #Flux inside disk mask: 310.55 mJy #Peak intensity of source: 27.27 mJy/beam #rms: 5.58e-02 mJy/beam #Peak SNR: 488.54 """ 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]) """ 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.12mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_ap+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias27_SB_contap.image #Beam 0.233 arcsec x 0.196 arcsec (63.57 deg) #Flux inside disk mask: 306.69 mJy #Peak intensity of source: 27.00 mJy/beam #rms: 4.57e-02 mJy/beam #Peak SNR: 591.08 """ 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.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', imsize=4000, uvtaper='.03arcsec', savemodel='modelcolumn') estimate_SNR(combined_cont_p0+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Beam 0.053 arcsec x 0.041 arcsec (73.57 deg) #Flux inside disk mask: 330.57 mJy #Peak intensity of source: 4.43 mJy/beam #rms: 1.58e-02 mJy/beam #Peak SNR: 280.14 """ Self-calibration parameters """ combined_contspws = '0~15' combined_refant = 'DV09@A007,DA61@A015,DA50@A108' combined_spwmap = [0,0,0,0,0,0,0,0,8,8,8,8,12,12,12,12] LB1_obs0_timerange = '2017/09/07/00~2017/09/08/00' LB2_obs1_timerange = '2017/10/02/00~2017/10/05/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, scales=LB_scales, threshold='0.05mJy', imsize=4000, uvtaper='.03arcsec', savemodel='modelcolumn') estimate_SNR(combined_cont_p1+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias27_combined_contp1.image #Beam 0.053 arcsec x 0.041 arcsec (73.57 deg) #Flux inside disk mask: 330.74 mJy #Peak intensity of source: 4.46 mJy/beam #rms: 1.55e-02 mJy/beam #Peak SNR: 286.95 """ 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, scales=LB_scales, threshold='0.05mJy', imsize=4000, uvtaper='.03arcsec', savemodel='modelcolumn') estimate_SNR(combined_cont_p2+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias27_combined_contp2.image #Beam 0.053 arcsec x 0.041 arcsec (73.57 deg) #Flux inside disk mask: 330.89 mJy #Peak intensity of source: 4.53 mJy/beam #rms: 1.55e-02 mJy/beam #Peak SNR: 293.27 """ 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, scales=LB_scales, threshold='0.05mJy', imsize=4000, uvtaper='.03arcsec', savemodel='modelcolumn') estimate_SNR(combined_cont_p3+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias27_combined_contp3.image #Beam 0.053 arcsec x 0.041 arcsec (73.57 deg) #Flux inside disk mask: 330.87 mJy #Peak intensity of source: 4.66 mJy/beam #rms: 1.54e-02 mJy/beam #Peak SNR: 303.14 """ 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]) 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.05mJy', imsize=4000, uvtaper='.03arcsec', savemodel='modelcolumn') estimate_SNR(combined_cont_p4+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias27_combined_contp4.image #Beam 0.053 arcsec x 0.041 arcsec (73.57 deg) #Flux inside disk mask: 331.03 mJy #Peak intensity of source: 4.74 mJy/beam #rms: 1.54e-02 mJy/beam #Peak SNR: 308.03 """ 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='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_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=common_mask, scales=LB_scales, threshold='0.05mJy', imsize=4000, uvtaper='.03arcsec', savemodel='modelcolumn') estimate_SNR(combined_cont_ap+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias27_combined_contap.image #Beam 0.054 arcsec x 0.041 arcsec (73.62 deg) #Flux inside disk mask: 329.29 mJy #Peak intensity of source: 4.64 mJy/beam #rms: 1.41e-02 mJy/beam #Peak SNR: 329.72 """ 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, 20, 50, 100, 200] tclean_wrapper(vis=combined_cont_ap+'.ms', imagename=prefix+'_continuum', mask=common_mask, scales=scales, threshold='0.08mJy', imsize=4000, uvtaper=['0.04arcsec','0.02arcsec', '173.2deg']) exportfits(prefix+'_continuum.image', prefix+'_continuum.fits')