#!/bin/csh -f
# The observation was made on 2004-3-29, SgrA*, in a 'continuum mode'
# lower spectral resolution 16 channels per chunk. PI. Zhao
# 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 rxif=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
# display and 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 
continue:
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
sleep 4
#
#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