""" 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 25 September 2017 and 03 October 2017 (2 execution blocks) reducer: S. Andrews """ """ Starting matter """ import os execfile('/pool/asha0/SCIENCE/p484/sa_work/reduction_utils.py') """ Input for loading data """ prefix = 'Elias24' 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': 'Elias_24', '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_24', '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='Elias24_SB1_initcont.ms', mode='save', versionname='init_cal_flags', comment='Flag states immediately after initial calibration') flagdata(vis='Elias24_SB1_initcont.ms', mode='manual', spw='3,4,6', flagbackup=False, field=data_params['SB1']['field'], timerange='2015/07/21/22:45:00~2015/07/21/23:10:00', antenna='DA41') flagdata(vis='Elias24_SB1_initcont.ms', mode='manual', spw='1,3,5', flagbackup=False, field=data_params['SB1']['field'], scan='30', antenna='DA46') flagdata(vis='Elias24_SB1_initcont.ms', mode='manual', spw='3,4,6', flagbackup=False, field=data_params['SB1']['field'], scan='30,35', antenna='DA59') flagdata(vis='Elias24_SB1_initcont.ms', mode='manual', spw='3,4', flagbackup=False, field=data_params['SB1']['field'], scan='30', antenna='DV08') flagdata(vis='Elias24_SB1_initcont.ms', mode='manual', spw='3', flagbackup=False, field=data_params['SB1']['field'], scan='30,35', antenna='DV18') """ Define simple masks and clean scales for imaging """ mask_PA = 52 # position angle of mask in degrees mask_maj = 1.6 # semimajor axis of mask in arcsec mask_min = 1.4 # semiminor axis of mask in arcsec mask_ra = '16h26m24.08s' mask_dec = '-24.16.13.89' 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] LB_scales = [0, 5, 20, 60, 180] 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.1mJy', interactive=False) """ Fit Gaussians to roughly estimate centers, inclinations, PAs """ fit_gaussian(prefix+'_SB1_initcont_exec0.image', region=SB1_mask) #Peak : J2000 16h26m24.078792s -24d16m13.84764s #PA of Gaussian component: 51.52 deg #Inclination of Gaussian component: 27.37 deg fit_gaussian(prefix+'_LB1_initcont_exec0.image', region=LB1_mask) #Peak : ICRS 16h26m24.077753s -24d16m13.88292s #PA of Gaussian component: 46.63 deg #Inclination of Gaussian component: 30.39 deg fit_gaussian(prefix+'_LB1_initcont_exec1.image', region=LB1_mask) #Peak : ICRS 16h26m24.078146s -24d16m13.88695s #PA of Gaussian component: 43.57 deg #Inclination of Gaussian component: 30.34 deg """ The emission appears to be aligned. But, we want to force everything to the same phase center to avoid confusion in the visibilities. """ """ Split out individual MSs for each execution """ split_all_obs(prefix+'_SB1_initcont.ms', prefix+'_SB1_initcont_exec') split_all_obs(prefix+'_LB1_initcont.ms', prefix+'_LB1_initcont_exec') """ Shift each MS to same phase center (LB1 EB1 as reference) """ SB0_shift = prefix+'_SB1_initcont_exec0_shift.ms' os.system('rm -rf '+SB0_shift+'*') fixvis(vis=prefix+'_SB1_initcont_exec0.ms', outputvis=SB0_shift, field=data_params['SB1']['field'], phasecenter='ICRS 16h26m24.070000s -24d16m13.50000s') """ Now that everything is aligned, we inspect the flux calibration. """ if not skip_plots: """ Assign rough emission geometry parameters and centers. """ PA, incl = 45., 30. phasecenter = au.radec2deg('16:26:24.07000, -24.16.13.50000') peakpos = au.radec2deg('16:26:24.078, -24.16.13.883') offsets = au.angularSeparation(peakpos[0], peakpos[1], phasecenter[0], phasecenter[1], True) offx, offy = 3600.*offsets[3], 3600.*offsets[2] """ Export MS contents into Numpy save files """ for msfile in [prefix+'_SB1_initcont_exec0_shift.ms', prefix+'_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_shift.vis.npz', prefix+'_LB1_initcont_exec0.vis.npz', prefix+'_LB1_initcont_exec1.vis.npz'], fluxscale=[1.0, 1.0, 1.0], offx=offx, offy=offy, PA=PA, incl=incl, show_err=False) """ Now inspect offsets by comparing against a reference """ estimate_flux_scale(reference=prefix+'_LB1_initcont_exec0_shift.vis.npz', comparison=prefix+'_SB1_initcont_exec0_shift.vis.npz', offx=offx, offy=offy, incl=incl, PA=PA) estimate_flux_scale(reference=prefix+'_LB1_initcont_exec0_shift.vis.npz', comparison=prefix+'_LB1_initcont_exec1.vis.npz', offx=offx, offy=offy, incl=incl, PA=PA) # No flux rescaling is necessary. """ 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_exec0_shift.ms '+SB_cont_p0+'.ms') """ Set up a clean mask """ mask_PA = 45. mask_ra = '16h26m24.078s' mask_dec = '-24.16.13.883' common_mask = 'ellipse[[%s, %s], [%.1farcsec, %.1farcsec], %.1fdeg]' % \ (mask_ra, mask_dec, mask_maj, mask_min, mask_PA) """ Initial clean """ tclean_wrapper(vis=SB_cont_p0+'.ms', imagename=SB_cont_p0, mask=common_mask, scales=SB_scales, threshold='0.15mJy', savemodel='modelcolumn') """ Define a noise annulus, measure the peak SNR in map """ noise_annulus = "annulus[[%s, %s],['%.2farcsec', '4.25arcsec']]" % \ (mask_ra, mask_dec, 1.1*mask_maj) estimate_SNR(SB_cont_p0+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias24_SB_contp0.image #Beam 0.234 arcsec x 0.194 arcsec (62.68 deg) #Flux inside disk mask: 346.55 mJy #Peak intensity of source: 39.54 mJy/beam #rms: 1.29e-01 mJy/beam #Peak SNR: 307.27 """ Self-calibration parameters """ SB_contspws = '0~7' SB_refant = 'DV09, DV08, DV06' 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) #Elias24_SB_contp1.image #Beam 0.234 arcsec x 0.194 arcsec (62.33 deg) #Flux inside disk mask: 360.15 mJy #Peak intensity of source: 43.23 mJy/beam #rms: 5.87e-02 mJy/beam #Peak SNR: 736.83 """ 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.1mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_p2+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias24_SB_contp2.image #Beam 0.235 arcsec x 0.194 arcsec (62.35 deg) #Flux inside disk mask: 360.78 mJy #Peak intensity of source: 43.93 mJy/beam #rms: 5.76e-02 mJy/beam #Peak SNR: 762.80 """ 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.07mJy', savemodel='modelcolumn') estimate_SNR(SB_cont_ap+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias24_SB_contap.image #Beam 0.232 arcsec x 0.192 arcsec (62.90 deg) #Flux inside disk mask: 360.31 mJy #Peak intensity of source: 43.16 mJy/beam #rms: 4.73e-02 mJy/beam #Peak SNR: 911.65 """ 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.ms', prefix+'_LB1_initcont_exec1.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.1mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p0+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias24_combined_contp0.image #Beam 0.045 arcsec x 0.030 arcsec (86.49 deg) #Flux inside disk mask: 367.59 mJy #Peak intensity of source: 4.67 mJy/beam #rms: 1.83e-02 mJy/beam #Peak SNR: 255.28 """ Self-calibration parameters """ combined_contspws = '0~15' combined_refant = 'DV24@A090,DA61@A015,DV09@A007,DV08@A008,DV06@A014' combined_spwmap = [0,0,0,0,0,0,0,0,8,8,8,8,12,12,12,12] LB1_obs0_timerange = '2017/09/24/00~2017/09/26/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='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.1mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p1+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias24_combined_contp1.image #Beam 0.045 arcsec x 0.030 arcsec (86.49 deg) #Flux inside disk mask: 367.09 mJy #Peak intensity of source: 4.58 mJy/beam #rms: 1.81e-02 mJy/beam #Peak SNR: 253.03 # ** map quality still 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.04mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p2+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias24_combined_contp2.image #Beam 0.045 arcsec x 0.030 arcsec (86.49 deg) #Flux inside disk mask: 367.30 mJy #Peak intensity of source: 4.63 mJy/beam #rms: 1.72e-02 mJy/beam #Peak SNR: 269.92 # ** 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='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.04mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p3+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias24_combined_contp3.image #Beam 0.045 arcsec x 0.030 arcsec (86.49 deg) #Flux inside disk mask: 367.44 mJy #Peak intensity of source: 4.70 mJy/beam #rms: 1.70e-02 mJy/beam #Peak SNR: 276.39 """ Fourth round of phase-only self-cal (all data) """ combined_p4 = prefix+'_combined.p4' os.system('rm -rf '+combined_p4) gaincal(vis=combined_cont_p3+'.ms', caltable=combined_p4, gaintype='T', combine='spw,scan', spw=combined_contspws, refant=combined_refant, calmode='p', solint='60s', minsnr=1.5, minblperant=4) if not skip_plots: """ Inspect gain tables """ plotcal(caltable=combined_p4, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=LB1_obs0_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=combined_p4, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=LB2_obs1_timerange, plotrange=[0,0,-180,180]) """ Apply the solutions """ applycal(vis=combined_cont_p3+'.ms', spw=combined_contspws, spwmap=combined_spwmap, gaintable=[combined_p4], interp='linearPD', calwt=True, applymode='calonly') """ Split off a corrected MS """ combined_cont_p4 = prefix+'_combined_contp4' os.system('rm -rf %s*' % combined_cont_p4) split(vis=combined_cont_p3+'.ms', outputvis=combined_cont_p4+'.ms', datacolumn='corrected') """ Image the results; check the resulting map """ tclean_wrapper(vis=combined_cont_p4+'.ms', imagename=combined_cont_p4, mask=common_mask, scales=LB_scales, threshold='0.03mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p4+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias24_combined_contp4.image #Beam 0.045 arcsec x 0.030 arcsec (86.49 deg) #Flux inside disk mask: 367.61 mJy #Peak intensity of source: 4.88 mJy/beam #rms: 1.69e-02 mJy/beam #Peak SNR: 288.42 """ Fifth round of phase-only self-cal (all data) """ combined_p5 = prefix+'_combined.p5' os.system('rm -rf '+combined_p5) gaincal(vis=combined_cont_p4+'.ms', caltable=combined_p5, gaintype='T', combine='spw,scan', spw=combined_contspws, refant=combined_refant, calmode='p', solint='30s', minsnr=1.5, minblperant=4) if not skip_plots: """ Inspect gain tables """ plotcal(caltable=combined_p5, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=LB1_obs0_timerange, plotrange=[0,0,-180,180]) plotcal(caltable=combined_p5, xaxis='time', yaxis='phase', subplot=441, iteration='antenna', timerange=LB2_obs1_timerange, plotrange=[0,0,-180,180]) """ Apply the solutions """ applycal(vis=combined_cont_p4+'.ms', spw=combined_contspws, spwmap=combined_spwmap, gaintable=[combined_p5], interp='linearPD', calwt=True, applymode='calonly') """ Split off a corrected MS """ combined_cont_p5 = prefix+'_combined_contp5' os.system('rm -rf %s*' % combined_cont_p5) split(vis=combined_cont_p4+'.ms', outputvis=combined_cont_p5+'.ms', datacolumn='corrected') """ Image the results; check the resulting map """ tclean_wrapper(vis=combined_cont_p5+'.ms', imagename=combined_cont_p5, mask=common_mask, scales=LB_scales, threshold='0.03mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_p5+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias24_combined_contp5.image #Beam 0.045 arcsec x 0.030 arcsec (86.49 deg) #Flux inside disk mask: 367.67 mJy #Peak intensity of source: 4.95 mJy/beam #rms: 1.67e-02 mJy/beam #Peak SNR: 295.82 """ 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.025mJy', savemodel='modelcolumn') estimate_SNR(combined_cont_ap+'.image', disk_mask=common_mask, noise_mask=noise_annulus) #Elias24_combined_contap.image #Beam 0.041 arcsec x 0.028 arcsec (84.09 deg) #Flux inside disk mask: 369.28 mJy #Peak intensity of source: 4.35 mJy/beam #rms: 1.49e-02 mJy/beam #Peak SNR: 290.99 """ 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', robust=0., uvtaper=['.035arcsec', '.01arcsec', '166deg']) exportfits(prefix+'_continuum.image', prefix+'_continuum.fits')