#!/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,swaptsys 
smavarplt vis=h2o.lsb.tsys device=/xs xaxis=antel \
          yaxis=wsystemp 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