Example 2

From CARMA
Jump to: navigation, search

#! /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:

Personal tools