#!/bin/csh -f
# End-2-end reduction of data observed on 2004-08-30: Cepha Star Formation 
# Regions at 321 GHz. PI: Nimesh Patel
# lsb continuum 
# this script was used to test the basic software components in
# miriad for the reduction of SMA data (on rglinux4).
# jhz 2005-02-02
# Miriad Tasks:
# smalod,
# prthd,
# uvindex,
# uvplt,
# smauvplt,
# uvflag,
# smafix,
# varplt,
# smavarplt,
# uvspec,
# smauvspec,
# mfcal,
# gpplt,
# smamfcal,
# smagpplt,
# blflag,
# smablflag,
# uvaver,
# uvsplit,
# invert,
# mfclean,
# restor,
# cgdisp,
# imstat,
# selfcal.
# 
goto continue
#
#run
#load data and convert mir format to miriad format.
#the spectral channel data are resampled to 16 channels per
#chunk across the 24 chunks.
#
\rm -r cepha.c_rx1.lsb
smalod in=/sma/040830_04:59:20/ \
        out=cepha.c rsnchan=16 rxif=1 sideband=0 nscans=4
prthd in=cepha.c_rx1.lsb
uvindex vis=cepha.c_rx1.lsb
#
#inspect data: uvplt and smauvplt
#
uvplt vis=cepha.c_rx1.lsb axis=time,ampl device=/xs nxy=1,3
sleep 3
smauvplt vis=cepha.c_rx1.lsb axis=time,ampl device=/xs nxy=1,3
#
#uvflag: flag selected source and time range.
#
uvflag vis=cepha.c_rx1.lsb select='source(cep*)','time(0:0,6:0)' \
         flagval=flag
#
# plot the data
#
smauvplt vis=cepha.c_rx1.lsb axis=time,ampl device=/xs nxy=1,3
#
#tsys correction
#
\rm -r cepha.c_rx1.lsb.tsys
smafix vis=cepha.c_rx1.lsb out=cepha.c_rx1.lsb.tsys device=/xs xaxis=time \
        yaxis=systemp bant=6 gant=1 rant=10 nxy=4,2 \
        yrange=0,1200 rmsflag=2 dofit=2 \
        options=tsyscorr,dosour
#
#check tsys with varplt and smavarplt
#
varplt vis=cepha.c_rx1.lsb.tsys device=/xs xaxis=time \
        yaxis=systmp nxy=4,2 yrange=0,1200
sleep 3
smavarplt vis=cepha.c_rx1.lsb.tsys device=/xs xaxis=time \
        yaxis=systmp nxy=4,2 yrange=0,1200
#
# inspect spectra with uvspec and smauvspec
#
uvspec vis=cepha.c_rx1.lsb.tsys interval=15 \
        hann=, axis=freq,ampl device=/xs nxy=1,5
sleep 3
smauvspec vis=cepha.c_rx1.lsb.tsys interval=15 \
        hann=, axis=freq,ampl device=/xs nxy=1,5
sleep 3
smauvspec vis=cepha.c_rx1.lsb.tsys select='source(sat*)' interval=15 \
        hann=, axis=freq,both device=/xs nxy=1,2
#
# bandpass calibration with Saturn
# with mfcal and smamfcal
\rm cepha.c_rx1.lsb.tsys/bandpass cepha.c_rx1.lsb.tsys/gains
cepha.c_rx1.lsb.tsys/freqs
mfcal vis=cepha.c_rx1.lsb.tsys select=source\(sat\*\) refant=1 \
        interval=160 
#check bandpass: ante 6 lower S/N
gpplt vis=cepha.c_rx1.lsb.tsys device=/xs yaxis=amp,phase \
        options=bandpass  \
        nxy=2,4
\rm cepha.c_rx1.lsb.tsys/bandpass cepha.c_rx1.lsb.tsys/gains
cepha.c_rx1.lsb.tsys/freqs
smamfcal vis=cepha.c_rx1.lsb.tsys select=source\(sat\*\) refant=1 \
	interval=160 weight=2
#check bandpass: ante 6 lower S/N
smagpplt vis=cepha.c_rx1.lsb.tsys device=/xs yaxis=amp,phase \
        options=bandpass,opolyfit  \
        polyfit=5 nxy=2,4
#
# flag ante 6
#
uvflag vis=cepha.c_rx1.lsb.tsys \
        select='ante(6)' flagval=flag
\rm cepha.c_rx1.lsb.tsys/bandpass cepha.c_rx1.lsb.tsys/gains
cepha.c_rx1.lsb.tsys/freqs
#
# re-do bandpass
#
smamfcal vis=cepha.c_rx1.lsb.tsys select=source\(sat\*\) refant=1 \
        interval=160 weight=2
#check bandpass: ante 6 lower S/N
smagpplt vis=cepha.c_rx1.lsb.tsys device=/xs yaxis=amp,phase \
        options=bandpass,opolyfit  \
        polyfit=5 nxy=2,4
#
#baseline-based interactive flagging with blflag and smablflag
#smablflag, to flag bad integrations on phase calibrators
#
blflag vis=cepha.c_rx1.lsb.tsys select='source(bll*,0102*)' \
        device=/xs axis=time,ampl options=nocal
smablflag vis=cepha.c_rx1.lsb.tsys select='source(bll*,0102*)' \
        device=/xs axis=time,ampl options=nocal
#
#flag the edge channels 
#
uvflag vis=cepha.c_rx1.lsb.tsys \
        edge=2,2 flagval=flag
#
#check the vis data of the calibrators
#
smauvplt vis=cepha.c_rx1.lsb.tsys axis=time,real \
        select='source(bll*,0102*)' \
        device=/xs nxy=2,2  options=nocal
smauvplt vis=cepha.c_rx1.lsb.tsys axis=time,imag select='source(bll*,0102*)' \
        device=/xs nxy=2,2  options=nocal
#
#calibrate the antenna gains
#
smamfcal vis=cepha.c_rx1.lsb.tsys select='source(bll*,0102*)' \
        refant=1 interval=2 weight=1 \
        options=nopassol flux=2.0 
#check the gain solutions
smagpplt vis=cepha.c_rx1.lsb.tsys device=/xs yaxis=amp,phase \
        options=gain,opolyfit  \
        polyfit=6 nxy=2,4
#check the calibrated data
smauvplt vis=cepha.c_rx1.lsb.tsys axis=time,ampl select='source(bll*,0102*)' \
        device=/xs nxy=2,2  options=, 
smauvplt vis=cepha.c_rx1.lsb.tsys axis=time,ampl select='source(cep*)' \
        device=/xs nxy=2,2  options=,
#flag bad integrations on cepha
smablflag vis=cepha.c_rx1.lsb.tsys select='source(cep*)' \
        device=/xs axis=time,ampl
#apply the gains and bandpass to the targets
\rm -r cepha.c_rx1.lsb.tsys.cal cepha.c_rx1.lsb.tsys.cal.ch0
uvaver vis=cepha.c_rx1.lsb.tsys out=cepha.c_rx1.lsb.tsys.cal 
#split cepha data from the multiple source files
\rm -r cepha.322272
uvsplit vis=cepha.c_rx1.lsb.tsys.cal select='source(cep*)' \
        options=nowindow
#check the spectra to find the channel number on the spectral
#lines.
smauvspec vis=cepha.322272 interval=1000 options=nobase,avall\
        hann=, axis=freq,ampl device=/xs nxy=1,1
sleep 5
smauvspec vis=cepha.322272 interval=1000 options=nobase,avall\
        hann=, axis=chan,ampl device=/xs nxy=1,1 \
        line=channel,44,300,1
#
#line filter: using the linetype in smauvflag to filter the lines
# this is not a good way to do.
#
uvflag vis=cepha.322272 \
        line=channel,3,19,1 flagval=flag
uvflag vis=cepha.322272 \
        line=channel,5,194,1 flagval=flag
uvflag vis=cepha.322272 \
        line=channel,6,214,1 flagval=flag
uvflag vis=cepha.322272 \
        line=channel,6,317,1 flagval=flag
#
#Make a pseudo continuum or ch0 in each chunk
#
\rm -rf cepha.c_rx1.lsbh0
uvaver vis=cepha.322272 out=cepha.c_rx1.lsbh0 \
        line=channel,24,2,13,16
continue:
set i=0
#make map
\rm -rf cepha.cal
\cp -r cepha.c_rx1.lsbh0 cepha.cal
mapping:
#make dirty image
\rm -rf cepha.map cepha.beam cepha.icmp
invert vis=cepha.cal map=cepha.map \
        beam=cepha.beam \
        imsize=512,512 cell=.125 sup=0 options=mfs,sdb,systemp \
        line=channel,8,1,3
#clean map
mfclean map=cepha.map beam=cepha.beam out=cepha.icmp gain=0.1 \
        cutoff=0 niters=500, region='boxes(200,200,312,312)'
\rm -rf cepha.icln
#restore the image with a gaussian beam
restor model=cepha.icmp beam=cepha.beam map=cepha.map \
        out=cepha.icln
#display the image
cgdisp in=cepha.icln type=pixel region='arcsec,boxes(-10,-10,10,10)' \
        xybin=2 device=/xs \
        nxy=1 labtyp=arcsec,arcsec \
	options=full,beambr,wedge,fiddle,trlab \
	range=0,0,lin,2 \
        lines=1,2 csize=1,1,0.8
#do statistics in an area off the source
imstat in=cepha.icln region=arcsec,'boxes(10,10,20,20)'  \
        beam=cepha.beam \
        axes=ra,dec 
#do statistics in an area on the source
imstat in=cepha.icln region=quarter  \
        beam=cepha.beam \
        axes=ra,dec 
#make a contore map
cgdisp in=cepha.icln type=c region=arcsec,'boxes(-8,-8,8,8)' \
        xybin=2 device=/xs nxy=1 labtyp=hms,dms \
        options=full,beambr,wedge,fiddle,trlab csize=1,1,0.8 \
        slev=a,0.01 levs1='-4,4,5,8,12,17,23,30,38,47,57,68,80,93'
#
#smaselfcal 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=cepha.cal interval=20 model=cepha.icmp \
        options=phase,mfs refant=1, line=channel,8,1,3
# using uvputhd to correct the j2000 coordinates in the uv data
# J2000
echo "&&&&&&&&&&&&&&&&&&&&&&&&&"
echo "& smaselfcal iteration="$i
echo "&&&&&&&&&&&&&&&&&&&&&&&&&"
@ i++
goto mapping
case '1':
selfcal vis=cepha.cal interval=15 model=cepha.icmp \
        options=phase,mfs refant=1, line=channel,6,1,4
echo "&&&&&&&&&&&&&&&&&&&&&&&&&&"
echo "& smaselfcal iteration="$i
echo "&&&&&&&&&&&&&&&&&&&&&&&&&&"
@ i++
goto mapping
case '2':
#Make hardcopy of the contour plots of the final image.
\rm cepha_cont.ps
cgdisp in=cepha.icln type=c region=arcsec,'boxes(-8,-8,8,8)' \
        xybin=2 device=cepha_cont.ps/vps nxy=1 labtyp=hms,dms \
        options=full,beambr,wedge,fiddle,trlab csize=1,1,0.8 \
        slev=a,0.08 levs1='-4,4,5,8,12,17,23,30,38,47,57,68,80,93'
#done!
exit
#option for self-calibration on amplitude       
selfcal vis=cepha.cal interval=30 model=cepha.icmp \
        options=amp,mfs refant=1, line=channel,6,1,4
echo "&&&&&&&&&&&&&&&&&&&&&&&"
echo "& smaselfcal iteration="$i
echo "&&&&&&&&&&&&&&&&&&&&&&&"
@ i++
goto mapping
case 3:
exit

