#!/bin/csh -f
#handling daul receivers non-polarization hybrid-resolution data 
#
#Responsible: Jun-Hui Zhao (SAO)
#Version 1.1
#jhz 2011-02-25 created for reduction of the USB data which 
#               contains the CO line in 340 receiver.
#jhz 2011-03-05 separated antenna-based tsys for different receivers.
#jhz 2011-03-10 tested the new methods for handling tsys.              
#jhz 2011-03-17 added CO identification with the JPL catalog.
#jhz 2011-03-25 added LSB data from both receivers for continuum.
#jhz 2011-04-01 using different color palette for illustration figure
#jhz 2011-04-25 updated the figures
#jhz 2011-05-06 updated the instruction lines
#jhz 2011-05-18 added goto tags, and tested smauvspec
#               correct typo in deleting file name
#
# testing data: 101007_06:42:16 from Ken (Taco) Young, PI: Charlie Qi
#               The data set has not been moved to the public domain yet.
#               Please get the PI's permission for use of the testing data
#               and the figures produced from this script.
# 
# demo illustration figures
# fig1.ps = continuum image of dmtau   
# fig2.ps = spectrum from the vis data of dmtau 
# fig3.ps = spectrum from CO(3-2) line image cube dmtau  
# fig4.ps = image of integrated CO(3-2) line flux dmtau  
#
set darea=/sma/SMAusers/jzhao/data/101007_06:42:16
set outf=101007
#goto continue
#
#########################################
# enable one of the next several goto 
# lines allow users to restart the script 
# from a working step.
#########################################
#
#
#goto SMALOD        #Load smadata
goto INSPECT       #Inspect visdata
#goto TSYSFIX       #Corrections for fringe amplitude
#goto BPCAL         #Bandpass calibration
#goto APPLYBP       #Apply bandpass solutions
#goto GCAL          #Gain calibration for target
#goto FLUX          #Determine flux density scale
#goto SPLIT         #Split the data into single-source file
#goto LCSEPARATE    #Line and continuum separation
#goto CONTMAP       #make continuum map
#goto LIDENTIFY     #molecular line identification
#goto LINECUBE      #make line image cube
#goto MOMENTMAP     #make moment map
#goto DIAPLAY       #display the result
#
SMALOD:
##############################################
echo "1. READ THE MIR DATA INTO MIRIAD"
##############################################
\rm -r ${outf}_rxl.*
smalod in=$darea out=$outf \
rxif=-1 sideband=2
\rm -r ${outf}_rxh.*
smalod in=$darea out=$outf \
rxif=-2 sideband=2
$<
INSPECT:
##############################################
echo "2. INSPECT THE DATA"
##############################################
uvindex vis=${outf}_rxl.usb
$<
uvindex vis=${outf}_rxh.usb
$<
uvindex vis=${outf}_rxl.lsb
$<
uvindex vis=${outf}_rxh.lsb
$<
smauvplt vis=${outf}_rxl.usb \
        axis=time,ampl device=/xs nxy=1,1
$<
smauvplt vis=${outf}_rxh.usb \
        axis=time,ampl device=/xs nxy=1,1
$<
smauvplt vis=${outf}_rxl.lsb \
        axis=time,ampl device=/xs nxy=1,1
$<
smauvplt vis=${outf}_rxh.lsb \
        axis=time,ampl device=/xs nxy=1,1
$<
smauvspec vis=${outf}_rxl.usb \
        select='source(dmtau)' \
        interval=1000 hann=1 \
        axis=chan,both device=/xs nxy=1,3
$<
smauvspec vis=${outf}_rxh.usb \
        select='source(dmtau)' \
        interval=1000 hann=1 \
        axis=chan,both device=/xs nxy=1,3

smauvspec vis=${outf}_rxl.lsb \
        select='source(dmtau)' \
        interval=1000 hann=1 \
        axis=chan,both device=/xs nxy=1,3
$<
smauvspec vis=${outf}_rxh.lsb \
        select='source(dmtau)' \
        interval=1000 hann=1 \
        axis=chan,both device=/xs nxy=1,3
$<
smavarplt vis=${outf}_rxl.usb device=/xs \
          xaxis=time \
          yaxis=systemp  nxy=2,4
$<
smavarplt vis=${outf}_rxh.usb device=/xs \
          xaxis=time \
          yaxis=systemp  nxy=2,4

$<
smavarplt vis=${outf}_rxl.lsb device=/xs \
          xaxis=time \
          yaxis=systemp  nxy=2,4
$<
smavarplt vis=${outf}_rxh.lsb device=/xs \
          xaxis=time \
          yaxis=systemp  nxy=2,4
$<
TSYSFIX:
###################################################
echo "3. CORRECTIONS FOR FRINGE AMPLITUDE"
###################################################
$<
\rm -r ${outf}_rxl.usb.tsys
smafix vis=${outf}_rxl.usb out=${outf}_rxl.usb.tsys \
        device=/xs xaxis=time \
        yaxis=systemp nxy=2,4   \
        options=tsyscorr 
$<
\rm -r ${outf}_rxh.usb.tsys
smafix vis=${outf}_rxh.usb out=${outf}_rxh.usb.tsys \
        device=/xs xaxis=time \
        yaxis=systemp nxy=2,4   \
        options=tsyscorr   

$<
\rm -r ${outf}_rxl.lsb.tsys
smafix vis=${outf}_rxl.lsb out=${outf}_rxl.lsb.tsys \
        device=/xs xaxis=time \
        yaxis=systemp nxy=2,4   \
        options=tsyscorr
$<
\rm -r ${outf}_rxh.lsb.tsys
smafix vis=${outf}_rxh.lsb out=${outf}_rxh.lsb.tsys \
        device=/xs xaxis=time \
        yaxis=systemp nxy=2,4   \
        options=tsyscorr
      
$<
smablflag vis=${outf}_rxl.usb.tsys \
        axis=time,ampl device=/xs 
$<
smablflag vis=${outf}_rxh.usb.tsys \
        axis=time,ampl device=/xs
$<
smablflag vis=${outf}_rxl.lsb.tsys \
        axis=time,ampl device=/xs
$<
smablflag vis=${outf}_rxh.lsb.tsys \
        axis=time,ampl device=/xs
$<
smauvspec vis=${outf}_rxl.usb.tsys \
        select='source(3c454.3)' \
        interval=1000  \
        axis=freq,both device=/xs nxy=1,5
$<
smauvspec vis=${outf}_rxh.usb.tsys \
        select='source(3c454.3)' \
        interval=1000  \
        axis=freq,both device=/xs nxy=1,5
$<
smauvspec vis=${outf}_rxl.lsb.tsys \
        select='source(3c454.3)' \
        interval=1000  \
        axis=freq,both device=/xs nxy=1,5
$<
smauvspec vis=${outf}_rxh.lsb.tsys \
        select='source(3c454.3)' \
        interval=1000  \
        axis=freq,both device=/xs nxy=1,5
$<
BPCAL:
#########################################
echo '4. BANDPASS CALIBRATION'
#########################################
smamfcal vis=${outf}_rxl.usb.tsys select='source(3c454.3)' edge=8,8  \
        weight=2 refant=3 interval=60000
$<
smagpplt vis=${outf}_rxl.usb.tsys  device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=1,8
smamfcal vis=${outf}_rxh.usb.tsys select='source(3c454.3)' edge=8,8  \
        weight=2 refant=3 interval=60000
$<
smagpplt vis=${outf}_rxh.usb.tsys  device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=1,8

smamfcal vis=${outf}_rxl.lsb.tsys select='source(3c454.3)' edge=8,8  \
        weight=2 refant=3 interval=60000
$<
smagpplt vis=${outf}_rxl.lsb.tsys  device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=1,8
smamfcal vis=${outf}_rxh.lsb.tsys select='source(3c454.3)' edge=8,8  \
        weight=2 refant=3 interval=60000
$<
smagpplt vis=${outf}_rxh.lsb.tsys  device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=1,8
$<
APPLYBP:
#####################################
# Apply bandpass calibration
#####################################
#
\rm -r ${outf}_rxl.usb.tsys.bp ${outf}_rxh.usb.tsys.bp
uvaver vis=${outf}_rxl.usb.tsys options=nocal \
       out=${outf}_rxl.usb.tsys.bp
$<
uvaver vis=${outf}_rxh.usb.tsys options=nocal \
       out=${outf}_rxh.usb.tsys.bp

\rm -r ${outf}_rxl.lsb.tsys.bp ${outf}_rxh.lsb.tsys.bp
uvaver vis=${outf}_rxl.lsb.tsys options=nocal \
       out=${outf}_rxl.lsb.tsys.bp
$<
uvaver vis=${outf}_rxh.lsb.tsys options=nocal \
       out=${outf}_rxh.lsb.tsys.bp
$<
GCAL:
#################################################
echo '5. GAIN CALIBRATION'
#################################################
$<
smamfcal vis=${outf}_rxl.usb.tsys.bp  \
select='source(0449+113)'  \
refant=3 interval=2 options=nopassol
$<
smagpplt vis=${outf}_rxl.usb.tsys.bp device=/xs \
        yaxis=amp,pha  \
        options=gains  polyfit=6 \
        nxy=2,4
smamfcal vis=${outf}_rxh.usb.tsys.bp  \
select='source(0449+113)'  \
refant=3 interval=2 options=nopassol
$<
smagpplt vis=${outf}_rxh.usb.tsys.bp device=/xs \
        yaxis=amp,pha  \
        options=gains  polyfit=6 \
        nxy=2,4

$<
smamfcal vis=${outf}_rxl.lsb.tsys.bp  \
select='source(0449+113)'  \
refant=3 interval=2 options=nopassol
$<
smagpplt vis=${outf}_rxl.lsb.tsys.bp device=/xs \
        yaxis=amp,pha  \
        options=gains  polyfit=6 \
        nxy=2,4
smamfcal vis=${outf}_rxh.lsb.tsys.bp  \
select='source(0449+113)'  \
refant=3 interval=2 options=nopassol
$<
smagpplt vis=${outf}_rxh.lsb.tsys.bp device=/xs \
        yaxis=amp,pha  \
        options=gains  polyfit=6 \
        nxy=2,4

\rm -r ${outf}_rxl.usb.tsys.bp.g ${outf}_rxh.usb.tsys.bp.g
uvaver vis=${outf}_rxl.usb.tsys.bp \
       out=${outf}_rxl.usb.tsys.bp.g
$<
uvaver vis=${outf}_rxh.usb.tsys.bp \
       out=${outf}_rxh.usb.tsys.bp.g

\rm -r ${outf}_rxl.lsb.tsys.bp.g ${outf}_rxh.lsb.tsys.bp.g
uvaver vis=${outf}_rxl.lsb.tsys.bp \
       out=${outf}_rxl.lsb.tsys.bp.g
$<
uvaver vis=${outf}_rxh.lsb.tsys.bp \
       out=${outf}_rxh.lsb.tsys.bp.g
$<
FLUX:
############################################
echo '6. BOOTSTRAPE THE FLUX DENSITY SCALE'
############################################
$<
smaflux vis=${outf}_rxl.usb.tsys.bp.g \
	select='source(callisto)' mirhome=$MIR
$<
smaflux vis=${outf}_rxh.usb.tsys.bp.g \
        select='source(callisto)' mirhome=$MIR

$<
smaflux vis=${outf}_rxl.lsb.tsys.bp.g \
        select='source(callisto)' mirhome=$MIR
$<
smaflux vis=${outf}_rxh.lsb.tsys.bp.g \
        select='source(callisto)' mirhome=$MIR
$<
\rm -r ${outf}_rxl.usb.tsys.bp.ga ${outf}_rxh.usb.tsys.bp.ga
uvaver vis=${outf}_rxl.usb.tsys.bp.g \
       out=${outf}_rxl.usb.tsys.bp.ga
$<
uvaver vis=${outf}_rxh.usb.tsys.bp.g \
       out=${outf}_rxh.usb.tsys.bp.ga
$<
\rm -r ${outf}_rxl.lsb.tsys.bp.ga ${outf}_rxh.lsb.tsys.bp.ga
uvaver vis=${outf}_rxl.lsb.tsys.bp.g \
       out=${outf}_rxl.lsb.tsys.bp.ga
$<
uvaver vis=${outf}_rxh.lsb.tsys.bp.g \
       out=${outf}_rxh.lsb.tsys.bp.ga
SPLIT:
###################################
# Split visibility on source based.  
###################################
$<
\rm -r dmtau.345233 dmtau.371859  dmtau.usb.350  dmtau.usb.400
uvsplit vis=${outf}_rxl.usb.tsys.bp.ga select='source(dmtau)' \
	options=nowindow
uvsplit vis=${outf}_rxh.usb.tsys.bp.ga select='source(dmtau)' \
        options=nowindow

\rm -r dmtau.337111 dmtau.363737  dmtau.lsb.350  dmtau.lsb.400
uvsplit vis=${outf}_rxl.lsb.tsys.bp.ga select='source(dmtau)' \
        options=nowindow
uvsplit vis=${outf}_rxh.lsb.tsys.bp.ga select='source(dmtau)' \
        options=nowindow
$<
mv dmtau.345233 dmtau.usb.350
mv dmtau.371859 dmtau.usb.400
mv dmtau.337111 dmtau.lsb.350
mv dmtau.363737 dmtau.lsb.400
$<
LCSEPARATE:
###################################
# Line and continuum separation  
###################################
$<
\rm -r dmtau.usb.350.ch0  dmtau.usb.400.ch0
uvlin vis=dmtau.usb.350 chans=1,3008 out=dmtau.usb.350.ch0 \
        order=1 mode=chan0 options=nowin
$<
uvlin vis=dmtau.usb.350 chans=1,3008 out=dmtau.usb.400.ch0 \
        order=1 mode=chan0 options=nowin
$<
\rm -r dmtau.lsb.350.ch0  dmtau.lsb.400.ch0
uvlin vis=dmtau.lsb.350 chans=1,3008 out=dmtau.lsb.350.ch0 \
        order=1 mode=chan0 options=nowin
$<
uvlin vis=dmtau.lsb.350 chans=1,3008 out=dmtau.lsb.400.ch0 \
        order=1 mode=chan0 options=nowin

CONTMAP:
#############################################
echo "7. CONTINUUM IMAGE"
#############################################
$<
\rm -r dmtau.map dmtau.beam dmtau.icmp
set visdata=dmtau.usb.400.ch0,dmtau.lsb.400.ch0,dmtau.usb.350.ch0,dmtau.lsb.350.ch0
invert vis=${visdata} map=dmtau.map beam=dmtau.beam \
        imsize=512,512 cell=0.3 robust=2 \
        options=systemp
$<
cgdisp in=dmtau.map type=pixel region=arcsec,'boxes(-15,-15,15,15)' \
        xybin=1,1 device=/xs nxy=1,1 options=full,beambr,wedge \
        labtyp=arcsec,arcsec range=0,0,lin,2 cols1=7 csize=1,1
$<
clean map=dmtau.map beam=dmtau.beam out=dmtau.icmp gain=0.08 \
        cutoff=0 niters=5000 \
        region=arcsec,'boxes(-1,-1,1,1)'
$<
\rm -r dmtau.icln
restor model=dmtau.icmp beam=dmtau.beam map=dmtau.map  \
        out=dmtau.icln
$<
cgdisp in=dmtau.icln type=pixel region=arcsec,'boxes(-15,-15,15,15)' \
        xybin=1,1 device=fig1.ps/cps nxy=1,1 options=beambr,wedge,blacklab \
        labtyp=arcsec,arcsec range=0.005,0.14,lin,2 cols1=7 csize=1,1
$<
imstat in=dmtau.icln region=arcsec,'boxes(-15,-15,-12,-12)' \
         axes=ra,dec
$<
\rm -r dmtau.345800
uvsplit vis=dmtau.usb.350 select='window(8)' 
$<
mv dmtau.345800 dmtau.co
$<
LIDENTIFY:
####################################################
echo "8. LINE IDENTIFICATION"
####################################################
$<
#...
#skip
#select  27002    HNC           (<type> y to select it or n to skip it or t to terminate)
#n
#skip
#select  27003    HCN-v2        (<type> y to select it or n to skip it or t to terminate)
#n
#skip
#select  27004    C-13-N        (<type> y to select it or n to skip it or t to terminate)
#n
#skip
#select  28001    CO            (<type> y to select it or n to skip it or t to terminate)
#y
#yes
#select  28002    HC-13-N       (<type> y to select it or n to skip it or t to terminate)
#t
#terminate selection
smauvspec vis=dmtau.usb.350 hann=15 select='source(dmtau)' \
        interval=1000  options=nobase,avall,jplcat,restfreq \
        axis=freq,both device=fig2.ps/cps nxy=1,1 \
        vsource=6.5 veltype=lsr catpath=$MIRCAT/jplcat veldef=radio
$<
\rm -r dmtau.co.lin
uvlin vis=dmtau.co chans=1,100,170,250 out=dmtau.co.lin \
        order=1 mode=line options=nowin
$<
smauvspec vis=dmtau.co.lin hann=7 \
        interval=1000  options=avall,nobase \
        axis=freq,both device=/xs nxy=1,1
$<
LINECUBE:
####################################################
echo '9. LINE IMAGING'
####################################################
$<
\rm -r dmtau.l.map dmtau.l.beam dmtau.l.icmp
invert vis=dmtau.co.lin map=dmtau.l.map beam=dmtau.l.beam \
        imsize=1024,1024 cell=0.2 robust=0 \
        options=systemp line=vel,100,-10,0.35 
$<
clean map=dmtau.l.map beam=dmtau.l.beam out=dmtau.l.icmp gain=0.08 \
        cutoff=0.05 niters=10000 \
        region=arcsec,'boxes(-20,-20,20,20)(1,100)'
$<
\rm -r dmtau.l.icln
restor model=dmtau.l.icmp beam=dmtau.l.beam map=dmtau.l.map  \
        out=dmtau.l.icln
$<
cgdisp in=dmtau.l.icln type=pixel \
	region=arcsec,'boxes(-15,-15,15,15)(1,100)' \
        xybin=1,1 device=/xs nxy=5,5 options=beambr,wedge,blacklab \
        labtyp=arcsec,arcsec range=0,0,lin,2 cols1=7 csize=1,1
$<
imstat in=dmtau.l.icln region=arcsec,'boxes(-25,-25,-22,-22)' \
         axes=ra,dec

$<
MOMENTMAP:
#####################################
echo "spectrum of the CO line"
#####################################
$<
imspect in=dmtau.l.icln region=arcsec,'boxes(-0.,-0.,0.,0.)(1,100)' \
	device=fig3.ps/ps yrange=-0.5,2.
$<
\rm -r dmtau.l.flux
moment in==dmtau.l.icln region='quarter(20,70)' out=dmtau.l.flux mom=0 \
clip=-10,0.25
$<
DISPLAY:
######################################
echo "10. DISPLAY THE LINE FLUX IMAGE"
######################################
$<
cgdisp in=dmtau.l.flux,dmtau.l.flux type=c,p \
        region=arcsec,'boxes(-8,-8.5,8,7.5)' \
        xybin=1,1 device=fig4.ps/cps nxy=1,1 options=beambr,wedge,blacklab \
        labtyp=arcsec,arcsec range=0.01,4,lin,8 cols1=3 csize=1.2,1.2 \
        slev=p,1 levs1=6,10,15,25,35,45,55,65,75,85,95	
exit