next up previous
Next: TEST THE BASIC Miriad Up: End-to-End Data Reduction Previous: MFS CONTINUUM IMAGE: contmscript.csh

MAKE A SPECTRAL LINE CUBE: maserscript.csh

This is a C-shell script to reduce the 658 GHz data observed on April 8/9, 2004. The script covers a procedure of reducing spectral line data observed with the SMA in Miriad :

#!/bin/csh -f
# End-2-end reduction of data observed on 2004-04-8/9:
# H2O masers at 658 GHz for several sources.
# P.I. 690GHz campaign in 2004 spring.
# The script was used for preparation of Miriad Users Guide
# for the reduction of SMA data.
# jhz 2005-02-02
#goto continue
#continue
#
#load data and convert mir format to miriad format
#the original data has 128 channels per chunk
#smalod resamples the data to 64 channels per chunk
#
\rm -r h2o.lsb h2o.lsb.tsys
smalod in=/home/miriad/SMAdata/040408:08:55:31/ \
         rsnchan=64 out=h2o sideband=0 nscans=10,1430
#inspect data
smauvplt vis=h2o.lsb axis=time,ampl device=/xs nxy=5,3
#tsys correction
\rm -r h2o.lsb.tsys
smafix vis=h2o.lsb out=h2o.lsb.tsys device=/xs xaxis=antel \
       yaxis=systemp nxy=3,2 yrange=0,15000 rmsflag=2 dofit=2 \
       options=tsyscorr,tsysswap
smavarplt vis=h2o.lsb.tsys device=/xs xaxis=antel \
       yaxis=systemp nxy=3,2
#inspect data again
smauvplt vis=h2o.lsb.tsys axis=time,ampl device=/xs nxy=5,3
#inspect spectral data
smauvspec vis=h2o.lsb.tsys \
	select='source(why*)','time(9:30,10:20)' \
        interval=100 \
        hann=7 axis=freq,ampl device=/xs nxy=5,3
#flag data: ant 2 are flagged
uvflag vis=h2o.lsb.tsys \
        select='ant(2)' flagval=flag
#determine bandpass.  
#Using weight=3 because calculation of 
#sigma from original systemp corrupts.
smamfcal vis=h2o.lsb.tsys select='source(ve*,ne*)'  \
        weight=3 refant=3  polyfit=0,0 interval=60
#
#check bandpass
#
smagpplt vis=h2o.lsb.tsys device=fig4.5.ps/vcps \
	yaxis=amp,phase \
        options=bandpass,opolyfit \
        polyfit=2 nxy=1,5
#
#flux scale and gains. The flux density of 53.4 Jy for 
#Callisto is assumed. Also, no phase calibrator is included 
#in the observation. 
mfcal vis=h2o.lsb.tsys select='source(cal*)' \
	refant=3 interval=1000 \
        options=nopassol flux=53.4
#
#smablflag, to flag bad integrations
smablflag vis=h2o.lsb.tsys device=/xs axis=time,ampl
#
#apply and select those spectral windows near/include 
#the masers only apply the gain scale not bpass and 
#time-dependent gains applied.
#
\rm -r h2o.lsb.tsys.cal
uvaver vis=h2o.lsb.tsys out=h2o.lsb.tsys.cal \
        select='window(11,12,13,14,15,16,17,18)'
#
#uvflag vis=h2o.lsb.tsys.cal  flagval=flag
#breaks the file into single source
#
\rm -r *.658128
uvsplit vis=h2o.lsb.tsys.cal options=nowindow
#
#plot the spectra of individual sources and make a hard copy
#
smauvspec vis=whya.658128 \
        hann=7 interval=1000 axis=freq,ampl \
	device=why.ps/vcps nxy=2,5
#
#plot the spectra of individual sources on an x window device.
#
smauvspec vis=whya.658128 \
        hann=7 interval=1000 axis=freq,ampl \
	device=/xs nxy=2,5
#
# cp whya data to a new one which will be further corrected 
# for the residual errors using selfcal. In fact, in this 
# observations, there were no adequate calibrator for antenna 
# gains. We use the maser spot for calculating the antenna 
# gains.
#
\rm -r whya.selfcal
\cp -r whya.658128 whya.selfcal
#
#selfcal phase with point source model
#
selfcal vis=whya.selfcal select='window(4)' \
        interval=10 options=pha \
        refant=3 line=channel,16,24,1
#
# make dirty map
#
\rm -r whya.map whya.beam whya.icmp
invert vis=whya.selfcal map=whya.map beam=whya.beam \
        imsize=512,512 cell=.2 sup=0 select='window(4)' \
        line=channel,16,24,1
# clean map
clean map=whya.map beam=whya.beam out=whya.icmp gain=0.1 \
        cutoff=0 niters=500 \
	region='boxes(225,225,287,287)(1,16)'
# restore image
\rm -r whya.icln
restor model=whya.icmp beam=whya.beam map=whya.map \
        out=whya.icln
cgdisp in=whya.icln type=c xybin=1,1 device=/xs \
        nxy=4,4 labtyp=arcsec,arcsec,absghz \
        options=full,beambr,wedge,trlab,3val \
        lines=1,2 csize=.8,.8,0.8 slev=a,200 \
        levs1=-20,20,30,40,50,60,70,80,90 \
        region='arcsec,boxes(-8,-8,8,8)(1,16)'
selfcal vis=whya.selfcal select='window(4)' \
        interval=7 options=pha \
        refant=3 line=channel,16,24,1 model=whya.icmp
# make dirty map again
\rm -r whya.map whya.beam whya.icmp
invert vis=whya.selfcal map=whya.map beam=whya.beam \
        imsize=512,512 cell=.2 sup=0 select='window(4)' \
        line=channel,16,24,1
# clean map
clean map=whya.map beam=whya.beam out=whya.icmp gain=0.1 \
        cutoff=0 niters=500 \
	region='boxes(225,225,287,287)(1,16)'
# restore image
\rm -r whya.icln
restor model=whya.icmp beam=whya.beam map=whya.map \
        out=whya.icln
cgdisp in=whya.icln type=c xybin=1,1 device=/xs \
        nxy=4,4 labtyp=arcsec,arcsec,absghz \
        options=full,beambr,wedge,trlab,3val \
        lines=1,2 csize=.8,0.8,0.8 slev=a,200 \
        levs1=-20,20,30,40,50,60,70,80,90 \
        region=arcsec,'boxes(-8,-8,8,8)(1,16)'

#selfcal amp with clean point model
selfcal vis=whya.selfcal select='window(4)' \
        interval=12 options=amp \
        refant=3 line=channel,16,24,1 model=whya.icmp
# make dirty map
\rm -r whya.map whya.beam whya.icmp
invert vis=whya.selfcal map=whya.map beam=whya.beam \
        imsize=512,512 cell=.1 sup=0 select='window(4)' \
        line=channel,16,24,1
# clean map
clean map=whya.map beam=whya.beam out=whya.icmp gain=0.1 \
        cutoff=0 niters=500 \
	region='boxes(225,225,287,287)(1,16)'
# restore image
\rm -r whya.icln
restor model=whya.icmp beam=whya.beam map=whya.map \
        out=whya.icln
# make contour images from the image cube
cgdisp in=whya.icln,whya.icln type=p,c \
        xybin=1,1 device=fig4.8.ps/vcps \
        nxy=4,4 labtyp=arcsec,arcsec \
        options=full,beambr,wedge,trlab,3val \
        lines=1,2 csize=.8,.8,0.8 slev=a,1000 \
        levs1=-30,30,40,50,60,70,80,90 \
        region=arcsec,'boxes(-8,-8,8,8)(1,16)' \
        range=20000,50000,lin,2 cols1=7
# do statistics on the image cubes
imstat in=whya.icln region='boxes(10,10,50,50)' \
        beam=whya.beam device=/xs
exit
The Miriad command lines can be written in other shell scripts.

Jun-Hui Zhao (miriad for SMA)
2012-07-09