From CARMA
#! /bin/csh
#########################################
# Original example script by A. Bolatto #
# please detail changes here #
#########################################
set FILE="c0001.n604_coC.1.miriad"
set SRC="NGC604"
set CAL1="0205+322"
set CAL2="0237+288"
set PBCAL="3C454.3"
set NOISE="NOISE"
set FLUX="URANUS"
set WIDE="channel,1,1,15,1"
set LINE="velocity,63,-317.521,2.54,2.54"
set CAL=$CAL2
set OCAL=$CAL1
set RESTFREQ="115.271202"
set REFA=9
# "uvlist vis=$FILE options=spectra" produces the following results
# rest frequency : 115.27120 115.27120 115.27120 115.27120 115.27120 115.27120
# starting channel : 1 16 79 142 157 220
# number of channels : 15 63 63 15 63 63
# starting frequency : 111.47148 111.08239 111.05307 114.93370 115.32280 115.35211
# frequency interval : -0.03125 -0.00049 -0.00049 0.03125 0.00049 0.00049
# starting velocity : 9854.752 10866.791 10943.046 849.513 -162.527 -238.782
# ending velocity : 10992.691 10945.532 11021.788 -288.426 -241.268 -317.523
# velocity interval : 81.274 1.270 1.270 -81.274 -1.270 -1.270
goto continue
# move the "continue:" label as the reduction progresses
continue:
# split data into a dataset per source
foreach i ($SRC $CAL1 $CAL2 $PBCAL $NOISE $FLUX)
rm -rf $i
uvcat vis=$FILE out=$i select="-auto,source("$i")"
end
########################
# look at data quality #
########################
# look at the phase vs time: it should have a smooth trend
uvplt vis=$CAL device=/xs line=$WIDE axis=time,phase
# look at the amplitude vs. time: it should be constant and the same
# for all antennas
uvplt vis=$CAL device=/xs line=$WIDE axis=time,amp
# look at the source phase vs time: if it has no continuum it
# should be just random. Any trends would be due to false fringes
uvplt vis=$SRC device=/xs line=$WIDE axis=time,phase
# flag some suspicious data
uvflag vis=$CAL,$PBCAL,$SRC flagval=flag \
select="ant(5),time(21:30:00,22:15:00)"
goto end
########################
# passband calibration #
########################
## erase product datasets and calibration tables
rm -rf $NOISE.conj $PBCAL.pb $CAL.pb
gpedit vis=$PBCAL options=flag
gpedit vis=$CAL options=flag
## conjugate noise source LSB and replace it in USB
## needs MARCH07 version of uvcal
#uvcal vis=$NOISE out=$NOISE.conj options=noisesrc
## find passband in IF, copy it to PBCAL
#mfcal vis=$NOISE.conj line="channel,282,1,1,1" interval=9999 refant=$REFA
#gpcopy vis=$NOISE.conj out=$PBCAL options=nocal
## this doesn't fix most problems
##uvspec vis=$PBCAL axis=chan,phase line="channel,282,1,1,1" device=/xs \
# interval=99999 yrange=-180,180
# use PBCAL to derive astronomical passband calibration
mfcal vis=$PBCAL line="channel,282,1,1,1" interval=9999 refant=$REFA
# this works out considerably better
uvspec vis=$PBCAL axis=chan,phase line="channel,282,1,1,1" device=/xs \
interval=99999 yrange=-180,180
# copy passband table from PBCAL to CAL
gpcopy vis=$PBCAL out=$CAL options=nocal
# create new dataset with calibration applied, otherwise linetype averaging
# does not work properly. Use all wideband channels.
uvcat vis=$CAL out=$CAL.pb select="window(1,4)"
goto end
###############################
# amplitude/phase calibration #
###############################
rm -rf $FLUX.pb
gpcopy vis=$PBCAL out=$FLUX options=nocal
uvcat vis=$FLUX out=$FLUX.pb select="window(1,4)"
# not clear that I need this step (bootflux probably already does it)
# it won't hurt
selfcal vis=$FLUX.pb options=apriori,amp,noscale interval=0.1 \
line="channel,1,1,30,1" refant=$REFA
bootflux vis=$FLUX.pb,$CAL.pb primary=$FLUX \
line="channel,1,1,30,1" taver=9999
goto end
# selfcalibrate phase cal, with passband calibration applied
# imposing flux found by bootflux solution
selfcal vis=$CAL.pb line="channel,1,1,30,1" options=amp,noscale,apriori \
flux=1.2 interval=20 refant=$REFA
# every channel should have zero phase
uvspec vis=$CAL.pb axis=chan,phase line="channel,272,1,10,1" device=/xs \
interval=99999 yrange=-180,180
# show time series of selfcalibrated wideband channels
uvplt vis=$CAL.pb axis=time,phase device=/xs line="channel,1,16,15,1"
uvplt vis=$CAL.pb axis=time,amp device=/xs line="channel,1,16,15,1"
goto end
##########################
# phase-only calibration #
##########################
# selfcalibrate phase cal, with passband calibration applied
# most of the time the online amplitude calibration seems very good...
# Note that we are averaging over all wideband channels.
# Channel linetype averaging does not weight by
# bandwidths and/or Tsys. This is why we split out only the continuum
# windows.
selfcal vis=$CAL.pb line="channel,1,1,30,1" interval=20 refant=$REFA
# every channel should have zero phase
uvspec vis=$CAL.pb axis=chan,phase line="channel,272,1,10,1" device=/xs \
interval=99999 yrange=-180,180
# show time series of selfcalibrated wideband channels
uvplt vis=$CAL.pb axis=time,phase device=/xs line="channel,1,16,15,1"
goto end
###################################
# sanity checks after calibration #
###################################
# plot antenna gain tables
# amplitude gain>1 indicates phasecal was weaker than expected
# probably caused by mispointing
gpplt vis=$CAL.pb device=/xs yaxis=amp
# phase gain should be smooth
gpplt vis=$CAL.pb device=/xs yaxis=phase
# look at phase vs time after selfcal: it should be centered around zero
uvplt vis=$CAL.pb device=/xs axis=time,phase line=$WIDE
# look at phase vs baseline length to assess atmospheric decorrelation
# it should be flaring at the longer baselines
uvplt vis=$CAL.pb device=/xs line=$WIDE axis=uvdist,phase options=nobase
# look at amplitude vs. time: it should be about was I set it to
# be in the selfcal if I did amplitude selfcal
uvplt vis=$CAL.pb device=/xs line=$WIDE axis=time,amp
goto end
#######################
# image the phase cal #
#######################
# invert the visibilities to obtain a map of the calibrator
rm -rf $CAL.map $CAL.beam
invert vis=$CAL.pb map=$CAL.map beam=$CAL.beam line=$WIDE \
sup=0 options=systemp
# report size of calibrator and beam
imfit in=$CAL.map object=beam region="arcsec,box(-15,-15,15,15)"
imfit in=$CAL.beam object=beam region="arcsec,box(-15,-15,15,15)"
# calibrator should be a point source with some sidelobes
cgdisp in=$CAL.map,$CAL.map device=/xs type=pixel,contour \
slev=p,10 levs2=1,2,3,4,5,6,7,8,9
# compare to the beam: it better look very similar
cgdisp in=$CAL.beam \
type=con slev=p,10 levs1=1,2,3,4,5,6,7,8,9 \
labtyp=arcsec options=full,mirror device=/xs
# clean the calibrator map
rm -rf $CAL.cc $CAL.cm
clean map=$CAL.map beam=$CAL.beam out=$CAL.cc gain=0.1 niters=250
restor map=$CAL.map model=$CAL.cc beam=$CAL.beam out=$CAL.cm
# now it should be a very good point source at the 1% level
cgdisp in=$CAL.cm,$CAL.cm device=/xs type=pixel,contour \
slev=p,1 levs2=1,5,10,20,50,90
goto end
#################################
# image the secondary phase cal #
#################################
# check on the second calibrator, to see if the
# passband/phase solution transfers correctly
rm -rf $OCAL.pb
gpedit vis=$OCAL options=flag
gpcopy vis=$CAL out=$OCAL mode=copy options=nocal
gpcopy vis=$CAL.pb out=$OCAL mode=apply options=nopass
uvcat vis=$OCAL out=$OCAL.pb
rm -rf $SRC.map $SRC.beam
invert vis=$OCAL.pb map=$OCAL.map beam=$OCAL.beam line=$WIDE \
sup=0 options=systemp
cgdisp in=$OCAL.map,$OCAL.map device=/xs type=pixel,contour \
slev=p,10 levs2=1,2,3,4,5,6,7,8,9
rm -rf $OCAL.cc $OCAL.cm
clean map=$OCAL.map beam=$OCAL.beam out=$OCAL.cc gain=0.1 niters=250
restor map=$OCAL.map model=$OCAL.cc beam=$OCAL.beam out=$OCAL.cm
cgdisp in=$OCAL.cm,$OCAL.cm device=/xs type=pixel,contour \
slev=p,1 levs2=20,50,90
fits in=$OCAL.cm out=$OCAL.cm.fits op=xyout
goto end
####################
# image the source #
####################
# correct header information
puthd in=$SRC/restfreq value=$RESTFREQ
# copy the passband and phase solutions
# produce new dataset with solutions applied so that
# linetype averaging works properly
rm -rf $SRC.pb
gpcopy vis=$CAL out=$SRC mode=copy options=nocal
gpcopy vis=$CAL.pb out=$SRC mode=copy options=nopass
uvcat vis=$SRC out=$SRC.pb
# there are some flakey channels in this dataset
uvflag vis=$SRC.pb line="channel,1,197,1,1" flagval=flag
uvflag vis=$SRC.pb line="channel,1,262,1,1" flagval=flag
uvflag vis=$SRC.pb line="channel,1,259,1,1" flagval=flag
# now invert the calibrated source visibilities
# for the line windows. Miriad easily handles combining many datasets.
# by default Miriad only images out to the half power point of
# each field. To make them stitch seamlessly I specify imsize and
# cellsize to make sure it images out to the 20% point.
rm -rf $SRC.map $SRC.beam
invert vis=../3/$SRC.pb,../1/$SRC.pb,../4/$SRC.pb map=$SRC.map \
beam=$SRC.beam line=$LINE imsize=160 cell=0.5 \
sup=0 options=systemp,mosaic select="window(5,6)" slop=0.1,interp
# make a continuum map too
rm -rf $SRC.cont
uvcat vis=$SRC.pb out=$SRC.cont select="window(1,4)"
rm -rf $SRC.cmap $SRC.cbeam
invert vis=../1/$SRC.cont,../3/$SRC.cont,../4/$SRC.cont map=$SRC.cmap \
beam=$SRC.cbeam line="channel,1,1,30,1" imsize=160 cell=0.5 \
sup=0 options=systemp,mosaic
fits in=$SRC.cmap out=$SRC.cmap.fits op=xyout
# figure out the average restoring beam for the mosaic
# by inverting the calibrated source visibilities
# without the mosaic option
rm -rf $SRC.bmap $SRC.abeam
invert vis=../3/$SRC.pb,../1/$SRC.pb,../4/$SRC.pb map=$SRC.bmap \
beam=$SRC.abeam line=$LINE imsize=160 cell=0.5 \
sup=0 options=systemp select="window(5,6)" slop=0.1,interp
rm -rf $SRC.bmap
# see what mospsf predicts (less accurate, presumably)
rm -rf $SRC.mospsf
mospsf beam=$SRC.beam out=$SRC.mospsf
goto end
##################################
# CLEANing and wrapping up tasks #
##################################
# standard CLEANing for the source map
# rm -rf $SRC.cc $SRC.cm
# mossdi map=$SRC.map beam=$SRC.beam out=$SRC.cc gain=0.1 \
# niters=100 cutoff=0.2
# restor map=$SRC.map model=$SRC.cc beam=$SRC.beam out=$SRC.cm
mkmask:
# masked CLEANing gives better results for faint sources
# make mask for CLEANing
rm -rf $SRC.conv $SRC.mask
convol map=$SRC.map fwhm=5 out=$SRC.conv
imstat in=$SRC.conv region="arcsec,box(-30,-30,30,30)" device=/xs
#fits in=$SRC.conv out=$SRC.conv.fits op=xyout
# I usually mask at the 2.5 sigma level in the convolved map
maths exp="(<"$SRC".conv>.gt.2.3).and.(abs(x).le.60).and.(abs(y).le.60)" \
out=$SRC.mask xrange=-100,100 yrange=-100,100
fits in=$SRC.mask out=$SRC.mask.fits op=xyout
goto end
# CLEAN the source map using a mask
rm -rf $SRC.cc
mossdi map=$SRC.map beam=$SRC.beam out=$SRC.cc gain=0.1 \
niters=100 cutoff=0.2 region="mask("$SRC.mask")"
rm -rf $SRC.cm
restor map=$SRC.map model=$SRC.cc beam=$SRC.abeam out=$SRC.cm
# Make masked moment map
rm -r $SRC.mom0 $SRC.mcm
maths exp="<"$SRC.cm">*<"$SRC.mask">" out=$SRC.mcm
moment in=$SRC.mcm \
out=$SRC.mom0 mom=0 clip=-1000,-1 region="images(20,45)"
mkfits:
# output it as a .fits file that I can peruse with other viewers
rm -rf $SRC.fits
fits in=$SRC.map out=$SRC.fits op=xyout
fits in=$SRC.cm out=$SRC.cm.fits op=xyout
fits in=$SRC.mom0 out=$SRC.mom0.fits op=xyout
# display/print an image of the channel maps in grayscale with
# the integrated intensity contours overlaid, beam in corner,
# velocity scale indicated, etc
cgdisp in=$SRC.cm,$SRC.mom0 device=/xs type=pix,con \
region="arcsec,box(-30,-30,30,30)(27,42)" labtyp=arcsec \
beamtyp=b,l nxy=2,2 range=0,0.3,lin,1 options=3value \
levs1=0.5,1.5,2.5,3.5,4.5 csize=1,1,0,0
goto end
end:
