""" This script was written for CASA 5.1.1 Datasets calibrated (in order of date observed): SB1: 2013.1.00226.S Observed 02 July 2014 (1 execution block) SB2: 2013.1.00226.S Observed 17 July 2014 (1 execution block) SB3: 2015.1.00486.S Observed 22 September 2016 (1 execution block) SB4: 2015.1.00486.S Observed 26 September 2016 (1 execution block) SB5: 2016.1.00484.L Observed 09 May 2017 (1 execution block) LB1: 2016.1.00484.L Observed 07 September 2017 and 19 September 2017 (2 execution blocks) reducer: V. Guzman """ """ Starting matter """ import os execfile('reduction_utils.py') """ Input for loading data """ prefix = 'AS209' SB1_path = '/full_path/to_calibrated/msfile.ms' # Note we only use the SB1 upper sideband, with frequencies near the LP data SB2_path = '/full_path/to_calibrated/msfile.ms' # Note we only use the SB2 lower sideband, with frequencies near the LP data SB3_path = '/full_path/to_calibrated/msfile.ms' # Note we only use the SB3 upper sideband, with frequencies near the LP data SB4_path = '/full_path/to_calibrated/msfile.ms' # Note we only use the SB4 upper sideband, with frequencies near the LP data SB5_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': 'as_209', 'line_spws': np.array([15]), # CO SPWs 'line_freqs': np.array([2.30538e11]), 'cont_spws': '14,15,16', 'width_array': [960,960,960], }, 'SB2': {'vis' : SB2_path, 'name' : 'SB2', 'field': 'as_209', 'line_spws': np.array([]), # CO SPWs 'line_freqs': np.array([]), 'cont_spws': '12,13,14,15,16,17', 'width_array': [1920,1920,960,960,960,960], }, 'SB3': {'vis' : SB3_path, 'name' : 'SB3', 'field': 'AS_209', 'line_spws': np.array([1]), # CO SPWs 'line_freqs': np.array([2.30538e11]), 'cont_spws': '0,1', 'width_array': [8,480], }, 'SB4': {'vis' : SB4_path, 'name' : 'SB4', 'field': 'AS_209', 'line_spws': np.array([1]), # CO SPWs 'line_freqs': np.array([2.30538e11]), 'cont_spws': '0,1', 'width_array': [8,480], }, 'SB5': {'vis' : SB5_path, 'name' : 'SB5', 'field': 'AS_209', 'line_spws': np.array([0]), # CO SPWs 'line_freqs': np.array([2.30538e11]), 'cont_spws': None, 'width_array': None, }, 'LB1': {'vis' : LB1_path, 'name' : 'LB1', 'field' : 'AS_209', 'line_spws': np.array([3,7]), # CO SPWs 'line_freqs': np.array([2.30538e11, 2.30538e11]), 'cont_spws': None, 'width_array': None, } } """ Split out the SB data, with the continuum SPWs averaged down """ os.system('rm -rf %s*' % prefix+'_SB1_lines') split(vis=data_params['SB1']['vis'], field=data_params['SB1']['field'], spw='14~16', outputvis=prefix+'_SB1_lines.ms', width=[960,1,960], datacolumn='data', intent='OBSERVE_TARGET#ON_SOURCE', keepflags=False) os.system('rm -rf %s*' % prefix+'_SB2_lines') split(vis=data_params['SB2']['vis'], field=data_params['SB2']['field'], spw='12~17', outputvis=prefix+'_SB2_lines.ms', width=[1920,1920,960,960,960,960], datacolumn='data', intent='OBSERVE_TARGET#ON_SOURCE', keepflags=False) os.system('rm -rf %s*' % prefix+'_SB3_lines') split(vis=data_params['SB3']['vis'], field=data_params['SB3']['field'], spw='0,1', outputvis=prefix+'_SB3_lines.ms', width=[128,1], datacolumn='data', intent='OBSERVE_TARGET#ON_SOURCE', keepflags=False) os.system('rm -rf %s*' % prefix+'_SB4_lines') split(vis=data_params['SB4']['vis'], field=data_params['SB4']['field'], spw='0,1', outputvis=prefix+'_SB4_lines.ms', width=[128,1], datacolumn='data', intent='OBSERVE_TARGET#ON_SOURCE', keepflags=False) os.system('rm -rf %s*' % prefix+'_SB5_lines') split(vis=data_params['SB5']['vis'], field=data_params['SB5']['field'], spw='0~3', outputvis=prefix+'_SB5_lines.ms', width=[1,128,128,128], datacolumn='data', intent='OBSERVE_TARGET#ON_SOURCE', keepflags=False) os.system('rm -rf %s*' % prefix+'_LB1_lines_exec0') split(vis=data_params['LB1']['vis'], field=data_params['LB1']['field'], spw='0~3', outputvis=prefix+'_LB1_lines_exec0.ms', width=[128,128,128,1], timebin='6s', datacolumn='data', intent='OBSERVE_TARGET#ON_SOURCE', keepflags=False) os.system('rm -rf %s*' % prefix+'_LB1_lines_exec1') split(vis=data_params['LB1']['vis'], field=data_params['LB1']['field'], spw='4~7', outputvis=prefix+'_LB1_lines_exec1.ms', width=[128,128,128,1], timebin='6s', datacolumn='data', intent='OBSERVE_TARGET#ON_SOURCE', keepflags=False) """ Apply same shifts and re-scalings as for continuum """ common_dir = 'J2000 16h49m15.29463s -014.22.09.048165' SB1_shift = prefix+'_SB1_lines_shift.ms' os.system('rm -rf '+SB1_shift+'*') fixvis(vis=prefix+'_SB1_lines.ms', outputvis=SB1_shift, field=data_params['SB1']['field'], phasecenter='J2000 16h49m15.299722s -14d22m08.85131s') fixplanets(vis=SB1_shift, field=data_params['SB1']['field'], direction=common_dir) SB2_shift = prefix+'_SB2_lines_shift.ms' os.system('rm -rf '+SB2_shift+'*') fixvis(vis=prefix+'_SB2_lines.ms', outputvis=SB2_shift, field=data_params['SB2']['field'], phasecenter='J2000 16h49m15.301179s -14d22m08.95623s') fixplanets(vis=SB2_shift, field=data_params['SB2']['field'], direction=common_dir) SB3_shift = prefix+'_SB3_lines_shift.ms' os.system('rm -rf '+SB3_shift+'*') fixvis(vis=prefix+'_SB3_lines.ms', outputvis=SB3_shift, field=data_params['SB3']['field'], phasecenter='ICRS 16h49m15.292292s -14d22m09.02254s') fixplanets(vis=SB3_shift, field=data_params['SB3']['field'], direction=common_dir) SB4_shift = prefix+'_SB4_lines_shift.ms' os.system('rm -rf '+SB4_shift+'*') fixvis(vis=prefix+'_SB4_lines.ms', outputvis=SB4_shift, field=data_params['SB4']['field'], phasecenter='ICRS 16h49m15.293687s -14d22m09.08240s') fixplanets(vis=SB4_shift, field=data_params['SB4']['field'], direction=common_dir) SB5_shift = prefix+'_SB5_lines_shift.ms' os.system('rm -rf '+SB5_shift+'*') fixvis(vis=prefix+'_SB5_lines.ms', outputvis=SB5_shift, field=data_params['SB5']['field'], phasecenter='ICRS 16h49m15.293254s -14d22m09.06211s') fixplanets(vis=SB5_shift, field=data_params['SB5']['field'], direction=common_dir) LB1_shift = prefix+'_LB1_lines_shift_exec0.ms' os.system('rm -rf '+LB1_shift+'*') fixvis(vis=prefix+'_LB1_lines_exec0.ms', outputvis=LB1_shift, field=data_params['LB1']['field'], phasecenter='ICRS 16h49m15.293513s -14d22m09.04925s') fixplanets(vis=LB1_shift, field=data_params['LB1']['field'], direction=common_dir) LB1_shift = prefix+'_LB1_lines_shift_exec1.ms' os.system('rm -rf '+LB1_shift+'*') fixvis(vis=prefix+'_LB1_lines_exec1.ms', outputvis=LB1_shift, field=data_params['LB1']['field'], phasecenter='ICRS 16h49m15.293891s -14d22m09.05974s') fixplanets(vis=LB1_shift, field=data_params['LB1']['field'], direction=common_dir) rescale_flux(prefix+'_SB1_lines_shift.ms', [0.919]) rescale_flux(prefix+'_SB2_lines_shift.ms', [1.072]) rescale_flux(prefix+'_SB4_lines_shift.ms', [0.965]) """ Apply the self-calibration solutions to the SB data """ SB_CO = prefix+'_SB_CO' os.system('rm -rf %s*' % SB_CO) concat(vis=[prefix+'_SB1_lines_shift_rescaled.ms', prefix+'_SB2_lines_shift_rescaled.ms', prefix+'_SB3_lines_shift.ms', prefix+'_SB4_lines_shift_rescaled.ms', prefix+'_SB5_lines_shift.ms'], concatvis=SB_CO+'.ms', dirtol='0.1arcsec', copypointing=False) applycal(vis=SB_CO+'.ms', spw='0~16', gaintable=[prefix+'_SB.p1', prefix+'_SB.p2', prefix+'_SB.ap'], interp='linearPD', calwt=True, flagbackup=False) os.system('rm -rf '+SB_CO+'_selfcal.ms*') split(vis=SB_CO+'.ms', outputvis=SB_CO+'_selfcal.ms', datacolumn='corrected', keepflags=False) # save some space by deleting intermediate files os.system('rm -rf '+SB_CO+'.ms*') """ Now apply the self-calibration solutions to the combined dataset """ combined_CO = prefix+'_combined_CO' os.system('rm -rf '+combined_CO+'.ms*') concat(vis=[SB_CO+'_selfcal.ms', prefix+'_LB1_lines_exec0_shift.ms', prefix+'_LB1_lines_exec1_shift.ms'], concatvis=combined_CO+'.ms', dirtol='0.1arcsec', copypointing=False) spwmap = [0,0,0,3,3,3,3,3,3,9,9,11,11,13,13,13,13,17,17,17,17,21,21,21,21] applycal(vis=combined_CO+'.ms', spw='0~24', spwmap=[spwmap]*5, gaintable=[prefix+'_combined.p1', prefix+'_combined.p2', prefix+'_combined.p3', prefix+'_combined.p4', prefix+'_combined.ap'], interp='linearPD', calwt=True, applymode='calonly') os.system('rm -rf '+combined_CO+'_selfcal.ms*') split(vis=combined_CO+'.ms', outputvis=combined_CO+'_selfcal.ms', field='', spw='1,10,12,13,20,24', datacolumn='corrected', keepflags=False) # save some space by deleting intermediate files os.system('rm -rf '+combined_CO+'.ms*') os.system('rm -rf '+prefix+'_LB1_lines_exec*ms*') """ Continuum subtraction """ fitspw = '0:0~460;840~959, 1~2:0~135;325~479, 3~5:0~1873;1969~3839' os.system('rm -rf '+combined_CO+'_selfcal.ms.contsub*') uvcontsub(vis=combined_CO+'_selfcal.ms', spw='0~5', fitspw=fitspw, excludechans=False, solint='int', fitorder=1, want_cont=False) """ Define the channels of interest """ chanstart = '-4.0km/s' chanwidth = '0.35km/s' nchan = 60 """ Split and regrid into the channels of interest """ # Continuum-subtracted SB_only = prefix+'_CO_SBonly.ms.contsub' os.system('rm -rf '+SB_only+'*') split(vis=combined_CO+'_selfcal.ms.contsub', outputvis=SB_only, spw='0,1,2', datacolumn='data') SB_cvel = SB_only+'.cvel' os.system('rm -rf '+SB_cvel+'*') mstransform(vis=SB_only, outputvis=SB_cvel, keepflags=False, datacolumn='data', regridms=True, mode='velocity', start=chanstart, width=chanwidth, nchan=nchan, outframe='LSRK', veltype='radio', restfreq='230.538GHz') LB_only = prefix+'_CO_LBonly.ms.contsub' os.system('rm -rf '+LB_only+'*') split(vis=combined_CO+'_selfcal.ms.contsub', outputvis=LB_only, spw='3,4,5', datacolumn='data') LB_cvel = LB_only+'.cvel' os.system('rm -rf '+LB_cvel+'*') mstransform(vis=LB_only, outputvis=LB_cvel, keepflags=False, datacolumn='data', regridms=True, mode='velocity', start=chanstart, width=chanwidth, nchan=nchan, outframe='LSRK', veltype='radio', restfreq='230.538GHz') contsub_concat = prefix+'_COcube.ms.contsub' os.system('rm -rf '+contsub_concat) concat(vis=[SB_cvel, LB_cvel], concatvis=contsub_concat, dirtol='0.1arcsec', copypointing=False) # With continuum cSB_only = prefix+'_CO_SBonly.ms' os.system('rm -rf '+cSB_only) split(vis=combined_CO+'_selfcal.ms', outputvis=cSB_only, spw='0,1,2', datacolumn='data') cSB_cvel = cSB_only+'.cvel' os.system('rm -rf '+cSB_cvel) mstransform(vis=cSB_only, outputvis=cSB_cvel, keepflags=False, datacolumn='data', regridms=True, mode='velocity', start=chanstart, width=chanwidth, nchan=nchan, outframe='LSRK', veltype='radio', restfreq='230.538GHz') cLB_only = prefix+'_CO_LBonly.ms' os.system('rm -rf '+cLB_only) split(vis=combined_CO+'_selfcal.ms', outputvis=cLB_only, spw='3,4,5', datacolumn='data') cLB_cvel = cLB_only+'.cvel' os.system('rm -rf '+cLB_cvel) mstransform(vis=cLB_only, outputvis=cLB_cvel, keepflags=False, datacolumn='data', regridms=True, mode='velocity', start=chanstart, width=chanwidth, nchan=nchan, outframe='LSRK', veltype='radio', restfreq='230.538GHz') cont_concat = prefix+'_COcube.ms' os.system('rm -rf '+cont_concat) concat(vis=[cSB_cvel, cLB_cvel], concatvis=cont_concat, dirtol='0.1arcsec', copypointing=False) """ Imaging (continuum-subtracted only) """ imagename = prefix+'_CO' for ext in ['.image','.mask','.model','.pb','.psf','.residual','.sumwt']: os.system('rm -rf '+ imagename + ext) tclean(vis=prefix+'_COcube.ms.contsub', imagename=imagename, specmode='cube', imsize=1500, deconvolver='multiscale', start=chanstart, width=chanwidth, nchan=nchan, outframe='LSRK', veltype='radio', restfreq='230.538GHz', cell='0.01arcsec', scales = [0,10,25,75,150,250], gain=0.1, niter=50000, weighting='briggs', robust=1.0, threshold='3mJy', interactive=False, nterms=1, restoringbeam='common', uvtaper=['.025arcsec','.010arcsec','10deg']) """ Final outputs """ """ Save the final MS files """ os.system('tar cvzf '+prefix+'_CO.ms.tgz '+prefix+'_COcube.ms.contsub') os.system('tar cvzf '+prefix+'_COcont.ms.tgz '+prefix+'_COcube.ms') """ Save the imaged datacube """ exportfits(imagename+'.image', imagename+'.fits')