#!/bin/csh -f
# handling SMA high-resolution (HiRes) spectral data
#
# Responsible: Jun-Hui Zhao (SAO)
# Version 2.0
# jhz 2006-05-1
# jhz 2011-02-17 update
# jhz 2011-02-28 update
# jhz 2011-05-06 update by organizing the instruction lines
# jhz 2011-06-09 changed the data path to CF where the test
#                data located.
#
# testing data: 060512_18:10:36 from Ken (Taco) Young 
#
# demo illustration figures
# fig1.ps = continuum image of R-Cas 
# fig2.ps = spectrum from the vis data of R-Cas
# fig3.ps = spectrum from line image cube R-Cas
# fig4.ps = image of integrated line flux R-Cas 
goto continue
#
#continue:
set darea=/sma/SMAusers/jzhao/data/060512_18:10:36
#########################################
echo "1. READ THE MIR DATA INTO MIRIAD"
#########################################
\rm -r  060512_rx0.lsb  060512_rx0.usb 
smalod in=$darea out=060512 rxif=0 \
       sideband=0 readant=8  nscans=1,285
$<
#########################################
echo "2. INSPECT THE RAW UV DATA"
echo "List the observing scans in the uv data"
#########################################
uvindex vis=060512_rx0.lsb
$<
#########################################
echo "Plot the visibilities of the raw data"
#########################################
smauvplt vis=060512_rx0.lsb  \
	axis=time,ampl device=/xs nxy=1,1 
$<
smauvspec vis=060512_rx0.lsb \
        select='source(R-Cas)' \
        interval=1000 hann=1 \
        axis=chan,both device=/xs nxy=1,3
$<
uvflag vis=060512_rx0.lsb select='window(15)' flagval=flag
uvflag vis=060512_rx0.lsb select='window(19)' flagval=flag 
#########################################
echo "2. CORRECTIONS FOR AMPLITUDE"
echo "Tsys correction"
#########################################
\rm -r 060512_rx0.lsb.tsys
smafix vis=060512_rx0.lsb out=060512_rx0.lsb.tsys \
        device=/xs xaxis=time   \
	yaxis=systemp nxy=2,4  options=tsyscorr 
#########################################
echo "Inspect vis data and flagging bad data"
#########################################
smablflag vis=060512_rx0.lsb axis=time,ampl device=/xs
#########################################
echo "Inspect spectral data again"
echo "for R-Cas"
#########################################
smauvspec vis=060512_rx0.lsb.tsys \
        select='source(R-Cas)' \
        interval=1000 hann=7  \
        axis=chan,both device=/xs nxy=1,3
#########################################
echo "neptune"
#########################################
smauvspec vis=060512_rx0.lsb.tsys \
	select='source(nep*)' \
        interval=1000  \
        axis=chan,both device=/xs nxy=1,3
#########################################
echo "ceres"
#########################################
smauvspec vis=060512_rx0.lsb.tsys \
        select='source(cer*)' \
        interval=1000  \
        axis=chan,both device=/xs nxy=1,3
##########################################
echo "4. BANDPASS CALIBRATIONS"
##########################################
$<
##########################################
echo "4.1. Correction for Possible Bandpass Ripples "
\rm -r 060512_rx0.lsb.tsys.lr 060512_rx0.lsb.tsys.hr  
##########################################
echo "separate the first 12 chunks to be used for continuum gain calibration"
##########################################
uvaver vis=060512_rx0.lsb.tsys select='window(1,2,3,4,5,6,7,8,9,10,11,12)' \
        options=nocal out=060512_rx0.lsb.tsys.lr
##########################################
echo "separate the 16th chunk for the high resolution"
##########################################
uvaver vis=060512_rx0.lsb.tsys select='window(16)' \
        options=nocal out=060512_rx0.lsb.tsys.hr
##########################################
echo "check for low resolution data"
##########################################
smauvspec vis=060512_rx0.lsb.tsys.lr \
        select='source(nept*)' \
        interval=1000  \
        axis=chan,both device=/xs nxy=1,3
$<
##########################################
echo "check for high resolution data"
##########################################
smauvspec vis=060512_rx0.lsb.tsys.hr \
        select='source(nept*)' \
        interval=1000  \
        axis=chan,both device=/xs nxy=1,3
##########################################
echo "Solve for the antenna-based bandpass"
##########################################
smamfcal vis=060512_rx0.lsb.tsys.lr select='source(neptune)' edge=3,3   \
         refant=3 interval=600
smamfcal vis=060512_rx0.lsb.tsys.hr select='source(neptune)'    \
         weight=2 refant=3 interval=600 smooth=7,1,0.9 options=msmooth
#########################################
echo "Check up the bandpass solutions"
echo "low resolution"
#########################################
smagpplt vis=060512_rx0.lsb.tsys.lr  device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=2,4
$<
#########################################
echo "high resolution ampl"
#########################################
smagpplt vis=060512_rx0.lsb.tsys.hr  device=/xs \
        yaxis=amp \
        options=bandpass  polyfit=5 \
        nxy=2,4 yrange=0.5,1.5
$<
#########################################
echo "high resolution phase"
#########################################
smagpplt vis=060512_rx0.lsb.tsys.hr  device=/xs \
        yaxis=phase \
        options=bandpass  polyfit=5 \
        nxy=2,4 yrange=-50,50
$<
\rm -r 060512_rx0.lsb.tsys.lr.bp 060512_rx0.lsb.tsys.hr.bp
#########################################
echo "apply the bandpass"
#########################################
uvaver  vis=060512_rx0.lsb.tsys.lr options=nocal \
        out=060512_rx0.lsb.tsys.lr.bp
uvaver  vis=060512_rx0.lsb.tsys.hr options=nocal \
        out=060512_rx0.lsb.tsys.hr.bp
#########################################
echo "Upto this point, the calibration is done for antenna-based."
echo "4.2. Corrections for Baseline-based Errors in Bandpass"
#########################################
smauvspec vis=060512_rx0.lsb.tsys.lr.bp \
        select='source(nep*)' \
        interval=1000   \
        axis=freq,both device=/xs nxy=2,2 
#########################################
echo "looks fine; skip baseline-based calibration"
echo "flag possible band vis  points"
#########################################
smablflag vis=060512_rx0.lsb.tsys.lr.bp \
        axis=time,ampl device=/xs 
#########################################
echo "5. ANTENNA GAIN CALIBRATION"
#########################################
echo "Up to this point, we have corrected for bandpass errors."
echo "For 0133+476"
#
\rm -r  0133+476.216816
uvsplit vis=060512_rx0.lsb.tsys.lr.bp     \
	select='source(0133+476)' options=nowindow
\rm -r 0133+476.lsb.tsys.lr.bp
cp -r 0133+476.216816 0133+476.lsb.tsys.lr.bp
selfcal vis=0133+476.lsb.tsys.lr.bp \
refant=5 interval=.5 options=pha
#########################################
echo "Plot the gain solutions"
#########################################
smagpplt vis=0133+476.lsb.tsys.lr.bp device=/xs \
        yaxis=amp,pha  \
        options=gains  polyfit=6 \
        nxy=2,4
$<
gpcopy vis=0133+476.lsb.tsys.lr.bp out=060512_rx0.lsb.tsys.lr.bp \
	mode=copy options=nopass
gpcopy vis=0133+476.lsb.tsys.lr.bp out=060512_rx0.lsb.tsys.hr.bp \
        mode=copy options=nopass
puthd in=060512_rx0.lsb.tsys.lr.bp/interval value=600 type=double
puthd in=060512_rx0.lsb.tsys.hr.bp/interval value=600 type=double
###########################################
echo "6. BOOTSTRAPE THE FLUX DENSITY SCALE"
###########################################
smaflux vis=060512_rx0.lsb.tsys.lr.bp \
	select='source(ceres)' mirhome=$MIR
gpcopy vis=060512_rx0.lsb.tsys.lr.bp \
        out=060512_rx0.lsb.tsys.hr.bp mode=copy \
        options=nopass
###########################################
echo "Apply the gains to the data"
###########################################
\rm -r 060512_rx0.lsb.tsys.hr.bp.av 060512_rx0.lsb.tsys.lr.bp.av 
uvaver vis=060512_rx0.lsb.tsys.hr.bp out=060512_rx0.lsb.tsys.hr.bp.av
uvaver vis=060512_rx0.lsb.tsys.lr.bp out=060512_rx0.lsb.tsys.lr.bp.av
$<
##########################################
echo "7. Continuum IMAGE"
##########################################
uvflag vis=060512_rx0.lsb.tsys.lr.bp.av edge=7,7 flagval=flag
smauvspec vis=060512_rx0.lsb.tsys.lr.bp.av \
        select='source(R*)' hann=9 \
        interval=1000   \
        axis=chan,both device=/xs nxy=1,1 options=nobase,avall
$<
\rm -r R-Cas.map R-Cas.beam R-Cas.icmp
invert vis=060512_rx0.lsb.tsys.lr.bp.av  \
        map=R-Cas.map beam=R-Cas.beam \
        imsize=128,128 cell=0.3 sup=0 \
        select='source(R-Cas)' \
        line=channel,1,1,1536 options=systemp,double
$<
clean map=R-Cas.map beam=R-Cas.beam out=R-Cas.icmp gain=0.1 \
        cutoff=0.01 niters=1000 \
        region=arcsec,'boxes(-10,-10,10,10)'
$<
selfcal vis=060512_rx0.lsb.tsys.lr.bp.av \
        select='source(R*)' model=R-Cas.icmp \
        refant=5 interval=.5 options=pha
$<
\rm -r R-Cas.map R-Cas.beam R-Cas.icmp
invert vis=060512_rx0.lsb.tsys.lr.bp.av  \
        map=R-Cas.map beam=R-Cas.beam \
        imsize=256,256 cell=0.3 sup=0 \
        select='source(R-Cas)' \
        line=channel,1,1,1536 options=systemp,double
$<
clean map=R-Cas.map beam=R-Cas.beam out=R-Cas.icmp gain=0.1 \
        cutoff=0.01 niters=1000 \
        region=arcsec,'boxes(-10,-10,10,10)'
echo "Restore cleaned image"
\rm -r R-Cas.icln
restor model=R-Cas.icmp beam=R-Cas.beam map=R-Cas.map  \
        out=R-Cas.icln
echo "Display the restored clean image"
cgdisp in=R-Cas.icln type=pixel \
        region=arcsec,'boxes(-25,-25,25,25)' \
        xybin=1,1 device=fig1.ps/cps nxy=1,1 options=beambr,wedge,blacklab \
        labtyp=hms,dms range=0,0,lin,10 cols1=7 csize=1.2,1.2
$<
###########################################
echo "8. LINE IMAGE CUBE"
###########################################
echo "R-Cas"
###########################################
uvflag vis=060512_rx0.lsb.tsys.hr.bp.av edge=100,100 flagval=flag
smauvspec vis=060512_rx0.lsb.tsys.hr.bp.av device=fig2.ps/cps hann=9 \
        axis=chan,both select='source(R-Cas)' interval=1000 nxy=1,1 \
        options=nobase,avall
$<
###########################################
echo "Make a dirty image"
###########################################
\rm -r R-CasLin.map R-CasLin.beam R-CasLin.icmp
invert vis=060512_rx0.lsb.tsys.hr.bp.av  \
	map=R-CasLin.map beam=R-CasLin.beam \
        imsize=256,256 cell=0.3 sup=0 \
	select='source(R-Cas)' \
        line=channel,400,900,1 options=systemp,double
###########################################
echo "Clean  the dirty image"
###########################################
clean map=R-CasLin.map beam=R-CasLin.beam out=R-CasLin.icmp gain=0.1 \
        cutoff=0.01 niters=1000 \
        region=arcsec,'boxes(-10,-10,10,10)'
###########################################
echo "Restore cleaned image"
###########################################
\rm -r R-CasLin.icln
restor model=R-CasLin.icmp beam=R-CasLin.beam map=R-CasLin.map  \
        out=R-CasLin.icln
$<
###########################################
echo "Display the restored clean image"
###########################################
cgdisp in=R-CasLin.icln type=pixel \
	region=arcsec,'boxes(-15,-15,15,15)(1,400)' \
        xybin=1,1 device=/xs nxy=10,10 options=full,beambr,wedge \
        labtyp=arcsec,arcsec range=0,0,lin,10 cols1=7 csize=1.2,1.2 
$<
###########################################
echo "spectrum of the central 2 arcsec"
###########################################
imspect in=R-CasLin.icln region=arcsec,'boxes(-1,-1,1,1)(1,400)' \
	device=fig3.ps/ps yrange=-1,12  
$<
\rm -r R-CasLin.flux
moment in=R-CasLin.icln region='quarter(101,300)' out=R-CasLin.flux
$<
##########################################
echo "Display the restored clean image"
##########################################
cgdisp in=R-CasLin.flux,R-CasLin.flux type=c,p \
        region=arcsec,'boxes(-10,-10,10,10)' \
        xybin=1,1 device=fig4.ps/ps nxy=1,1 options=beambr,wedge,blacklab \
        labtyp=arcsec,arcsec range=0,0,lin,2 cols1=7 csize=1.2,1.2 \
        slev=p,1 levs1=10,15,25,35,45,55,65,75,85,95	
exit
Updated 2010
Updated 2011