########################################################################### # C. Brogan script for ALMA1.1 test: N1333 ########################################################################### ## Include relevant tools if not already done: ## include 'vlafiller.g'; include 'ms.g'; include 'flagger.g'; include 'msplot.g'; include 'calibrater.g'; include 'viewer.g'; include 'imager.g'; include 'image.g'; include 'imagepol.g'; system.output.olog := 'ALMATST1_n1333.log' ########################################################### ## Fill Data: ## mset:=fitstoms(fitsfile='N1333_1.UVFITS',msfile='N1333_1.uv.ms'); # Obtain summary listing mset.summary(verbose=T) mset.done() ######################################################## ##### Begin data editing ## First Flag Autocorrelations af:=autoflag('N1333_1.uv.ms'); af.setdata(); af.setselect(autocorr=T); # Select autocorrelations af.run(); af.done(); ## Interactive flagging plot:=msplot(msfile='N1333_1.uv.ms', edit=T) ## Flag bad correlations af:=autoflag('N1333_1.uv.ms') af.setdata() af.setselect(timerng="2003/05/02/21:38:30 2003/05/02/22:13:30") af.run() af.done() ## Flag edge channels af:=autoflag('N1333_1.uv.ms') af.setdata() af.setselect(chan=[1,3,61,63]) af.run() af.done() ## Flag antenna 9. It has low amplitudes. ## Note this is ONLY apparent if you look at an average of the ## channels af:=autoflag('N1333_1.uv.ms') af.setdata() af.setselect(ant=9) af.run() af.done() ## Flag antenna 14 for a certain timerange for low amplitudes ## Note this is ONLY apparent if you look at an average of the ## channels af:=autoflag('N1333_1.uv.ms') af.setdata() af.setselect(timerng="2003/05/02/18:45:00 2003/05/02/19:15:00",ant=14) af.run() af.done() ## Bad data on 0319+415 ## spid 1 all antennas af:=autoflag('N1333_1.uv.ms') af.setdata() af.setselect(timerng="2003/05/02/20:12:00 2003/05/02/20:20:00") af.run() af.done() ## Bad data on 0319+415 ## spid 2 antennas 8 and 22 af:=autoflag('N1333_1.uv.ms') af.setdata() af.setselect(timerng="2003/05/02/20:20:00 2003/05/02/20:31:00",ant=22) af.run() af.done() af:=autoflag('N1333_1.uv.ms') af.setdata() af.setselect(timerng="2003/05/02/20:20:00 2003/05/02/20:31:00",ant=8) af.run() af.done() ########################################################### ## Set the flux density of 3C~147 (fields 13&14): ## imgr:=imager(filename='N1333_1.uv.ms'); # Start the imager tool imgr.setjy(fieldid=13); # Calculate & set flux density for imgr.setjy(fieldid=14); # 0542+498 (3c147) ########################################################### ## Derive gain calibration solutions: ## The details of the spectral windows can be worked out during the ## apply phase, unnecessary for solving. It is necessary to split ## up solving for gain/bp and flux cal because flux cal is slightly ## resolved cal:=calibrater(filename='N1333_1.uv.ms'); # Create calibrater tool cal.reset(); # Reset apply/solve state # of the calibrater tool cal.setdata(msselect='FIELD_ID IN [1,7,15,16]') # Note already flagged end channels # Do flux cal. separately cal.setsolve(type='G', # Solve for G for each scan t=0, # thus t=0 refant=26, # Choose reference antenna 26: table='N1333_1.uv.gcal'); # a well-behaved antenna # near the center of the array. cal.state(); # Review setsolve settings cal.solve(); # Solve for the net complex gains # and write solutions to the table cal.reset(); cal.setdata(msselect='FIELD_ID IN [13,14]', uvrange=[0,50]); # Flux cal is resolved cal.setsolve(type='G', # Solve for G for each scan t=0, # thus t=0 refant=26, # Choose reference antenna 26: table='N1333_1.uv.gcal', # Need append=T to add to table append=T); # instead of overwritting cal.state(); # Review setsolve settings cal.solve(); # Solve for the net complex gains # and write solutions to the table cal.plotcal(tablename='N1333_1.uv.gcal'); # Inspect solutions ########################################################### ## Derive bandpass calibration solutions: ## cal.reset(); # Reset apply/solve state # cal.setdata(msselect='FIELD_ID IN [15,16]'); # Select bandpass calibrator # (0319+415 = 3C84) cal.setapply(type='G', # Apply G solutions before determing t=0.0, # bpass table='N1333_1.uv.gcal'); cal.setsolve(type='B', # Arrange to solve for a single t=100000, # bandpass solution for whole refant=26, # observation. i.e. very large time table='N1333_1.uv.bcal'); # cal.state(); # Review setapply/setsolve settings cal.solve(); # Solve for the bandpass solutions # and write them to the table # N1333_1.uv.bcal located on disk. cal.plotcal(tablename='N1333_1.uv.bcal', # Inspect the solutions plottype='AMP', multiplot=T); ########################################################### ## Transfer the flux density scale: ## cal.fluxscale(tablein='N1333_1.uv.gcal', # Transfer the flux density scale tableout='N1333_1.uv.fluxcal', # from 3c147 flux cal reference='0542+498', # to 0336+323 and 0319+415 transfer=['0336+323','0319+415']); cal.plotcal(tablename='N1333_1.uv.fluxcal', # Inspect the solutions plottype='AMP', multiplot=T); ########################################################### ## Correct the target source data: ## cal.reset(); # Reset setapply/setsolve cal.setdata(msselect='FIELD_ID IN [2,3,4,5,6]'); # Select spwd=1 target sources # 2,3,4,5,6 cal.setapply(type='G', # Apply flux-scaled G for t=0.0, # spwd=1 gain cal table='N1333_1.uv.fluxcal', # which is fieldid=1 select='FIELD_ID IN [1]'); cal.setapply(type='B', # Arrange to apply B solutions t=0.0, # from the N1333_1.uv.bcal table. table='N1333_1.uv.bcal', select='FIELD_ID IN [15]'); # Select spwd=1 bp cal cal.state(); # Review setapply settings cal.correct(); # Apply solutions and write the # CORRECTED_DATA column in the MS. # Repeat for spwd=2 cal.reset(); # Reset setapply/setsolve cal.setdata(msselect='FIELD_ID IN [8,9,10,11,12]'); # Select spwd=2 target sources # cal.setapply(type='G', # Apply flux-scaled G for t=0.0, # spwd=1 gain cal table='N1333_1.uv.fluxcal', # which is fieldid=7 select='FIELD_ID IN [7]'); cal.setapply(type='B', # Arrange to apply B solutions t=0.0, # from the N1333_1.uv.bcal table. table='N1333_1.uv.bcal', select='FIELD_ID IN [16]'); # Select spwd=1 bp cal cal.state(); # Review setapply settings cal.correct(); # Apply solutions and write the # CORRECTED_DATA column in the MS. cal.done(); ########################################################### ## Split out spwd=1 and 2 calibrated target source data: ## mset:= ms(filename="N1333_1.uv.ms", readonly=F); # Create the ms tool mset.split(outputms="N1333_1.split.ms", # Select output file name, fieldids=[2,3,4,5,6,8,9,10,11,12], # write out target fields whichcol="CORRECTED_DATA", # write out CORRECTED_DATA spwids=[1,2], nchan=56, #Leave out flagged end channels start=4, step=1); ## Get a summary mset.summary(); # Get a new summary so you know # the new field_id & spwid mset.done(); ############################################################# ## Plot target data plot:=msplot(msfile='N1333_1.split.ms', edit=T) ## Antenna 18 badly behaved at beginning af:=autoflag('N1333_1.split.ms') af.setdata() af.setselect(timerng="2003/05/02/18:30:00 2003/05/02/19:30:00",ant=18) af.run() af.done() ## Antenna 22 af:=autoflag('N1333_1.split.ms') af.setdata() af.setselect(timerng="2003/05/02/18:45:00 2003/05/02/18:51:00",ant=22) af.run() af.done() ## Antenna 15 af:=autoflag('N1333_1.split.ms') af.setdata() af.setselect(timerng="2003/05/02/22:38:00 2003/05/02/22:47:00",ant=15) af.run() af.done() ## Antenna 1 af:=autoflag('N1333_1.split.ms') af.setdata() af.setselect(timerng="2003/05/02/18:51:00 2003/05/02/19:08:30",ant=1) af.run() af.done() ## Antenna 20 af:=autoflag('N1333_1.split.ms') af.setdata() af.setselect(timerng="2003/05/02/23:50:00 2003/05/03/00:02:00",ant=20) af.run() af.done() ########################################################### ## Image the target source data: ## imgr:=imager(filename='N1333_1.split.ms'); # Create imager tool # Use split out calibrated data. imgr.setdata(mode='channel', nchan=[56,56], start=[1,1], step=[1,1], fieldid=[1:10], spwid=[1,2]); imgr.setimage(nx=640, # Set imaging parameters ny=256, cellx='0.5arcsec', celly='0.5arcsec', stokes='I', mode='channel', # 63+63-30(overlap)=96. Removed 4 end nchan=22, # channels for total of ~88 start=1, # End channels were excluded in split step=4, fieldid=1, spwid=[1,2]); imgr.weight(type="natural",mosaic=T); imgr.setvp(dovp=T); imgr.setoptions(padding=1.0,ftmachine='mosaic') imgr.clean(algorithm='mfclark', # Image and deconvolve niter=100, # with Hogbom CLEAN down to a threshold gain=0.1, threshold='0.02Jy', model='N1333_1.split.mod', image='N1333_1.split.im', # residual='N1333_1.split.resid', displayprogress=T); # ## Display the image img := image('N1333_1.split.im'); img.view(); img.done(); ################################################### ### Make moment maps ## im:=image('N1333_1.split.im') im_mom0:=im.moments(outfile='N1333_1.mom0.im', moments=0,axis=4, mask='indexin(4,[2:21])', #This will limit channels includepix=[0.01,1]); #Give cutoff values im_mom1:=im.moments(outfile='N1333_1.mom1.im', moments=1,axis=4, mask='indexin(4,[2:21])', #This will limit channels includepix=[0.01,1]); #Give cutoff values im.mom0.done() im.mom1.done() # Total intensity with no cutoff im_mom0:=im.moments(outfile='N1333_1.mom0NO.im', moments=0,axis=4, mask='indexin(4,[2:21])'); #This will limit channels im.done() dv.gui() ## Write images to a fits file: ## im:=image(infile="N1333_1.split.im"); im.tofits(outfile="N1333_1.split.fits"); im.done(); im:=image(infile="N1333_1.mom0.im"); im.tofits(outfile="N1333_1.mom0.fits"); im.done(); im:=image(infile="N1333_1.mom1.im") im.tofits(outfile="N1333_1.mom1.fits"); im.done(); im:=image(infile="N1333_1.mom0NO.im") im.tofits(outfile="N1333_1.mom0NO.fits"); im.done();