#! /bin/csh -f
# handling standard 2 GHz polarization data 
# 
# Responsible: Josep Miquel Girart (CSIC-IEEC) & Jun-Hui Zhao (SAO)
# Version 1.2 
# testing data: 070728_02:39:09 from SMA archive
#
# demo illustration figures
# fig1.ps = image of stokes I, total intensity
# fig2.ps = image of stokes Q 
# fig3.ps = image of stokes U
# fig4.ps = image of polarization vector (linear lines)
#           overlaid on polarization intensity (greyscaled)
#
##################################################################
echo ' '
echo '************************************************************'
echo 'Acknowledgement:'
echo 'The reduction procedure for SMA polarization data written in'
echo 'the following script was based on the communications with'  
echo 'Josep Miquel Girart, Meredith Hughes, Dan Marrone, Ram Rao,'
echo 'Bob Sault and Mel Wright.'
echo ' '
echo 'For further information on the SMA polarimeter'
echo 'Please read the article written by Dan Marrone and Ram Rao'
echo 'http://adsabs.harvard.edu/abs/2008SPIE.7020E..60M'
echo 'and the Ph.D. Thesis of Dan Marrone.'
echo '************************************************************'
echo ' '
################################################################## 
#
#jhz 2011-04-04 start the script (Ver1.0) for reduction of polarization data
#jhz 2011-04-05 check the polarization sampling pattern
#               check the time steps on the sources
#               check the e-vector and mount type 
#               check the parallactic angle coverage of
#               pol calibarators
#jhz 2011-04-07 check gain calibrators
#               check the algorithm in solving for D-terms
#               and the error assessments
#jhz 2011-04-08 worked on the pol data editing tool programs
#               fixed a few bugs
#jhz 2011-04-11 test the calibration with the data from the GC pol
#               program. 
#               .......
#jhz 2011-04-20 contact SMA users for scripts of polarization data
#               reductions. Get positive responses.
#               .......
#
#jhz 2011-05-04 receive a pol-data reduction script from Josep Miquel Girart
#
#jhz 2011-05-04 merge the two scripts to Ver 1.1
#               change the number of argv to 2
#               make an e2e script
#               check variable setups
#               comment out puthd for setup of QWP evector and type of mount
#               comment out uvredo for recalculate parallatic angles
#               which have been coded in both Miriad-CVS and SMA beta Miriad
#               since MIRIAD CVS 4.0.3
#
#jhz 2011-05-05 correct fcal input
#               correct pcal input 
#               correct fname setup
#               define variable dt
#               comment out foreach loop which was failed
#               split the individual sources and calibrators separately
#               change variable name ant to refant
#               add uvflag to flag pol(i)
#               add variable prcal for reference point pol source
#               remove interval=5 in the uvaver in splitting the sources 
#               before gain calibration
#               change to line=ch,24,7,114,128
#               add the variable avetime.
#
#jhz 2011-05-06 reset the interval 30 to 3000 in smamfcal in 
#               solving for bandpass 
#               add selfcal to calibrate the gains before solving for
#               D-terms.
#               add uvflux to calculate the flux density
#               of pol-calibrator that is used in GPCAL.
#               test the polarization calibration with a reference 
#               point source.
#               receive a part of polarization imaging.
#
#jhz 2011-05-09 add the instruction lines for polimage
#               work out imaging for four stokes.
#               add the statistics for stokes images
#               derive the rms for i and sqrt(q*u) which
#               are used in impol.
#               work out imaging for pol intensity
#                                and pol vector. 
#               create four postscript figures for illustration.
#               make Ver1.2
#
#jhz 2011-05-10 send to pol experts/users for review and comments
#               make an e2e test on jupiter
#               complete mfclean setup
#               add an if-then in task tags
#
#jhz 2011-05-11 re-organize the instruction lines
#               submit to sma beta Miriad site
#
#jhz 2011-05-12 change interval to =${avetime} in POLCAL section
#               update the script
#
#jhz 2011-05-13 update the instruction lines 
#
#jhz 2011-06-09 changed the data path to CF where the test data located
#
#jhz 2011-06-13 Following the suggestions from software meeting today,
#               commented out blflag and add a note for instruction. 
$<
echo 'Press return to continue'
$<
###########################################################
#enable the following line: 'goto SMATEST'
#allow users to select their own visibility data
#for the reduction of polarization data with
#this script.
###########################################################
goto SMATEST
if($#argv<2) then
echo "Incorrect usage: Needs Task Date Time"
echo "example: sma_pol_cal.csh 070728 02:39:09"
goto terminate
endif
##########################################
#--> Define date, time ofthe observations
##########################################
set date=$1
set hour=$2
goto USERSREDUCTION
SMATEST:
##########################################
# for SMA software test
##########################################
set date=070728
set hour=02:39:09
USERSREDUCTION:
##########################################
#--> Define side band (calibration is done 
#    separately for the USB and LSB)
#    nn=0 for LSB, nn=1 for USB
##########################################
set sb     = usb 
set nn    = 1
set rx    = rx1
set fname = "$date""_""$rx".$sb
echo $fname
##########################################
#-->  Define the Flux, Gain, Bandpass and 
#     Polarization Calibrators used
##########################################
set fcal  = uranus
set gcal  = 1749+096 
set bcal  = 3c454.3  
set pcal  = 3c454.3
set prcal = 3c273
#########################################
#-->  Define the Source name
#########################################
set font  = g3141
set source = ''
#########################################
#--> Define reference antenna
#           average time interval
#########################################
set refant  = 2
set avetime = 5
#########################################
#--> Other definitions
set region  = quart
set regdisp = 'arcsec,box(-12,-12,12,12)'
set dt      = flx
#
#########################################
#goto $1
#goto continue
##################################################################
# For more detailed and complete information on standard SMA 
# calibration, 
# please check: 
# http://www.cfa.harvard.edu/sma/miriad/miriad4sma/stepbystep.html 
################################################################## 
#
echo;echo
#
#echo "=====================================================================" 
#echo "     Script for SMA polarization calibration (noted by Girart)       "
#echo "=====================================================================" 
#echo " script based on the pol cal done of NGC 1333 IRAS 4A & G31.41+0.31: "
#echo " Girart, Rao, Marrone 2006, Science, 313, 812                        "
#echo " Girart, et al. 2009, Science, 324, 1408                             "
#echo "----------------------> version 03may11 <----------------------------"
#
#           Warning in the original script from Girart 
#
echo "WARNING: This is a simplified script.  It is recommended to check all"
echo "different calibration steps using the appropiate  visualization tools"
echo "for the visibility data. "
echo "=====================================================================" 
echo;echo
#
#########################################
# enable one of the next several goto 
# lines allow users to restart the script 
# from a working step.
#########################################
#
#goto SMALOD        #Load SMA data
#goto INSPECT       #Inspect visibility data
#goto IPOLFLAG      #Flag data with incorrect pol tags
#goto TSYSFIX       #Corrections for fringe amplitude
#goto FLAG          #Flag bad visibilities
#goto BPCAL         #Bandpass calibration
#goto APPLYBP       #Apply bandpass solutions
#goto FLUX          #Determine flux density scale
#goto SPLIT         #Split the data into single-source file
#goto GETFLUX       #Get the flux density of the pol-calibrator
#goto POLCAL        #Do polarization calibration
#goto GCAL          #Gain calibration for target
#goto AVER          #Make averaging visibility to a longer time interval
#goto POLMAP        #Start polarization mapping
#goto DIRTYMAP      #Create dirty images for stokes
#goto CLEANMAP      #Clean the stokes' images
#goto RESTOR        #Restore the stokes' images from clean-pixel images
#goto DISPLAYSTOKES #Display the stokes' images
#goto STATISTICS    #Do statistics for the stokes'images
#goto POLVEC        #Construct polarization intensity and vector images 
#
SMALOD:
#####################################
# Convert raw data into miriad format
#####################################
echo $fname
\rm -fr $fname
smalod in="/sma/SMAusers/jzhao/data/"$date"_"$hour out=$date rxif=1 \
       options=circular sideband=$nn refant=${refant} nscans=5,
#goto terminate
$<
echo 'Press return to continue'
$<
INSPECT:
#####################################
# Inspect at 
#####################################
uvindex  vis=$fname
$<
smacheck vis=$fname var=systemp flagval=flag range=0,1000
$<
smauvplt vis=$fname axis=time,ampl device=/xs nxy=1,1
$<
IPOLFLAG:
#####################################
# Flag data with incorrect pol tags
#####################################
uvflag vis=$fname select='pol(i)' flagval=flag
\rm -fr    $fname.flgd
uvflag vis=$fname select='-pol(ll,lr,rl,rr)' flagval=f
uvcal  vis=$fname options=unflagged out=$fname.flgd
#goto terminate

TSYSFIX:
#####################################
# Tsystem correction:
#####################################
\rm -fr    $fname.tsys
smafix vis=$fname.flgd out=$fname.tsys device=/xw nxy=1,1 \
       yrange=100,800 options=tsyscorr,dosour
#goto terminate

FLAG:
#####################################
# Flag corrupted data: check your 
# visibility data to look for bad data
#####################################
uvflag vis=$fname.tsys flagval=f edge=7 
$<
uvflag vis=$fname.tsys flagval=f select='time(13:59,14:11),sou('$bcal')'
$<
#
# users can enable blflag for flagging the bad visibilities.
#
#blflag vis=$fname.tsys axis=time,ampl device=/xs select='sou('$bcal')'
#
#goto terminate

BPCAL:
#####################################
# Bandpass calibration
#####################################
set options=averrll
smamfcal select='pol(ll,rr),sou('$bcal')' refant=$refant \
	interval=3000 weight=2 \
        vis=$fname.tsys  options=$options edge=7 
#goto terminate

APPLYBP:
#####################################
# Apply bandpass calibration
#####################################
\rm -fr    $fname.tsys.bp
uvaver vis=$fname.tsys out=$fname.tsys.bp options=nocal
#goto terminate

#################################################################
#
# Define the correct quarter wave plate paramenters.
# these parameters have been coded in both Miriad-CVS 
# and SMA beta Miriad
# since Miriad-CVS 4.0.3
#
#EVEC: 
#puthd in=$fname.tsys.bp/evector value=0.7853
# which has been coded in Miriad
#puthd in=$fname.tsys.bp/mount value=4
# which has been coded in Miriad
#\rm -fr    $fname.tsys.bp.chi
#uvredo vis=$fname.tsys.bp options=chi out=$fname.tsys.bp.chi
# which has been coded in Miriad
#
#goto terminate
#################################################################


FLUX:
#################################################################
# Flux Calibration
#################################################################

echo $fcal
smaflux vis=$fname.tsys.bp select='pol(ll,rr),so('$fcal')' mirhome=$MIR
#goto terminate
$<
echo 'Press return to continue'
$<

SPLIT:
#################################################################
# Split visibility on source based.  
#################################################################
#goto terminate

set so=$fcal
\rm -fr "$so""_""$dt".$sb
uvaver vis=$fname.tsys.bp select='sou('$so')' \
       out="$so""_""$dt".$sb line=ch,24,7,114,128
set so=$gcal
\rm -fr "$so""_""$dt".$sb
uvaver vis=$fname.tsys.bp select='sou('$so')' \
       out="$so""_""$dt".$sb line=ch,24,7,114,128
set so=$bcal
\rm -fr "$so""_""$dt".$sb
uvaver vis=$fname.tsys.bp select='sou('$so')' \
       out="$so""_""$dt".$sb line=ch,24,7,114,128
set so=$pcal
\rm -fr p"$so""_""$dt".$sb
uvaver vis=$fname.tsys.bp select='sou('$so')' \
       out=p"$so""_""$dt".$sb line=ch,24,7,114,128
set so=$prcal
\rm -fr "$so""_""$dt".$sb
uvaver vis=$fname.tsys.bp select='sou('$so')' \
       out="$so""_""$dt".$sb line=ch,24,7,114,128
set so=$font
\rm -fr "$so""_""$dt".$sb
uvaver vis=$fname.tsys.bp select='sou('$so')' \
       out="$so""_""$dt".$sb line=ch,24,7,114,128
#goto terminate

GETFLUX:
###################################################
# get the flux density for the pol calibrator
###################################################
uvflux vis=p$pcal"_"$dt.$sb stokes=i
###################################################
# Remove the phase error
###################################################
selfcal vis=p$pcal"_"$dt.$sb refant=${refant} \
	options=pha select='pol(ll,rr)' \
	interval=${avetime} 
###################################################
# Averaging visibility to a long time interval 
###################################################
\rm -fr pp$pcal"_"$dt.$sb
uvaver vis=p$pcal"_"$dt.$sb out=pp$pcal"_"$dt.$sb \
	interval=${avetime}

POLCAL:
###################################################
# Polarization calibration
###################################################
$<
# users can enable blflag for flagging the bad visibilities.
#
#blflag vis=pp$pcal"_"$dt.$sb axis=real,imag options=nobase \
#	device=/xs select='sou('$bcal')'
$<
gpcal vis=pp$pcal"_"$dt.$sb refant=${refant} \
	options=noxy,circular,qusolve \
	interval=${avetime} flux=14.47
#goto terminate

GCAL:
###################################################
# Gain Calibation for target. 
# "smamfcal options=nopassol" can also be used 
# instead of selfcal
###################################################
$<
selfcal vis=$gcal"_"$dt.$sb refant=${refant} \
	options=amp select='pol(ll,rr)' \
	interval=${avetime}
$<
smagpplt vis=$gcal"_"$dt.$sb device=/xw yaxis=phas \
	nxy=2,4
#goto terminate

AVER:
###################################################
# Averaging & Applying Gain & Pol calibration to the target
###################################################
gpcopy  vis=$gcal"_"$dt.$sb out=$font"_"$dt.$sb mode=copy
$<
\rm -rf $font"_"$dt.$sb.g
uvaver  vis=$font"_"$dt.$sb out=$font"_"$dt.$sb.g interval=${avetime}
$<
gpcopy  vis=pp$pcal"_"$dt.$sb out=$font"_"$dt.$sb.g options=nocal
#goto terminate

\rm -fr    p$font"_"$dt.$sb.g
uvaver vis=$font"_"$dt.$sb.g  out=p$font"_"$dt.$sb.g
$<
#goto terminate

###############################################################
# Gain calibration for a point source and apply pol calibration 
###############################################################
selfcal vis=$prcal"_"$dt.$sb refant=${refant} options=pha \
	select='pol(ll,rr)' interval=${avetime}
$<
\rm -fr $prcal"_"$dt.$sb.g
uvaver vis=$prcal"_"$dt.$sb out=$prcal"_"$dt.$sb.g \
	interval=${avetime}
gpcopy  vis=pp$pcal"_"$dt.$sb out=$prcal"_"$dt.$sb.g \
	options=nocal
$<
\rm -fr p$prcal"_"$dt.$sb.g
uvaver vis=$prcal"_"$dt.$sb.g out=p$prcal"_"$dt.$sb.g \
	interval=${avetime}
#goto terminate

POLMAP:
##################################################
echo ' Polarization Imaging'
##################################################
# select source for polarization mapping
# in next a few set up lines.
##################################################
$<
#########################
# select source
#########################
set source=$pcal

#########################
# select pol visdata set 
#########################
set polvis=pp$pcal"_"$dt.$sb

#########################
# setup the imaging size
#########################
set imsize=512

#########################
# setup the cell size
#########################
set cell=0.1

#########################
# setup clean region
#########################
set clreg='arcsec,box(-3,-3,3,3)'

#########################
# setup display region
#########################
set disreg='arcsec,box(-8,-8,8,8)'

#########################
# setup pol display region
#########################
set pdisreg='arcsec,box(-4,-4,4,4)'

#########################
# setup the area for statistics
#########################
set statreg='arcsec,box(-10,-10,-8,-8)'
#goto terminate


echo 'POLARIZATION IMAGING'
$<
echo 'Press return to continue'
$<

DIRTYMAP:
########################################
#make sure to have a right vis file
#corresponding to the polarization source 
#for imaging ->  invert vis=????????????
########################################
if ($source == '') then
goto POLMAP
endif
\rm -fr $source.$sb.i $source.$sb.q $source.$sb.u $source.$sb.v $source.$sb.b
#
#######################################
echo ' Make dirty images for four stokes'
#######################################
$<
invert vis=$polvis beam=$source.$sb.b \
	imsize=$imsize cell=$cell \
	map=$source.$sb.i,$source.$sb.q,$source.$sb.u,$source.$sb.v \
	options=mfs,sdb,systemp sup=0 \
	stokes=i,q,u,v robust=2 line=chan,24,1,1
$<
cgdisp in=$source.$sb.i,$source.$sb.i type=p,c \
	device=dirty.map/cps \
	labtyp=arcsec nxy=1,1 slev=p,10 \
	levs1=2,3,4,5,6,7,8,9 \
	options=wedge,mirror region=$disreg
#goto terminate

CLEANMAP:
#######################################
echo ' Clean dirty images for four stokes'
#######################################
$<
if ($source == '') then
goto POLMAP
endif
\rm -fr $source.$sb.icln $source.$sb.qcln $source.$sb.ucln $source.$sb.vcln 
\rm -fr $source.$sb.icmp $source.$sb.qcmp $source.$sb.ucmp $source.$sb.vcmp
mfclean map=$source.$sb.i beam=$source.$sb.b out=$source.$sb.icmp \
	region=$clreg gain=0.05 niters=150
$<
mfclean map=$source.$sb.q beam=$source.$sb.b out=$source.$sb.qcmp \
        region=$clreg gain=0.05 niters=50
$<
mfclean map=$source.$sb.u beam=$source.$sb.b out=$source.$sb.ucmp \
        region=$clreg gain=0.05 niters=50
$<
mfclean map=$source.$sb.v beam=$source.$sb.b out=$source.$sb.vcmp \
        region=$clreg gain=0.05 niters=50
#goto terminate

RESTOR:
#######################################
echo ' Restore cleaned pixel images for' 
echo ' four the stokes i, q, u, and v  '
#######################################
$<
if ($source == '') then
goto POLMAP
endif
$<
restor map=$source.$sb.i beam=$source.$sb.b out=$source.$sb.icln \
	model=$source.$sb.icmp
$<
restor map=$source.$sb.q beam=$source.$sb.b out=$source.$sb.qcln \
        model=$source.$sb.qcmp
$<
restor map=$source.$sb.u beam=$source.$sb.b out=$source.$sb.ucln \
        model=$source.$sb.ucmp
$<
restor map=$source.$sb.v beam=$source.$sb.b out=$source.$sb.vcln \
        model=$source.$sb.vcmp
$<
#goto terminate


DISPLAYSTOKES:
#################################################
# Display images of Stokes I, Q, U, V
#################################################
if ($source == '') then
goto POLMAP
endif
cgdisp in=$source.$sb.icln,$source.$sb.icln type=p,c device=fig1.ps/cps \
        labtyp=arcsec nxy=1,1 slev=p,10 levs1=2,3,4,5,6,7,8,9 \
        options=wedge,mirror,blacklab cols1=7 \
	region=$disreg range=0.14,14,lin,2
$<
cgdisp in=$source.$sb.qcln,$source.$sb.icln type=p,c device=fig2.ps/cps \
        labtyp=arcsec nxy=1,1 slev=p,10 levs1=2,3,4,5,6,7,8,9 \
        options=wedge,mirror,blacklab cols1=7 \
 	region=$disreg range=-0.05,-0.02,lin,-2
$<
cgdisp in=$source.$sb.ucln,$source.$sb.icln type=p,c device=fig3.ps/cps \
        labtyp=arcsec nxy=1,1 slev=p,10 levs1=2,3,4,5,6,7,8,9 \
        options=wedge,mirror,blacklab cols1=7 \
	region=$disreg range=-0.09,-0.02,lin,-2
$<
cgdisp in=$source.$sb.vcln,$source.$sb.icln type=p,c device=v.map/cps \
        labtyp=arcsec nxy=1,1 slev=p,10 levs1=2,3,4,5,6,7,8,9 \
        options=wedge,mirror region=$disreg
#goto terminate

STATISTICS:
########################################################
echo 'Make statistics from the images of stokes i, q, and u'
########################################################
$<
if ($source == '') then
goto POLMAP
endif
imstat in=$source.$sb.icln region=$statreg 
$<
imstat in=$source.$sb.qcln region=$statreg
$<
imstat in=$source.$sb.ucln region=$statreg
$<
imstat in=$source.$sb.vcln region=$statreg
#goto terminate

POLVEC:
########################################################
echo ' Make polarization intensity and vector images'
echo ' from the cleaned images of stokes i, q, and u'
########################################################
$<
if ($source == '') then
goto POLMAP
endif
$<
set irms=0.03
set qurms=0.006
\rm -fr $source.$sb.pi $source.$sb.pie
\rm -fr $source.$sb.pm $source.$sb.pme
\rm -fr $source.$sb.pa $source.$sb.pae
$<
impol in=$source.$sb.qcln,$source.$sb.ucln,$source.$sb.icln \
	poli=$source.$sb.pi,$source.$sb.pie \
	polm=$source.$sb.pm,$source.$sb.pme \
	pa=$source.$sb.pa,$source.$sb.pae \
        sigma=$qurms,$irms sncut=3,10
$<
cgdisp in=$source.$sb.pi,$source.$sb.pm,$source.$sb.pa \
	type=pix,amp,ang device=fig4.ps/cps \
        labtyp=arcsec options=mirror,blacklab \
        vecfac=2.0,3.0 cols1=0.5 \
        range=0,0,lin,1 \
        lines=5,8,8 region=$pdisreg
terminate:
echo 'end'
exit