#!/bin/csh -f
# for the reduction of SMA DB data (standard 4GHz mode).
#
# Responsible: Jun-Hui Zhao (SAO)
# Version 2.3 
# jhz 2010-01-20 for e2e test 
# jhz 2010-03-11 updated
# jhz 2010-05-07 updated
# jhz 2010-08-21 updated
# jhz 2010-10-18 updated
# jhz 2010-12-11 updated
# jhz 2011-01-30 updated
# jhz 2011-03-09 updated add lsb data for continuum
#                separate high- and low- chunks n calibration 
#                for each sideband
# jhz 2011-03-14 updated with identification of SiS
# jhz 2011-03-25 updated for automatically labelling SiS line
# jhz 2011-08-08 added goto TAG
# testing data: 090218_04:25:36 from Nimesh Patel 
#
# demo illustration figures
# fig1.ps = continuum image of vycmasma 
# fig2.ps = spectrum from the vis data of vycmasma
# fig3.ps = spectrum from line image cube vycmasma (SiS)
# fig4.ps = image of integrated line flux vycmasma (SiS)
#
set darea=/sma/SMAusers/jzhao/data/090218_04:25:36
set out=090218
goto SMALOD
#goto INSPECT
#goto SMAFIX
#goto BPASS
#goto GCAL
#goto SMAFLUX
#goto UVLIN0
#goto CMAP
#goto UVLINL
#goto LMCUBE
#goto MOMENT
SMALOD:
echo "1. READ THE MIR DATA INTO MIRIAD"
\rm -r  ${out}_rx0.usb ${out}_rx0.lsb
smalod in=$darea out=${out} rxif=0 sideband=2  
$<
INSPECT:
echo "2. INSPECT THE RAW UV DATA"
uvindex vis=${out}_rx0.usb
uvindex vis=${out}_rx0.lsb
$<
smacheck vis=${out}_rx0.usb var=systemp flagval=flag \
        range=0,500
smacheck vis=${out}_rx0.lsb var=systemp flagval=flag \
        range=0,500
$<
SMAFIX:
echo "3. FRINGE AMPLITUDE CORRECTION"
\rm -r ${out}_rx0.usb.tsys
smafix vis=${out}_rx0.usb out=${out}_rx0.usb.tsys \
	device=/xs xaxis=time \
       	yaxis=systemp nxy=2,4   \
       	options=tsyscorr
\rm -r ${out}_rx0.lsb.tsys
smafix vis=${out}_rx0.lsb out=${out}_rx0.lsb.tsys \
        device=/xs xaxis=time \
        yaxis=systemp nxy=2,4   \
        options=tsyscorr
#editing and inspection
smablflag vis=${out}_rx0.usb.tsys axis=time,ampl device=/xs options=nobase 
smablflag vis=${out}_rx0.lsb.tsys axis=time,ampl device=/xs options=nobase
$<
smauvspec vis=${out}_rx0.usb.tsys \
	select='source(3c273)' \
        interval=1000  \
        axis=freq,both device=/xs nxy=1,5
$<
smauvspec vis=${out}_rx0.usb.tsys \
        select='source(titan)' \
        interval=1000  \
        axis=freq,both device=/xs nxy=1,5
BPASS:
echo "4. BANDPASS CALIBRATIONS"
$<
\rm -r ${out}_rx0.usb.tsys.if1 ${out}_rx0.usb.tsys.if2
smasplt vis=${out}_rx0.usb.tsys  out=${out}_rx0.usb.tsys.if1 options=if1
smasplt vis=${out}_rx0.usb.tsys  out=${out}_rx0.usb.tsys.if2 options=if2
\rm -r ${out}_rx0.lsb.tsys.if1 ${out}_rx0.lsb.tsys.if2
smasplt vis=${out}_rx0.lsb.tsys  out=${out}_rx0.lsb.tsys.if1 options=if1
smasplt vis=${out}_rx0.lsb.tsys  out=${out}_rx0.lsb.tsys.if2 options=if2 
$<
smamfcal vis=${out}_rx0.usb.tsys.if1 select='source(3c273)' edge=5,5  \
        weight=2 refant=3 interval=60000
$<
smagpplt vis=${out}_rx0.usb.tsys.if1  device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=1,7
$<
smamfcal vis=${out}_rx0.usb.tsys.if2 select='source(3c273)' edge=5,5  \
        weight=2 refant=3 interval=60000
$<
smagpplt vis=${out}_rx0.usb.tsys.if2  device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=1,7
$<
smamfcal vis=${out}_rx0.lsb.tsys.if1 select='source(3c273)' edge=5,5  \
        weight=2 refant=3 interval=60000
$<
smagpplt vis=${out}_rx0.lsb.tsys.if1  device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=1,7
$<
smamfcal vis=${out}_rx0.lsb.tsys.if2 select='source(3c273)' edge=5,5  \
        weight=2 refant=3 interval=60000
$<
smagpplt vis=${out}_rx0.lsb.tsys.if2  device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=1,7
$<
\rm -r ${out}_rx0.usb.tsys.if1.bp ${out}_rx0.usb.tsys.if2.bp
uvaver vis=${out}_rx0.usb.tsys.if1 options=nocal \
        out=${out}_rx0.usb.tsys.if1.bp
uvaver vis=${out}_rx0.usb.tsys.if2 options=nocal \
        out=${out}_rx0.usb.tsys.if2.bp
$<
\rm -r ${out}_rx0.lsb.tsys.if1.bp ${out}_rx0.lsb.tsys.if2.bp
uvaver vis=${out}_rx0.lsb.tsys.if1 options=nocal \
        out=${out}_rx0.lsb.tsys.if1.bp
uvaver vis=${out}_rx0.lsb.tsys.if2 options=nocal \
        out=${out}_rx0.lsb.tsys.if2.bp
$<
uvflag vis=${out}_rx0.usb.tsys.if2.bp \
        line=chan,1,2337,1 \
        flagval=flag
uvflag vis=${out}_rx0.lsb.tsys.if2.bp \
        line=chan,1,2337,1 \
        flagval=flag
GCAL:
echo "5. GAIN CALIBRATION"
smamfcal vis=${out}_rx0.usb.tsys.if1.bp  \
select='source(0730-116)'  \
refant=3 interval=2 options=nopassol
$<
smagpplt vis=${out}_rx0.usb.tsys.if1.bp device=/xs \
        yaxis=amp,pha  \
        options=gains  polyfit=6 \
        nxy=2,4
$<
smamfcal vis=${out}_rx0.usb.tsys.if2.bp  \
select='source(0730-116)'  \
refant=3 interval=2 options=nopassol
$<
smagpplt vis=${out}_rx0.usb.tsys.if2.bp device=/xs \
        yaxis=amp,pha  \
        options=gains  polyfit=6 \
        nxy=2,4
$<
smamfcal vis=${out}_rx0.lsb.tsys.if1.bp  \
select='source(0730-116)'  \
refant=3 interval=2 options=nopassol
$<
smagpplt vis=${out}_rx0.lsb.tsys.if1.bp device=/xs \
        yaxis=amp,pha  \
        options=gains  polyfit=6 \
        nxy=2,4
$<
smamfcal vis=${out}_rx0.lsb.tsys.if2.bp  \
select='source(0730-116)'  \
refant=3 interval=2 options=nopassol
$<
smagpplt vis=${out}_rx0.lsb.tsys.if2.bp device=/xs \
        yaxis=amp,pha  \
        options=gains  polyfit=6 \
        nxy=2,4
$<
\rm -r ${out}_rx0.usb.tsys.if1.bp.g ${out}_rx0.usb.tsys.if2.bp.g
uvaver vis=${out}_rx0.usb.tsys.if1.bp out=${out}_rx0.usb.tsys.if1.bp.g
uvaver vis=${out}_rx0.usb.tsys.if2.bp out=${out}_rx0.usb.tsys.if2.bp.g
\rm -r ${out}_rx0.lsb.tsys.if1.bp.g ${out}_rx0.lsb.tsys.if2.bp.g
uvaver vis=${out}_rx0.lsb.tsys.if1.bp out=${out}_rx0.lsb.tsys.if1.bp.g
uvaver vis=${out}_rx0.lsb.tsys.if2.bp out=${out}_rx0.lsb.tsys.if2.bp.g
$<
SMAFLUX:
echo "6. BOOTSTRAPE THE FLUX DENSITY SCALE"
smaflux vis=${out}_rx0.usb.tsys.if1.bp.g \
	select='source(titan)' mirhome=$MIR  
$<
smaflux vis=${out}_rx0.usb.tsys.if2.bp.g \
        select='source(titan)' mirhome=$MIR
$<
smaflux vis=${out}_rx0.lsb.tsys.if1.bp.g \
        select='source(titan)' mirhome=$MIR
$<
smaflux vis=${out}_rx0.lsb.tsys.if2.bp.g \
        select='source(titan)' mirhome=$MIR
$<
\rm -r ${out}_rx0.usb.tsys.if1.bp.g.av
uvaver vis=${out}_rx0.usb.tsys.if1.bp.g out=${out}_rx0.usb.tsys.if1.bp.g.av 
$<
\rm -r ${out}_rx0.usb.tsys.if2.bp.g.av
uvaver vis=${out}_rx0.usb.tsys.if2.bp.g out=${out}_rx0.usb.tsys.if2.bp.g.av
$<
\rm -r ${out}_rx0.lsb.tsys.if1.bp.g.av
uvaver vis=${out}_rx0.lsb.tsys.if1.bp.g out=${out}_rx0.lsb.tsys.if1.bp.g.av
$<
\rm -r ${out}_rx0.lsb.tsys.if2.bp.g.av
uvaver vis=${out}_rx0.lsb.tsys.if2.bp.g out=${out}_rx0.lsb.tsys.if2.bp.g.av
$<
UVLIN0:
echo "7. LINE-FREE CONTINUUM"
$<
\rm -r ${out}_rx0.usb.tsys.if1.bp.g.av.ch0 
uvlin vis=${out}_rx0.usb.tsys.if1.bp.g.av chans=1,3072 out=${out}_rx0.usb.tsys.if1.bp.g.av.ch0 \
        order=1 mode=chan0 options=nowin
$<
\rm -r ${out}_rx0.usb.tsys.if2.bp.g.av.ch0
uvlin vis=${out}_rx0.usb.tsys.if2.bp.g.av chans=1,3072 out=${out}_rx0.usb.tsys.if2.bp.g.av.ch0 \
        order=1 mode=chan0 options=nowin
$<
\rm -r ${out}_rx0.lsb.tsys.if1.bp.g.av.ch0
uvlin vis=${out}_rx0.lsb.tsys.if1.bp.g.av chans=1,3072 out=${out}_rx0.lsb.tsys.if1.bp.g.av.ch0 \
        order=1 mode=chan0 options=nowin
$<
\rm -r ${out}_rx0.lsb.tsys.if2.bp.g.av.ch0
uvlin vis=${out}_rx0.lsb.tsys.if2.bp.g.av chans=1,3072 out=${out}_rx0.lsb.tsys.if2.bp.g.av.ch0 \
        order=1 mode=chan0 options=nowin
CMAP:
echo "8. CONTINUUM IMAGE"
$<
\rm -r vycmasma.map vycmasma.beam vycmasma.icmp
invert vis=${out}_rx0.usb.tsys.if1.bp.g.av.ch0,${out}_rx0.usb.tsys.if2.bp.g.av.ch0,${out}_rx0.lsb.tsys.if1.bp.g.av.ch0,${out}_rx0.lsb.tsys.if2.bp.g.av.ch0  map=vycmasma.map beam=vycmasma.beam \
        imsize=512,512 cell=0.15 robust=0 \
        options=systemp select='source(vycmasma)'
$<
cgdisp in=vycmasma.map type=pixel region=arcsec,'boxes(-15,-15,15,15)' \
        xybin=1,1 device=/xs nxy=1,1 options=full,beambr,wedge,trlab,3val \
        labtyp=arcsec,arcsec range=0,0,lin,2 cols1=7 csize=0.5,0.5,0.5
$<
clean map=vycmasma.map beam=vycmasma.beam out=vycmasma.icmp gain=0.08 \
        cutoff=0.001 niters=5000 \
        region=arcsec,'boxes(-1,-1,1,1)'
$<
\rm -r vycmasma.icln
restor model=vycmasma.icmp beam=vycmasma.beam map=vycmasma.map  \
        out=vycmasma.icln
$<
cgdisp in=vycmasma.icln type=pixel region=arcsec,'boxes(-4,-4.5,5,4.5)' \
        xybin=1,1 device=fig1.ps/cps nxy=1,1 options=beambr,wedge,blacklab \
        labtyp=arcsec range=0.01,0.18,lin,2 cols1=7 csize=1,1
$<
#...
#skip
#select  58001    CCS-34        (<type> y to select it or n to skip it or t to terminate)
#n
#skip
#select  58002    NaCl          (<type> y to select it or n to skip it or t to terminate)
#n
#skip
#select  60001    OCS           (<type> y to select it or n to skip it or t to terminate)
#n
#skip
#select  60002    SiS           (<type> y to select it or n to skip it or t to terminate)
#y
#yes
#select  60003    CH3OCHO-A     (<type> y to select it or n to skip it or t to terminate)
#t
#terminate selection
#...
smauvspec vis=${out}_rx0.lsb.tsys.if1.bp.g.av hann=11 \
        select='source(vycmasma)' \
        interval=1000  \
        axis=freq,both device=fig2.ps/cps nxy=1,1 options=nobase,avall,jplcat,restfreq \
	vsource=30 veltype=lsr catpath=$MIRCAT/jplcat veldef=radio
$<
\rm -r vycmasma.SiS
uvaver vis=${out}_rx0.lsb.tsys.if1.bp.g.av select='source(vycmasma),window(4,5,6,7)' \
        out=vycmasma.SiS
uvflag vis=vycmasma.SiS \
        edge=9,9 \
        flagval=flag
$<
UVLINL
echo "9. L/C SEPARATION"
$<
\rm -r vycmasma.SiS.lin
uvlin vis=vycmasma.SiS chans=1,150,350,500 out=vycmasma.SiS.lin \
        order=1 mode=line options=nowin
#
#SiS at 235.9611 GHz  13=>12 
puthd in=vycmasma.SiS.lin/restfreq value=235.9611 type=double
smauvspec vis=vycmasma.SiS.lin hann=11 \
        select='source(vycmasma)' \
        interval=1000  \
        axis=vel,both device=/xs nxy=1,1 options=nobase,avall
#
#the source velocity is incorrect.
#need to identify the line and reset the rest frequency.
#
LMCUBE:
echo "9. CONSTRUCT LINE IMAGE CUBE" 
\rm -r vycmasma.SiS.map vycmasma.SiS.beam vycmasma.SiS.icmp
invert vis=vycmasma.SiS.lin map=vycmasma.SiS.map beam=vycmasma.SiS.beam \
        imsize=512,512 cell=0.15 robust=2 \
        options=systemp line=vel,60,-100,4
$<
clean map=vycmasma.SiS.map beam=vycmasma.SiS.beam out=vycmasma.SiS.icmp gain=0.08 \
        cutoff=0 niters=3500 \
        region=arcsec,'boxes(-5,-5,5,5)(1,60)'
$<
\rm -r vycmasma.SiS.icln
restor model=vycmasma.SiS.icmp beam=vycmasma.SiS.beam map=vycmasma.SiS.map  \
        out=vycmasma.SiS.icln
$<
imspect in=vycmasma.SiS.icln region=arcsec,'boxes(-0.,-0.,0.,0.)(1,60)' \
	device=fig3.ps/ps  
$<
MOMENT:
\rm -r vycmasma.SiS.flux
moment in==vycmasma.SiS.icln region='quarter(10,45)' out=vycmasma.SiS.flux mom=0 \
clip=-10,0.06
$<
echo "Display the restored clean image"
$<
cgdisp in=vycmasma.SiS.flux,vycmasma.SiS.flux type=c,p \
        region=arcsec,'boxes(-5,-5,5,5)' \
        xybin=1,1 device=fig4.ps/cps nxy=1,1 options=beambr,wedge,blacklab \
        labtyp=arcsec,arcsec range=5,40,lin,2 cols1=7 csize=1.2,1.2 \
        slev=p,1 levs1=7.5,15,25,35,45,55,65,75,85,95
echo "The end"
exit
