MFS CONTINUUM IMAGE: contmscript.csh

This is a C-shell script to reduce the 340 GHz data of SgrA* observed on March 29, 2004. This is a good exmaple to show a complete procedure in construction of an high dynamic range continuum image. The lower side band data is loaded into Miriad using smalod and Tsys corrections are made using smafix. The numerous problems in the data are identified with the graphic display tools smauvspec and smauvplt. The bad data points are flagged using either uvflag or smablflag. The bandpass solutions based on Callisto are determined using smamfcal and is smoothed using smagpplt. Also the antenna gains are solved based on NRAO 530. Both the bandpass and antenna gain solutions are applied to the visibility data and a `continuum' uv file is created with 'channel-0' (or averaged chunk data) in each of the 24 spectral windows. Using the multi-frequency synthesis (MSF) technique, a final self-calibrated image is created.

In addition, in the observation on March 29, 2004 the position of SgrA* was offset by $\sim1.3$ arcsec from the phase center. This script provides a way to shift Sgr A* to the phase center using selfcal and correct the J2000 coordinates in the uv data using uvputhd.

Comments are interleaved with the command lines in this scripts.

#!/bin/csh -f
# The observation was made on 2004-3-29, SgrA*, in a 
# 'continuum mode' or lower spectral resolution of 16 channels per 
# chunk. Observers: P.I. Zhao and his company.
# The script was used to prepare the Users Guide for the
# Reduction of SMA Data:
# loading, inspection, editing flagging, calibrate and image
# for continuum observations.
# Miriad Tasks :
# smalod,
# prthd,
# uvindex,
# smafix,
# smavarplt,
# smauvspec,
# uvflag,
# smablflag,
# qvack,
# smamfcal,
# smagpplt,
# uvaver,
# uvsplit,
# invert,
# mfclean,
# restor,
# imstat,
# cgdisp,
# puthd,
# selfcal
# jhz: 2005-02-02
goto continue
#continue:
#
#loading the data and converting to miriad format
#
\rm -r gc.lsb
smalod in=/home/miriad/SMAdata/sma/040329_09:37:37/ out=gc \
        ifsel=1 sideband=0 nscans=35,950
# print header
prthd in=gc.lsb
# scans uv data
uvindex vis=gc.lsb
#Tsys corrections
\rm -r gc.lsb.tsys
smafix vis=gc.lsb out=gc.lsb.tsys device=fig2.1.ps/vcps \
        xaxis=time yaxis=systemp nxy=2,4 \
        rmsflag=2 dofit=2 options=tsyscorr,dosour \
        yrange=100,1200
# plot and check variables
smavarplt vis=gc.lsb.tsys device=fig2.2.ps/vcps xaxis=time \
        yaxis=systmp nxy=2,4
# inspect spectral data
smauvspec vis=gc.lsb.tsys interval=30 axis=chan,ampl \
        device=fig2.3.ps/vcps \
        nxy=4,7 select='source(jup*)'
smauvspec vis=gc.lsb.tsys interval=30 axis=chan,phas \
        device=fig2.4.ps/vcps \
        nxy=4,7 select='source(jup*)'
smauvspec vis=gc.lsb.tsys interval=100 axis=freq,phas \
        device=fig2.5.ps/vcps nxy=2,4 \
        select='ant(4)(1,2,3,4,5,6,7,8)','source(cal*)' \
        options='nopass,nocal' select='source(cal*)'
# flag bad data points s1,2,3,4 on baseline 5-8
uvflag vis=gc.lsb.tsys select='window(1,2,3,4)','ant(5)(8)' \
	flagval=flag
# inspect spectral data again
smauvspec vis=gc.lsb.tsys interval=60 axis=freq,ampl \
	device=/xs nxy=2,2 select='source(cal*)'
# inspect vis vs time
smauvplt vis=gc.lsb.tsys stokes=xx axis=time,ampl \
	device=fig2.6.ps/vcps \
        select='ant(4)' yrange='-10,250' nxy=2,4
smauvplt vis=gc.lsb.tsys axis=time,ampl device=/xs \
        nxy=2,4 yrange='-10,250'
# flag sgra-star before ut=10hr
uvflag vis=gc.lsb.tsys select='source(sgra*)','time(0:0,10:0)' \
        flagval=flag
# flag ante 2
uvflag vis=gc.lsb.tsys select='ant(2)' flagval=flag
# plot data
smauvplt vis=gc.lsb.tsys stokes=xx axis=time,ampl \
	device=/xs nxy=2,2
#delete the beginning integration in each scan
qvack vis=gc.lsb.tsys mode=source interval=1
#check the uvplot
smauvplt vis=gc.lsb.tsys stokes=xx axis=time,ampl \
	device=/xs nxy=2,2
# baseline based flag
echo "Please use mouse to  flag the  bad data point now."
echo "left button to mark the bad points, and "
echo "right button to change baseline."
smablflag vis=gc.lsb.tsys device=/xs axis=time,ampl
#
# calculate antenna-based bandpass from Callisto data
#
smamfcal vis=gc.lsb.tsys select='source(cal*)' \
	refant=3 interval=100 \
    	weight=2
# display, interpolate the bandpass
smagpplt vis=gc.lsb.tsys device=/xs yaxis=a,p \
        options=bandpass,opolyfit \
        polyfit=5 nxy=2,4
sleep 5
#
# calculate antenna gains from NRAO 530
#
smamfcal vis=gc.lsb.tsys select='source(nra*)' \
	refant=3 weight=1 \
        interval=1 options=nopassol \
        flux=2
# interpolate gains
smagpplt vis=gc.lsb.tsys device=/xs yaxis=real,imag \
        options=gains,opolyfit \
        polyfit=5 nxy=2,4
# apply bandpass and gains
\rm -r gc.lsb.tsys.cal gc.lsb.tsys.cal.ch0
uvaver vis=gc.lsb.tsys out=gc.lsb.tsys.cal
# average and make ch0 in each spectral windows
uvaver vis=gc.lsb.tsys.cal out=gc.lsb.tsys.cal.ch0 \
        line=channel,24,2,13,16
# split sgra-star data from multiple-source file
\rm -r sgra-star.335866
uvsplit vis=gc.lsb.tsys.cal.ch0 \
        select='source(sgra*)' options=nowindow
\rm -r gc.cal
cp -r sgra-star.335866 gc.cal
# show the uv coverage
uvflag vis=gc.cal select='time(12:35:30,12:40)' flagval=flag
smauvplt vis=gc.cal axis=uu,vv options=nobase device=/xs nxy=1,1
set i=0
sleep 5
mapping:
\rm -r sgra-star.map sgra-star.beam sgra-star.icmp
#
# make dirty map
#
invert vis=gc.cal map=sgra-star.map \
        beam=sgra-star.beam select='time(12:00,19:00)' \
        imsize=512,512 cell=.15 sup=0 options=mfs,sdb,systemp \
        line=channel,12,1,2
#clean map
mfclean map=sgra-star.map beam=sgra-star.beam \
        out=sgra-star.icmp gain=0.1 cutoff=0 \
        niters=500, region='boxes(200,200,312,312)'
\rm -r sgra-star.icln
#restore the image with a gaussian beam
restor model=sgra-star.icmp beam=sgra-star.beam \
	map=sgra-star.map \
        out=sgra-star.icln
#display the image
cgdisp in=sgra-star.icln type=pixel region=quarter xybin=2 \
        device=fig4.3.ps/vcps \
        nxy=1 labtyp=arcsec,arcsec \
	options=full,beambr,wedge,trlab \
        lines=1,2 range=0,0,lin,2 csize=1,1,0.8
#do statistics in an area off the source
imstat in=sgra-star.icln region='boxes(10,10,30,30)' \
        beam=sgra-star.beam \
        axes=ra,dec
#do statistics in an area on the source
imstat in=sgra-star.icln region=quarter \
        beam=sgra-star.beam \
        axes=ra,dec
#make a contore map
cgdisp in=sgra-star.icln type=c \
        region='arcsec,boxes(-15,-15,15,15)' \
        xybin=2 device=/xs nxy=1 labtyp=arcsec,arcsec \
        options=full,beambr csize=1,1,0.8 \
        slev=a,0.02 levs1=-3,3,5,8,12,17,23,30,38, \
	47,57,68,80,93,107,122,138,155,172,191
#selfcal phase only
echo $i
switch ($i)
case '0':
#
# SgrA*'s position was offset by 1.3 arcsec w.r.t. the phase 
# center shift the source to the phase center with a point 
# source model in selfcal
#
selfcal vis=gc.cal interval=5  \
        options=phase,mfs refant=3, line=channel,8,1,3
# using uvputhd to correct the j2000 coordinates in the uv data
# J2000
uvputhd vis=gc.cal hdvar=ra type=d length=1 varval=4.64985164 \
        out=gc.cal
echo "&&&&&&&&&&&&&&&&&&&&&&&&&"
echo "& smaselfcal iteration="$i
echo "&&&&&&&&&&&&&&&&&&&&&&&&&"
@ i++
goto mapping
case '1':
selfcal vis=gc.cal interval=5 model=sgra-star.icmp \
        options=phase,mfs refant=3, line=channel,8,1,3
echo "&&&&&&&&&&&&&&&&&&&&&&&&&&"
echo "& smaselfcal iteration="$i
echo "&&&&&&&&&&&&&&&&&&&&&&&&&&"
@ i++
goto mapping
case '2':
selfcal vis=gc.cal interval=5 model=sgra-star.icmp \
        options=phase,mfs refant=3, line=channel,2,1,12
echo "&&&&&&&&&&&&&&&&&&&&&&&"
echo "& smaselfcal iteration="$i
echo "&&&&&&&&&&&&&&&&&&&&&&&"
@ i++
goto mapping
case 3:
selfcal vis=gc.cal interval=3 model=sgra-star.icmp \
        options=phase,mfs refant=3, line=channel,2,1,12
echo "&&&&&&&&&&&&&&&&&&&&&&&"
echo "& selfcal iteration="$i
echo "&&&&&&&&&&&&&&&&&&&&&&&"
@ i++
goto mapping
case 4:
selfcal vis=gc.cal interval=2 model=sgra-star.icmp \
        options=phase,mfs refant=3, line=channel,1,1,24
echo "&&&&&&&&&&&&&&&&&&&&&&&"
echo "& selfcal iteration="$i
echo "&&&&&&&&&&&&&&&&&&&&&&&"
@ i++
goto mapping
case 5:
selfcal vis=gc.cal interval=10 model=sgra-star.icmp \
        options=amp,mfs refant=3 line=channel,1,1,24
@ i++
smablflag vis=gc.cal device=/xs axis=time,ampl
goto mapping
case 6:
continue:
cgdisp in=sgra-star.icln type=c
        region='arcsec,boxes(-15,-15,15,15)' \
        xybin=1,1 slev=a,0.008 \
        levs1=-5,3,4,6,9,13,18,24,31,41,51,71,91,100,150,300 \
        device=/xs nxy=1,1 labtyp=hms,dms \
        options=full,beambr,wedge,trlab,3val
    exit
endsw