### ### Script for PdBI project L02D data: ### A Keplerian disk around the protostar IRAS 20126+4104 ### Single field, 1 & 3mm continuum DSB ### many lines in each sub-band with strong continuum emission. ### Day 2 - 26 Jan 2002 data ### ### Steven T. Myers (NRAO) Version 13-Jan-2004 ### ########################################################### ## Tested version stable_ss3 v1.9 Build 275 ## Expected run time 9 minutes ########################################################### ## Include relevant tools if not already done: ## include 'almati2ms.g'; include 'ms.g'; include 'msplot.g'; include 'flagger.g' include 'imager.g'; include 'calibrater.g'; include 'iramcalibrater.g'; include 'viewer.g'; include 'image.g'; include 'imagepol.g'; ########################################################### ## Fill Data: (unless already in ms) ## almatifiller(msfile='26-jan-2002-l02d.ms', fitsdir='/home/cluster4/jmcmulli/almatest/L102D/data', pattern='26-jan-2002-l02d-*', append=F, compress=F, obsmode='CORR', chanzero='TIME_AVG'); ########################################################### ## Obtain a detailed summary (with default catalog tool): ## dc.summary('26-jan-2002-l02d.ms'); ########################################################### ## Flag shadowed data: ## shadow(msname="26-jan-2002-l02d.ms", trial=F, minsep=15.0); ########################################################### ## Swap out IRAM PdBI on-line phase corrected data if not good: ## # Phase correction selection of all data: ical:=iramcalibrater(msname="26-jan-2002-l02d.ms", initcal=T); ical.phcor(trial=F); #Ignore error: #: index (= 47) out of range, array length = 45 # File: iramcalutil.g, Line 691 # Stack: .(), iramcalibrater.g line 200 # .() ical.done(); ########################################################### ## Manually plot the data for visual inspection, ## interactively flag bad data: ## ## see the data in x-y plot format: #plot := msplot(msfile='26-jan-2002-l02d.ms', edit=T); # Allow editing as well as inspection. See cookbook for editing steps. ########################################################### ## Flag bad and irrelevant data on first day ## - the real data for this file is on 28-Jan ## Flag spectral windows of same shape together ## - note flagger cannot handle different shaped windows ## - could use autoflag instead ## fg:=flagger('26-jan-2002-l02d.ms'); fg.settimerange(starttime='27-Jan-2002/05:00:00.0', endtime='27-Jan-2002/07:00:00.0'); # First the 1-channel windows: fg.setids(spectralwindowid=[1,2,5,6,9,10,13,14,17,18,21,22,25,26,29,30]); fg.flag() #Should see: Flagging 2400 rows # Next the 128-channel windows: fg.setids(spectralwindowid=[7,8,15,16,19,20,27,28,31,32]); fg.flag() #Should see: Flagging 1500 rows # Next the 512-channel windows: fg.setids(spectralwindowid=[3,4,11,12,23,24]); fg.flag() #Should see: Flagging 900 rows # Now the Gibbs channels for 128-ch 1.25MHz wins fg.reset(); fg.setids(spectralwindowid=[7,8,19,20]); fg.setchan(chan=[1,2,127,128]); # NOTE - ignore the errors from setchan #: different-length operands in expression (4 vs. 32): # (chan > private.nchan) # File: flagger.g, Line 365 # Stack: any(), flagger.g line 365 # .() fg.flag(); #Should see: Flagging 4440 rows # Now the Gibbs channels for 128-ch 2.5MHz wins fg.reset(); fg.setids(spectralwindowid=[15,16,27,28,31,32]); fg.setchan(chan=[1,2,64,65,127,128]); # NOTE - ignore the errors from setchan #: different-length operands in expression (6 vs. 32): # (chan > private.nchan) # File: flagger.g, Line 365 # Stack: any(), flagger.g line 365 # .() fg.flag(); #Should see: Flagging 6500 rows # Now the Gibbs channels for 512-ch wins fg.reset(); fg.setids(spectralwindowid=[3,4,11,12,23,24]); fg.setchan(chan=[1,2,3,4]); # NOTE - ignore the errors from setchan #: different-length operands in expression (4 vs. 32): # (chan > private.nchan) # File: flagger.g, Line 365 # Stack: any(), flagger.g line 365 # .() fg.flag(); #Should see: Flagging 6580 rows fg.done(); ########################################################### ## Set flux density scale for the flux density calibrators ## imgr:=imager(filename='26-jan-2002-l02d.ms'); # Set for the 3mm windows for(spw in 1:20){ imgr.setjy(fieldid=1,spwid=spw,fluxdensity=[10.87,0,0,0]); } for(spw in 1:20){ imgr.setjy(fieldid=3,spwid=spw,fluxdensity=[2.79,0,0,0]); } # Now set for the 1mm windows for(spw in 21:32){ imgr.setjy(fieldid=1,spwid=spw,fluxdensity=[6.00,0,0,0]); } for(spw in 21:32){ imgr.setjy(fieldid=2,spwid=spw,fluxdensity=[1.75,0,0,0]); } for(spw in 21:32){ imgr.setjy(fieldid=3,spwid=spw,fluxdensity=[1.80,0,0,0]); } ########################################################### ########## 3 mm Calibration ########## ########################################################### ########################################################### ## Get first cut phase solutions or BPass determination: ## cal:=calibrater(filename='26-jan-2002-l02d.ms'); ## Do 128-channel windows [7,8,15,16,19,20] ## Leave out first and last 2 unflagged channels ## First the LSB [7,15,19] cal.setdata(mode='channel', start=4, step=1, nchan=122, msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID IN [7,15,19]'); cal.setsolvegainspline(table='26-jan-2002-l02d.3mmLSB.gcal0', mode='PHAS', preavg=0, # no preavg too short, splinetime=10000, # spline fit timescale (~3hrs), refant=2, # refant = 2, npointaver=10, # minimize phase wrap pblms phasewrap=260, append=F); # create new table cal.solve(); shell('mv 26-jan-2002-l02d.ms.PHAS.ps 26-jan-2002-l02d.3mmLSB.gcal0.phase.ps') shell('mv 26-jan-2002-l02d.ms.phase.log 26-jan-2002-l02d.3mmLSB.gcal0.phase.log') ## Now the USB [8,16,20] cal.reset(); cal.setdata(mode='channel', start=4, step=1, nchan=122, msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID IN [8,16,20]'); cal.setsolvegainspline(table='26-jan-2002-l02d.3mmUSB.gcal0', mode='PHAS', preavg=0, # no preavg too short, splinetime=10000, # spline fit timescale (~3hrs), refant=2, # refant = 2, npointaver=10, # minimize phase wrap pblms phasewrap=260, append=F); # create new table cal.solve(); shell('mv 26-jan-2002-l02d.ms.PHAS.ps 26-jan-2002-l02d.3mmUSB.gcal0.phase.ps') shell('mv 26-jan-2002-l02d.ms.phase.log 26-jan-2002-l02d.3mmUSB.gcal0.phase.log') ########################################################### ## Derive BPass calibration: ## ## First for 3mm LSB continuum phase solns cal.reset(); cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID IN [7,15,19]'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mmUSB.gcal0'); # Note we are applying the USB phases to the LSB to preserve LSB/USB offsets cal.setsolvebandpoly(table='26-jan-2002-l02d.3mm.bpoly', degamp=12, degphase=24, visnorm=F, bpnorm=F, refant=2, append=F); cal.solve(); shell('mv 26-jan-2002-l02d.ms.RF_PHASE.ps 26-jan-2002-l02d.3mmLSB.BP_phase.ps') shell('mv 26-jan-2002-l02d.ms.RF_AMP.ps 26-jan-2002-l02d.3mmLSB.BP_amp.ps') ## Now for 3mm USB continuum phase solns... cal.reset(); cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID IN [8,16,20]'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mmUSB.gcal0'); ## NOTE append=T in caltable cal.setsolvebandpoly(table='26-jan-2002-l02d.3mm.bpoly', degamp=18, degphase=24, visnorm=F, bpnorm=F, refant=2, append=T); cal.solve(); shell('mv 26-jan-2002-l02d.ms.RF_PHASE.ps 26-jan-2002-l02d.3mmUSB.BP_phase.ps') shell('mv 26-jan-2002-l02d.ms.RF_AMP.ps 26-jan-2002-l02d.3mmUSB.BP_amp.ps') ########################################################### ## Phase calibrate 3mm LSB & USB 128-ch data: ## all cals ## cal.reset(); cal.setdata(mode='channel', start=4, step=1, nchan=122, msselect='FIELD_ID IN [1,2,3] && SPECTRAL_WINDOW_ID IN [7,8,15,16,19,20]'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mm.bpoly'); cal.setsolvegainspline(table='26-jan-2002-l02d.3mm.gcal', mode='PHAS', preavg=0, splinetime=10000, refant=2, npointaver=10, phasewrap=260, append=F); cal.solve(); shell('mv 26-jan-2002-l02d.ms.PHAS.ps 26-jan-2002-l02d.3mm.gcal.phase.ps') shell('mv 26-jan-2002-l02d.ms.phase.log 26-jan-2002-l02d.3mm.gcal.phase.log') ########################################################### ## Do a residual T solution to allow fluxscale ## cal.reset(); # Reset calibrator. cal.setdata(mode='channel', start=4, step=1, nchan=122, msselect='FIELD_ID IN [1,2,3] && SPECTRAL_WINDOW_ID IN [7,8,15,16,19,20]'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mm.bpoly'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal'); cal.setsolve(type='T', # Solve for any residual phase table='26-jan-2002-l02d.3mm.tcal', # errors with a timescale of t=300, # 5 min. refant=2); # Use ant 2 as reference. cal.solve(); # Solve for solutions. ########################################################### ## Get the fluxscale for MWC349 ## cal.fluxscale(tablein='26-jan-2002-l02d.3mm.tcal', tableout='26-jan-2002-l02d.3mm.fcal', reference="3C273 2013+370", transfer='MWC349'); ########################################################### ## Apply average fluxscale for MWC349 ## for(spw in 1:20){ imgr.setjy(fieldid=2,spwid=spw,fluxdensity=[1.06,0,0,0]); } ########################################################### ## Amplitude calibration of 3mm 128-ch data: ## all cals ## note that the phase solution in the table will be combined with the ## amplitude solution...doing effectively an incremental combination ## cal.reset(); cal.setdata(mode='channel', start=4, step=1, nchan=122, msselect='FIELD_ID IN [1,2,3] && SPECTRAL_WINDOW_ID IN [7,8,15,16,19,20]'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mm.bpoly'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal'); cal.setsolvegainspline(table='26-jan-2002-l02d.3mm.gcal', mode='AMP', preavg=0.0, splinetime=10000, refant=2, npointaver=10, phasewrap=260); cal.solve(); shell('mv 26-jan-2002-l02d.ms.AMP.ps 26-jan-2002-l02d.3mm.gcal.amp.ps') shell('mv 26-jan-2002-l02d.ms.amp.log 26-jan-2002-l02d.3mm.gcal.amp.log') ########################################################### ## Correct all of the 3mm 128ch data ## cal.reset(); cal.setdata(msselect='FIELD_ID IN [1,2,3,4] && SPECTRAL_WINDOW_ID IN [7,8,15,16,19,20]'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mm.bpoly'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal'); cal.correct(); ########################################################### ## Do separate bandpasses for the 512ch windows ## cal.reset(); cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID==3'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mm.bpoly'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal'); cal.setsolvebandpoly(table='26-jan-2002-l02d.3mmLSBa.bpoly', degamp=3, degphase=3, visnorm=F, bpnorm=F, refant=2, append=F); cal.solve(); shell('mv 26-jan-2002-l02d.ms.RF_PHASE.ps 26-jan-2002-l02d.3mmLSBa.BP_phase.ps') shell('mv 26-jan-2002-l02d.ms.RF_AMP.ps 26-jan-2002-l02d.3mmLSBa.BP_amp.ps') cal.reset(); cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID==11'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mm.bpoly'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal'); cal.setsolvebandpoly(table='26-jan-2002-l02d.3mmLSBb.bpoly', degamp=3, degphase=3, visnorm=F, bpnorm=F, refant=2, append=F); cal.solve(); shell('mv 26-jan-2002-l02d.ms.RF_PHASE.ps 26-jan-2002-l02d.3mmLSBb.BP_phase.ps') shell('mv 26-jan-2002-l02d.ms.RF_AMP.ps 26-jan-2002-l02d.3mmLSBb.BP_amp.ps') cal.reset(); cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID==4'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mm.bpoly'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal'); cal.setsolvebandpoly(table='26-jan-2002-l02d.3mmUSBa.bpoly', degamp=3, degphase=3, visnorm=F, bpnorm=F, refant=2, append=F); cal.solve(); shell('mv 26-jan-2002-l02d.ms.RF_PHASE.ps 26-jan-2002-l02d.3mmUSBa.BP_phase.ps') shell('mv 26-jan-2002-l02d.ms.RF_AMP.ps 26-jan-2002-l02d.3mmUSBa.BP_amp.ps') cal.reset(); cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID==12'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mm.bpoly'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal'); cal.setsolvebandpoly(table='26-jan-2002-l02d.3mmUSBb.bpoly', degamp=3, degphase=3, visnorm=F, bpnorm=F, refant=2, append=F); cal.solve(); shell('mv 26-jan-2002-l02d.ms.RF_PHASE.ps 26-jan-2002-l02d.3mmUSBb.BP_phase.ps') shell('mv 26-jan-2002-l02d.ms.RF_AMP.ps 26-jan-2002-l02d.3mmUSBb.BP_amp.ps') ########################################################### ## Correct the 3mm 512ch data ## cal.reset(); cal.setdata(msselect='FIELD_ID IN [1,2,3,4] && SPECTRAL_WINDOW_ID IN [3]'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mmLSBa.bpoly'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal'); cal.correct(); cal.reset(); cal.setdata(msselect='FIELD_ID IN [1,2,3,4] && SPECTRAL_WINDOW_ID IN [11]'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mmLSBb.bpoly'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal'); cal.correct(); cal.reset(); cal.setdata(msselect='FIELD_ID IN [1,2,3,4] && SPECTRAL_WINDOW_ID IN [4]'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mmUSBa.bpoly'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal'); cal.correct(); cal.reset(); cal.setdata(msselect='FIELD_ID IN [1,2,3,4] && SPECTRAL_WINDOW_ID IN [12]'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mmUSBb.bpoly'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal'); cal.correct(); ########################################################### ########## 1 mm Calibration ########## ########################################################### ########################################################### ## Initial phase solutions: ## Use the 128ch windows [27,28,31,32] ## First the LSB channels [27,31] ## cal.reset(); cal.setdata(mode='channel', start=4, step=1, nchan=122, msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID IN [27,31]'); cal.setsolvegainspline(table='26-jan-2002-l02d.1mmLSB.gcal0', mode='PHAS', preavg=0, splinetime=10000, refant=2, npointaver=10, phasewrap=260, append=F); cal.solve(); shell('mv 26-jan-2002-l02d.ms.PHAS.ps 26-jan-2002-l02d.1mmLSB.gcal0.phase.ps') shell('mv 26-jan-2002-l02d.ms.phase.log 26-jan-2002-l02d.1mmLSB.gcal0.phase.log') ## Now the USB channels [28,32] cal.reset(); cal.setdata(mode='channel', start=4, step=1, nchan=122, msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID IN [28,32]'); cal.setsolvegainspline(table='26-jan-2002-l02d.1mmUSB.gcal0', mode='PHAS', preavg=0, splinetime=10000, refant=2, npointaver=10, phasewrap=260, append=F); cal.solve(); shell('mv 26-jan-2002-l02d.ms.PHAS.ps 26-jan-2002-l02d.1mmUSB.gcal0.phase.ps') shell('mv 26-jan-2002-l02d.ms.phase.log 26-jan-2002-l02d.1mmUSB.gcal0.phase.log') ########################################################### ## Bandpass calibration: ## First for 1mm LSB ## cal.reset(); cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID IN [27,31]'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.1mmUSB.gcal0'); cal.setsolvebandpoly(table='26-jan-2002-l02d.1mm.bpoly', degamp=6, degphase=12, visnorm=F, bpnorm=F, refant=2, append=F); cal.solve(); shell('mv 26-jan-2002-l02d.ms.RF_PHASE.ps 26-jan-2002-l02d.1mmLSB.BP_phase.ps') shell('mv 26-jan-2002-l02d.ms.RF_AMP.ps 26-jan-2002-l02d.1mmLSB.BP_amp.ps') ## Now for 1mm USB cal.reset(); cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID IN [28,32]'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.1mmUSB.gcal0'); cal.setsolvebandpoly(table='26-jan-2002-l02d.1mm.bpoly', degamp=6, degphase=12, visnorm=F, bpnorm=F, refant=2, append=T); cal.solve(); shell('mv 26-jan-2002-l02d.ms.RF_PHASE.ps 26-jan-2002-l02d.1mmUSB.BP_phase.ps') shell('mv 26-jan-2002-l02d.ms.RF_AMP.ps 26-jan-2002-l02d.1mmUSB.BP_amp.ps') ########################################################### ## Phase calibrate 1mm LSB & USB: ## cal.reset(); cal.setdata(mode='channel', start=4, step=1, nchan=122, msselect='FIELD_ID IN [1,2,3] && SPECTRAL_WINDOW_ID IN [27,28,31,32]'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.1mm.bpoly'); ## ## Doing a curve transfer of 3mm solution to 1mm to reduce phase wrapping ## incremental solution with cumulative output table ## cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal'); cal.setsolvegainspline(table='26-jan-2002-l02d.1mm.gcal', mode='PHAS', preavg=0, splinetime=10000, refant=2, npointaver=10, phasewrap=260, append=F); cal.solve(); shell('mv 26-jan-2002-l02d.ms.PHAS.ps 26-jan-2002-l02d.1mm.gcal.phase.ps') shell('mv 26-jan-2002-l02d.ms.phase.log 26-jan-2002-l02d.1mm.gcal.phase.log') ########################################################### ## Amplitude calibration of 1mm : ## Note, we had used setjy for all 1mm calibrators ## so no need for fluxscale here ## cal.reset(); cal.setdata(mode='channel', start=4, step=1, nchan=122, msselect='FIELD_ID IN [1,2,3] && SPECTRAL_WINDOW_ID IN [27,28,31,32]'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.1mm.bpoly'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.1mm.gcal'); cal.setsolvegainspline(table='26-jan-2002-l02d.1mm.gcal', mode='AMP', preavg=0.0, splinetime=10000, refant=2, npointaver=10, phasewrap=260); cal.solve(); shell('mv 26-jan-2002-l02d.ms.AMP.ps 26-jan-2002-l02d.1mm.gcal.amp.ps') shell('mv 26-jan-2002-l02d.ms.amp.log 26-jan-2002-l02d.1mm.gcal.amp.log') ########################################################### ## Correct all of the 1mm 128ch data ## cal.reset(); cal.setdata(msselect='FIELD_ID IN [1,2,3,4] && SPECTRAL_WINDOW_ID IN [27,28,31,32]'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.1mm.bpoly'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.1mm.gcal'); cal.correct(); ########################################################### ## Do separate bandpasses for the 512ch windows ## cal.reset(); cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID==23'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.1mm.bpoly'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.1mm.gcal'); cal.setsolvebandpoly(table='26-jan-2002-l02d.1mmLSBa.bpoly', degamp=3, degphase=3, visnorm=F, bpnorm=F, refant=2, append=F); cal.solve(); shell('mv 26-jan-2002-l02d.ms.RF_PHASE.ps 26-jan-2002-l02d.1mmLSBa.BP_phase.ps') shell('mv 26-jan-2002-l02d.ms.RF_AMP.ps 26-jan-2002-l02d.1mmLSBa.BP_amp.ps') cal.reset(); cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID==24'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.1mm.bpoly'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.1mm.gcal'); cal.setsolvebandpoly(table='26-jan-2002-l02d.1mmUSBa.bpoly', degamp=3, degphase=3, visnorm=F, bpnorm=F, refant=2, append=F); cal.solve(); shell('mv 26-jan-2002-l02d.ms.RF_PHASE.ps 26-jan-2002-l02d.1mmUSBa.BP_phase.ps') shell('mv 26-jan-2002-l02d.ms.RF_AMP.ps 26-jan-2002-l02d.1mmUSBa.BP_amp.ps') ########################################################### ## Correct the 1mm 512ch data ## cal.reset(); cal.setdata(msselect='FIELD_ID IN [1,2,3,4] && SPECTRAL_WINDOW_ID IN [23]'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.1mmLSBa.bpoly'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.1mm.gcal'); cal.correct(); cal.reset(); cal.setdata(msselect='FIELD_ID IN [1,2,3,4] && SPECTRAL_WINDOW_ID IN [24]'); cal.setapply(type='BPOLY', table='26-jan-2002-l02d.1mmUSBa.bpoly'); cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.1mm.gcal'); cal.correct(); ########################################################### ## Done with calibration ## cal.done(); imgr.done(); ########################################################### ## Split out continuum, 128ch line, and 512ch line data ## myms:=ms('26-jan-2002-l02d.ms', readonly=F); ## 3mm ## averaging the 122 middle channel of the 128ch continuum spw to 1 channel myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmCONT.ms', fieldids=[4], spwids=[7,8,15,16,19,20], nchan=[1], start=[4], step=[122], whichcol='CORRECTED_DATA'); ## Now the 6 128ch bands... myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmLSB128a.ms', fieldids=[4], spwids=[7], nchan=[128], start=[1], step=[1], whichcol='CORRECTED_DATA'); myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmLSB128b.ms', fieldids=[4], spwids=[15], nchan=[128], start=[1], step=[1], whichcol='CORRECTED_DATA'); myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmLSB128c.ms', fieldids=[4], spwids=[19], nchan=[128], start=[1], step=[1], whichcol='CORRECTED_DATA'); myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmUSB128a.ms', fieldids=[4], spwids=[8], nchan=[128], start=[1], step=[1], whichcol='CORRECTED_DATA'); myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmUSB128b.ms', fieldids=[4], spwids=[16], nchan=[128], start=[1], step=[1], whichcol='CORRECTED_DATA'); myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmUSB128c.ms', fieldids=[4], spwids=[20], nchan=[128], start=[1], step=[1], whichcol='CORRECTED_DATA'); ## Now the 4 512ch line bands... myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmLSB512a.ms', fieldids=[4], spwids=[3], nchan=[512], start=[1], step=[1], whichcol='CORRECTED_DATA'); myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmLSB512b.ms', fieldids=[4], spwids=[11], nchan=[512], start=[1], step=[1], whichcol='CORRECTED_DATA'); myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmUSB512a.ms', fieldids=[4], spwids=[4], nchan=[512], start=[1], step=[1], whichcol='CORRECTED_DATA'); myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmUSB512b.ms', fieldids=[4], spwids=[12], nchan=[512], start=[1], step=[1], whichcol='CORRECTED_DATA'); ## 1mm ## averaging the 122 middle channel of the 128ch continuum spw to 1 channel myms.split(outputms='20126+4104.26-jan-2002-l02d.1mmCONT.ms', fieldids=[4], spwids=[27,28,31,32], nchan=[1], start=[4], step=[122], whichcol='CORRECTED_DATA'); ## Now the 4 128ch line bands... myms.split(outputms='20126+4104.26-jan-2002-l02d.1mmLSB128a.ms', fieldids=[4], spwids=[27], nchan=[128], start=[1], step=[1], whichcol='CORRECTED_DATA'); myms.split(outputms='20126+4104.26-jan-2002-l02d.1mmLSB128b.ms', fieldids=[4], spwids=[31], nchan=[128], start=[1], step=[1], whichcol='CORRECTED_DATA'); myms.split(outputms='20126+4104.26-jan-2002-l02d.1mmUSB128a.ms', fieldids=[4], spwids=[28], nchan=[128], start=[1], step=[1], whichcol='CORRECTED_DATA'); myms.split(outputms='20126+4104.26-jan-2002-l02d.1mmUSB128b.ms', fieldids=[4], spwids=[32], nchan=[128], start=[1], step=[1], whichcol='CORRECTED_DATA'); ## Now the 2 512ch line bands... myms.split(outputms='20126+4104.26-jan-2002-l02d.1mmLSB512a.ms', fieldids=[4], spwids=[23], nchan=[512], start=[1], step=[1], whichcol='CORRECTED_DATA'); myms.split(outputms='20126+4104.26-jan-2002-l02d.1mmUSB512a.ms', fieldids=[4], spwids=[24], nchan=[512], start=[1], step=[1], whichcol='CORRECTED_DATA'); myms.done(); ########################################################### ########################################################### ########################################################### ########################################################### ########## 3 mm Continuum Imaging ########## ########################################################### imgr:=imager(filename='20126+4104.26-jan-2002-l02d.3mmCONT.ms'); imgr.setdata(mode='channel', nchan=1, start=1, step=1, spwid=1, fieldid=1); ## ## IRAM 15m FWHM should be about 50" at 96GHz ## imgr.setimage(nx=256, ny=256, # Set image plane parameters cellx='0.25arcsec', celly='0.25arcsec', stokes='I', # Image total intensity mode='mfs', # Doesnt matter here nchan=1, start=1, # step=1, spwid=1, fieldid=1); imgr.weight(type='briggs', # Set up Robust weighting rmode='norm', robust=0.5); ## ## Make dirty map ## imgr.makeimage(type='corrected',image='20126+4104.3mmCONT.dirty'); ## ## Find peak and box around it ## im:=image('20126+4104.3mmCONT.dirty'); im.statistics(statsout=statrec); im.done(); print 'Processed 3mm continuum dirty image'; print 'Pixel position of maximum should be [129 128 1 1]:'; print statrec.maxpos; # Should be [129 128 1 1] print 'Maximum dirty flux should be 0.02082Jy:'; print statrec.max; # Should be 0.0208230969 maxblc:= [statrec.maxpos[1]-10, statrec.maxpos[2]-10, 1, 1]; maxtrc:= [statrec.maxpos[1]+10, statrec.maxpos[2]+10, 1, 1]; imgr.boxmask(mask='maxbox', blc=maxblc, trc=maxtrc); ## ## Clean 50 iterations around peak ## imgr.clean(algorithm='clark', # use the Clark CLEAN algorithm niter=50, # 50 iterations gain=0.1, model='20126+4104.3mmCONT.model', image='20126+4104.3mmCONT.clean', residual='20126+4104.3mmCONT.resid', mask='maxbox', # Use our box interactive=F); ## ## Clean 50 more iterations in inner quarter ## shell('cp -r 20126+4104.3mmCONT.model 20126+4104.3mmCONT.model2') imgr.clean(algorithm='clark', niter=50, # 50 more iterations gain=0.1, model='20126+4104.3mmCONT.model2', image='20126+4104.3mmCONT.clean2', residual='20126+4104.3mmCONT.resid2', mask='', # Inner quarter (no mask set) interactive=F); imgr.done(); ## ## Check results ## im:=image('20126+4104.3mmCONT.clean2'); im.statistics(statsout=statrec); im.done(); print 'Processed 3mm continuum clean image'; print 'Pixel position of maximum should be [129 128 1 1]:'; print statrec.maxpos; # Should be [129 128 1 1] print 'Maximum clean flux should be 0.01807Jy:'; print statrec.max; # Should be 0.0180715639 ########################################################### ########## 1 mm Continuum Imaging ########## ########################################################### imgr:=imager(filename='20126+4104.26-jan-2002-l02d.1mmCONT.ms'); imgr.setdata(mode='channel', nchan=1, start=1, step=1, spwid=1, fieldid=1); imgr.setimage(nx=256, ny=256, # Set image plane parameters cellx='0.1arcsec', celly='0.1arcsec', stokes='I', # Image total intensity mode='mfs', # Doesnt matter here nchan=1, start=1, # step=1, spwid=1, fieldid=1); imgr.weight(type='briggs', # Set up Robust weighting rmode='norm', robust=0.5); ## ## Dirty image, find peak and box for clean ## imgr.makeimage(type='corrected',image='20126+4104.1mmCONT.dirty'); im:=image('20126+4104.1mmCONT.dirty'); im.statistics(statsout=statrec); im.done(); print 'Processed 1mm continuum dirty image'; print 'Pixel position of maximum should be [129 129 1 1]:'; print statrec.maxpos; # Should be [129 129 1 1] print 'Maximum dirty flux should be 0.1386Jy:'; print statrec.max; # Should be 0.138645902 maxblc:= [statrec.maxpos[1]-10, statrec.maxpos[2]-10, 1, 1]; maxtrc:= [statrec.maxpos[1]+10, statrec.maxpos[2]+10, 1, 1]; imgr.boxmask(mask='maxbox', blc=maxblc, trc=maxtrc); ## ## Clean in box ## imgr.clean(algorithm='clark', # use the Clark CLEAN algorithm niter=50, # 50 iterations gain=0.1, model='20126+4104.1mmCONT.model', image='20126+4104.1mmCONT.clean', residual='20126+4104.1mmCONT.resid', mask='maxbox', # Use our box interactive=F); ## ## Further clean in inner quarter ## shell('cp -r 20126+4104.1mmCONT.model 20126+4104.1mmCONT.model2') imgr.clean(algorithm='clark', niter=50, # 50 more iterations gain=0.1, model='20126+4104.1mmCONT.model2', image='20126+4104.1mmCONT.clean2', residual='20126+4104.1mmCONT.resid2', mask='', # Inner quarter (no mask set) interactive=F); imgr.done(); ## ## Check results ## im:=image('20126+4104.1mmCONT.clean2'); im.statistics(statsout=statrec); im.done(); print 'Processed 1mm continuum clean image'; print 'Pixel position of maximum should be [129 129 1 1]:'; print statrec.maxpos; # Should be [129 129 1 1] print 'Maximum clean flux should be 0.1280Jy:'; print statrec.max; # Should be 0.128039777 ########################################################### ########## 1 mm Line cubes (512ch) ########## ########################################################### ## ## LSB ## imgr:=imager(filename='20126+4104.26-jan-2002-l02d.1mmLSB512a.ms'); imgr.setdata(mode='channel', nchan=512, start=1, step=1, spwid=1, fieldid=1); imgr.setimage(nx=128, ny=128, # Set image plane parameters cellx='0.15arcsec', celly='0.15arcsec', stokes='I', # Image total intensity mode='channel', # Cube nchan=256, # Average together 2 chans start=1, step=2, spwid=1, fieldid=1); imgr.weight(type='briggs', # Set up Robust weighting rmode='norm', robust=0.5); imgr.makeimage(type='corrected',image='20126+4104.1mmLSB512a.dirty'); imgr.done(); ## ## USB ## imgr:=imager(filename='20126+4104.26-jan-2002-l02d.1mmUSB512a.ms'); imgr.setdata(mode='channel', nchan=512, start=1, step=1, spwid=1, fieldid=1); imgr.setimage(nx=128, ny=128, # Set image plane parameters cellx='0.15arcsec', celly='0.15arcsec', stokes='I', # Image total intensity mode='channel', # Cube nchan=256, # Average together 2 chans start=1, step=2, spwid=1, fieldid=1); imgr.weight(type='briggs', # Set up Robust weighting rmode='norm', robust=0.5); imgr.makeimage(type='corrected',image='20126+4104.1mmUSB512a.dirty'); imgr.done(); ########################################################### ########## 3 mm Line cubes (128ch x 1.25MHz) ########## ########################################################### ## ## LSB ## imgr:=imager(filename='20126+4104.26-jan-2002-l02d.3mmLSB128a.ms'); imgr.setdata(mode='channel', nchan=128, start=1, step=1, spwid=1, fieldid=1); imgr.setimage(nx=256, ny=256, # Set image plane parameters cellx='0.25arcsec', celly='0.25arcsec', stokes='I', # Image total intensity mode='channel', # Cube nchan=128, # each channel start=1, step=1, spwid=1, fieldid=1); imgr.weight(type='briggs', # Set up Robust weighting rmode='norm', robust=0.5); imgr.makeimage(type='corrected',image='20126+4104.3mmLSB128a.dirty'); imgr.done(); imgr:=imager(filename='20126+4104.26-jan-2002-l02d.3mmLSB128c.ms'); imgr.setdata(mode='channel', nchan=128, start=1, step=1, spwid=1, fieldid=1); imgr.setimage(nx=256, ny=256, # Set image plane parameters cellx='0.25arcsec', celly='0.25arcsec', stokes='I', # Image total intensity mode='channel', # Cube nchan=128, # each channel start=1, step=1, spwid=1, fieldid=1); imgr.weight(type='briggs', # Set up Robust weighting rmode='norm', robust=0.5); imgr.makeimage(type='corrected',image='20126+4104.3mmLSB128c.dirty'); imgr.done(); ## ## USB ## imgr:=imager(filename='20126+4104.26-jan-2002-l02d.3mmUSB128a.ms'); imgr.setdata(mode='channel', nchan=128, start=1, step=1, spwid=1, fieldid=1); imgr.setimage(nx=256, ny=256, # Set image plane parameters cellx='0.25arcsec', celly='0.25arcsec', stokes='I', # Image total intensity mode='channel', # Cube nchan=128, # each channel start=1, step=1, spwid=1, fieldid=1); imgr.weight(type='briggs', # Set up Robust weighting rmode='norm', robust=0.5); imgr.makeimage(type='corrected',image='20126+4104.3mmUSB128a.dirty'); imgr.done(); imgr:=imager(filename='20126+4104.26-jan-2002-l02d.3mmUSB128c.ms'); imgr.setdata(mode='channel', nchan=128, start=1, step=1, spwid=1, fieldid=1); imgr.setimage(nx=256, ny=256, # Set image plane parameters cellx='0.25arcsec', celly='0.25arcsec', stokes='I', # Image total intensity mode='channel', # Cube nchan=128, # each channel start=1, step=1, spwid=1, fieldid=1); imgr.weight(type='briggs', # Set up Robust weighting rmode='norm', robust=0.5); imgr.makeimage(type='corrected',image='20126+4104.3mmUSB128c.dirty'); imgr.done(); ########################################################### ########## 3 mm Line cubes (512ch) ########## ########################################################### ## ## LSB ## imgr:=imager(filename='20126+4104.26-jan-2002-l02d.3mmLSB512a.ms'); imgr.setdata(mode='channel', nchan=512, start=1, step=1, spwid=1, fieldid=1); imgr.setimage(nx=128, ny=128, # Set image plane parameters cellx='0.4arcsec', celly='0.4arcsec', stokes='I', # Image total intensity mode='channel', # Cube nchan=256, # Average together 2 chans start=1, step=2, spwid=1, fieldid=1); imgr.weight(type='briggs', # Set up Robust weighting rmode='norm', robust=0.5); imgr.makeimage(type='corrected',image='20126+4104.3mmLSB512a.dirty'); imgr.done(); imgr:=imager(filename='20126+4104.26-jan-2002-l02d.3mmLSB512b.ms'); imgr.setdata(mode='channel', nchan=512, start=1, step=1, spwid=1, fieldid=1); imgr.setimage(nx=128, ny=128, # Set image plane parameters cellx='0.4arcsec', celly='0.4arcsec', stokes='I', # Image total intensity mode='channel', # Cube nchan=256, # Average together 2 chans start=1, step=2, spwid=1, fieldid=1); imgr.weight(type='briggs', # Set up Robust weighting rmode='norm', robust=0.5); imgr.makeimage(type='corrected',image='20126+4104.3mmLSB512b.dirty'); imgr.done(); ## ## USB ## imgr:=imager(filename='20126+4104.26-jan-2002-l02d.3mmUSB512a.ms'); imgr.setdata(mode='channel', nchan=512, start=1, step=1, spwid=1, fieldid=1); imgr.setimage(nx=128, ny=128, # Set image plane parameters cellx='0.4arcsec', celly='0.4arcsec', stokes='I', # Image total intensity mode='channel', # Cube nchan=256, # Average together 2 chans start=1, step=2, spwid=1, fieldid=1); imgr.weight(type='briggs', # Set up Robust weighting rmode='norm', robust=0.5); imgr.makeimage(type='corrected',image='20126+4104.3mmUSB512a.dirty'); imgr.done(); imgr:=imager(filename='20126+4104.26-jan-2002-l02d.3mmUSB512b.ms'); imgr.setdata(mode='channel', nchan=512, start=1, step=1, spwid=1, fieldid=1); imgr.setimage(nx=128, ny=128, # Set image plane parameters cellx='0.4arcsec', celly='0.4arcsec', stokes='I', # Image total intensity mode='channel', # Cube nchan=256, # Average together 2 chans start=1, step=2, spwid=1, fieldid=1); imgr.weight(type='briggs', # Set up Robust weighting rmode='norm', robust=0.5); imgr.makeimage(type='corrected',image='20126+4104.3mmUSB512b.dirty'); imgr.done(); ## ## These cubes can be examined in the viewer #dv.gui(); ## ## I did not see any sign of line emission in this single day of data ## ## ## Done ##