#! /bin/csh -f
# handling standard 4 GHz polarization data 
# 
# Responsible: testing version to be used for 4GHz polarization system
# Version 1.0 
# testing data: 100505_04:45:58 from SMA Ken (Taco) Young, PI: Keping Qiu
#               The testing data set has not been available to the public yet.
#               Please acquire the PI's permission for use of the testing data
#               and the figures produced from this script.
#
# demo illustration figures (under construction)
# 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)
#
#################################################################################### 
#
#jhz 2011-02-17 received the test data from Taco
#jhz 2011-02-17 received the c-code of extract baseline-based 
#               information regarding to the polarization visibility
#               
#jhz 2011-02-17 communicated with users for information on the 4GHz pol data and test
#
#jhz 2011-02-21 look at the data with mirREAD routines outside Miriad
#jhz 2011-03-02 come back to look into the polarization section with mirRead 
#               routines for 4 GHz polarization data; pol states and pol 
#               sampling patterns are difficult to understand; the documents
#               available appeared to inadequate for disentangle them.
#
#jhz 2011-03-03 start the script (Ver0.1) for testing the 4GHz polarization 
#              system 
#
#jhz 2011-03-03 work 4GHz polarization codes in Miriad
#jhz 2011-03-23 go through compiling
#jhz 2011-03-25 smalod reads all 4GHz polarization data
#jhz 2011-03-29 check the variables for polarization calibration
#               realized a few peculiar problems.
#jhz 2011-03-30 parallactic angles for the 4GHz data are inconsistent between 
#               the values from smalod and uvredo ....
#jhz 2011-03-31 polarization labelling for the 4GHz data does not make sense.              
#               check evector OK
#               check mount type OK
#jhz 2011-04-04 install SMA Miriad beta for 4 GHz pol test code on Jupiter
#jhz 2011-04-05 make a tarball for the 4GHz pol test Miriad to SMA Miriad
#               web site for testers
#
#test with 2GHz 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-28 received 7 variable tables of the 4GHz pol test data set
#               100505_04:45:58 from online software crew:
#               ....
#               Attached are seven tables.   The first gives the scan number, 
#               source name, parallactic angle (degrees), integration time 
#               (seconds) and UT (hours) for each scan. Each of the other 
#               6 tables 
#               gives the polarization state for one baseline to antenna 1 
#               (table2 has 1-2 data, table 3 has 1-3 data etc., there is 
#               no table for 1-7 as 7 was not in the array). For the
#               last 6 tables the format is scan number, source name, 
#               baseline name ("1-2" etc.) and polarization code. The 
#               polarization state for any baseline can be derived from 
#               the information in these tables. All of these tables
#               contain entries for every scan, including the scans 
#               which are flagged bad on one or more baselines.
#               Taco
#               
#jhz 2011-05-10 verified with Ken Young, bl.ipol, the baseline pol ids in 
#               the Mir formatted data for 4 GHz polarization mode with 
#               QWP on single receivers are related to circular pol states:
#               bl.ipol  pol states
#               0        Unknown
#               1        RR
#               2        RL
#               3        LR
#               4        LL
#
#jhz 2011-05-11 merge 2GHz script Ver1.2 with 4GHz Ver0.1 to Ver0.2
#               correct pol labelling and several other minor bugs.
#               recompile the package and upgrate to beta 1.3.0
#               install Ver 1.3.0 on Jupiter (centOS) and bigfoot (Redhat)
#               test smalod on both Jupiter bigfoot successfully.
#               update SMA Miriad web for testers downloading
#
#jhz 2011-05-12 add back chi2 in variables 
#               add lineset, script variable for linetype
#               change stokes = i to =ii in uvflux
#               add smasplt
#               add a few script variables, nw, spltopt
#               add goto tag SMASPLT for 4GH data
#               add gcal2 and fcal2
#
#jhz 2011-05-12 add goto tag BLCORR to correct for baseline errors
#               work out setup to correct for baseline errors.    
#
#jhz 2011-05-13 add goto tag CHI
#jhz 2011-05-18 update ephemeris for the solar objects used in smaflux 
#jhz 2011-05-19 update sma/mount = NASMYTH
#jhz 2011-06-03 implement the imaging part
#jhz 2011-06-09 changed the data path to CF where the test data located.
#jhz 2011-06-13 following the suggestions from the software meeting
#               today, comment out the 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 100505 04:45:58"
goto terminate
endif
##########################################
#--> Define date, time ofthe observations
##########################################
set date=$1
set hour=$2
goto USERSREDUCTION
SMATEST:
##########################################
# for SMA software test
##########################################
set date=100505
set hour=04:45:58
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  = titan
set fcal2 = neptune
set gcal  = 1830+063
set gcal2 = 1751+096 
set bcal  = 3c273  
set pcal  = 3c273
set prcal = 3c454.3
#########################################
#-->  Define the Source name
#########################################
set font  = I18360 
# select source for pol mapping
set source = ''       
#########################################
#--> Define reference antenna
#           average time interval
#########################################
set refant  = 5
set avetime = 6
#########################################
#--> Other definitions
set region  = quart
set regdisp = 'arcsec,box(-12,-12,12,12)'
set dt      = flx
set spltops = all 
set nw      = 48
set lineset = 'chan,48,6,116,128'
set edg     = '6,6'
set bc      = blc
#
#########################################
#
#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 IPOLFLAG      #Flag data with incorrect pol tags
#goto TSYSFIX       #Corrections for fringe amplitude
#goto SMASPLT       #Separate the ch0 data from spectral chunk data 
#goto FLAG          #Flag bad visibilities
#goto BPCAL         #Bandpass calibration
#goto APPLYBP       #Apply bandpass solutions
#goto CHI           #re-calculate chi with uvredo
#goto FLUX          #Determine flux density scale
#goto SPLIT         #Split the data into single-source file
#goto BLCORR        #Corrections for baseline errors
#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
#
#users: please change the data path below
set DATPATH = /sma/SMAusers/jzhao/data/
#
#####################################
# copy the antenna file to working area
#####################################
#
cp ${DATPATH}/$date"_"$hour/antennas .
#
#####################################
# Convert raw data into miriad format
#####################################
#
SMALOD:
echo $fname
\rm -fr $fname
smalod in=${DATPATH}/$date"_"$hour out=$date rxif=1 \
       options=circular sideband=$nn refant=${refant} nscans=5,
#goto terminate
$<
echo 'Press return to continue'
$<
#
#####################################
# Inspect at 
#####################################
INSPECT:
uvindex  vis=$fname
$<
smacheck vis=$fname var=systemp flagval=flag range=0,1000
$<
smauvplt vis=$fname axis=time,ampl device=/xs nxy=1,1
$<
#goto terminate
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

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

#####################################
#Separate the ch0 from spectral windows
#####################################
SMASPLT:
\rm -r $fname.tsys.$nw
smasplt vis=$fname.tsys  out=$fname.tsys.$nw options=$spltops
#####################################
# Flag corrupted data: check your 
# visibility data to look for bad data
#####################################
FLAG:
#####################################
#flag edge channels 
#####################################
uvflag vis=$fname.tsys.$nw flagval=f edge=$edg 
$<
#####################################
#flag bad time interval for selected source
#####################################
uvflag vis=$fname.tsys.$nw flagval=f select='time(13:59,14:11),sou('$prcal')'
uvflag vis=$fname.tsys.$nw flagval=f select='time(06:00,06:25),sou('$bcal')'
$<
#####################################
#flag bad  antennas
#####################################
uvflag vis=$fname select='ant(8),time(12:30:00,13:20:00)' flagval=f
$<
#
# users can enable blflag for flagging the bad visibilities.
#
#blflag vis=$fname.tsys.$nw axis=time,ampl device=/xs select='sou('$bcal')'
#
#goto terminate

#####################################
# Bandpass calibration
#####################################
BPCAL:
set options=averrll
smamfcal select='pol(ll,rr),sou('$bcal')' refant=$refant \
	interval=3000 weight=2 \
        vis=$fname.tsys.$nw options=$options edge=$edg
$<
smagpplt vis=$fname.tsys.$nw  device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=1,8

#goto terminate

#####################################
# Apply bandpass calibration
#####################################
APPLYBP:
\rm -fr    $fname.tsys.$nw.bp
uvaver vis=$fname.tsys.$nw out=$fname.tsys.$nw.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.$nw.bp/evector value=0.7853
# which has been coded in Miriad
#puthd in=$fname.tsys.$nw.bp/mount value=4
# which has been coded in Miriad
CHI:
################################
# re-calculate chi with uv-redo
################################
set chi = ''
#set chi = '.chi'
#\rm -fr    $fname.tsys.$nw.bp$chi
#uvredo vis=$fname.tsys.$nw.bp options=chi out=$fname.tsys.$nw.bp$chi
# which has been coded in Miriad
#
#goto terminate
#################################################################
#
#################################################################
# Flux Calibration
#################################################################
#
FLUX:
echo 'FLux calibrator =' $fcal
smaflux vis=$fname.tsys.$nw.bp$chi select='pol(ll,rr),so('$fcal')' mirhome=$MIR
#goto terminate
#
$<
echo 'Press return to continue'
$<
#
#################################################################
# Split visibility on source based.  
#################################################################
#goto terminate
#
SPLIT:
set so=$fcal
\rm -fr "$so""_""$dt".$sb
uvaver vis=$fname.tsys.$nw.bp$chi select='sou('$so')' \
       out="$so""_""$dt".$sb line=$lineset 
set so=$fcal2
\rm -fr "$so""_""$dt".$sb
uvaver vis=$fname.tsys.$nw.bp$chi select='sou('$so')' \
       out="$so""_""$dt".$sb line=$lineset
set so=$gcal
\rm -fr "$so""_""$dt".$sb
uvaver vis=$fname.tsys.$nw.bp$chi select='sou('$so')' \
       out="$so""_""$dt".$sb line=$lineset
set so=$gcal2
\rm -fr "$so""_""$dt".$sb
uvaver vis=$fname.tsys.$nw.bp$chi select='sou('$so')' \
       out="$so""_""$dt".$sb line=$lineset
set so=$bcal
\rm -fr "$so""_""$dt".$sb
uvaver vis=$fname.tsys.$nw.bp$chi select='sou('$so')' \
       out="$so""_""$dt".$sb line=$lineset
set so=$pcal
\rm -fr p"$so""_""$dt".$sb
uvaver vis=$fname.tsys.$nw.bp$chi select='sou('$so')' \
       out=p"$so""_""$dt".$sb line=$lineset
set so=$prcal
\rm -fr "$so""_""$dt".$sb
uvaver vis=$fname.tsys.$nw.bp$chi select='sou('$so')' \
       out="$so""_""$dt".$sb line=$lineset
set so=$font
\rm -fr "$so""_""$dt".$sb
uvaver vis=$fname.tsys.$nw.bp$chi select='sou('$so')' \
       out="$so""_""$dt".$sb line=$lineset
#goto terminate

BLCORR:
\rm -rf p$pcal"_"$dt.$sb.$bc
uvedit vis=p$pcal"_"$dt.$sb source=$pcal apfile=antennas \
	smaoffset=1,-0.085,0.036,0.046,2,0.110,-0.092,0.049,3,-0.067,0.049,0.046,4,-0.088,-0.025,0.032,5,0.088,-0.032,0.064 \
	out=p$pcal"_"$dt.$sb.$bc options=sma
#goto terminate

GETFLUX:
###################################################
# Remove the phase error
###################################################
selfcal vis=p$pcal"_"$dt.$sb.$bc refant=${refant} \
	options=pha select='pol(ll,rr)' \
	interval=${avetime} 

###################################################
# get the flux density for the pol calibrator
###################################################
uvflux vis=p$pcal"_"$dt.$sb.$bc stokes=ii

###################################################
# Averaging visibility to a long time interval 
###################################################
\rm -fr pp$pcal"_"$dt.$sb.$bc
uvaver vis=p$pcal"_"$dt.$sb.$bc out=pp$pcal"_"$dt.$sb.$bc \
	interval=${avetime}
#goto terminate

###################################################
# Polarization calibration
###################################################
POLCAL:
$<
gpcal vis=pp$pcal"_"$dt.$sb.$bc refant=${refant} \
	options=noxy,circular,qusolve \
	interval=${avetime} flux=5.5
#goto terminate

###################################################
# Gain Calibation for target. 
# "smamfcal options=nopassol" can also be used 
# instead of selfcal
###################################################
GCAL:
$<
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

###################################################
# 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.$bc out=$font"_"$dt.$sb.g options=nocal
#goto terminate

AVER:
\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=0.5
$<
\rm -fr $prcal"_"$dt.$sb.g
uvaver vis=$prcal"_"$dt.$sb out=$prcal"_"$dt.$sb.g \
	interval=${avetime}
gpcopy  vis=pp$pcal"_"$dt.$sb.$bc 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
##################################################
echo ' Polarization Imaging'
##################################################
POLMAP:
#########################
# select source for mapping
#########################
set source=$pcal

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

#########################
# 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)'
set statcet='arcsec,box(-2,-2,2,2)'
#goto terminate


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


DIRTYMAP:
if ($source == '') then
goto POLMAP
endif

\rm -fr $source.$sb.i $source.$sb.q $source.$sb.u $source.$sb.v $source.$sb.b
########################################
#make sure to have a right vis file
#corresponding to the polarization source 
#for imaging ->  invert vis=????????????
########################################
#
#######################################
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,$nw,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

#######################################
echo ' Clean dirty images for four stokes'
#######################################
$<
CLEANMAP:
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=100
$<
mfclean map=$source.$sb.u beam=$source.$sb.b out=$source.$sb.ucmp \
        region=$clreg gain=0.05 niters=100
$<
mfclean map=$source.$sb.v beam=$source.$sb.b out=$source.$sb.vcmp \
        region=$clreg gain=0.05 niters=100
#goto terminate

#######################################
echo ' Restore cleaned pixel images for'
echo ' four the stokes i, q, u, and v'
#######################################
$<
RESTOR:
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:
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.1,4,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.003,0.024,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.003,-0.024,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



########################################################
echo 'Make statistics from the images of stokes i, q, and u'
########################################################
$<
STATISTICS:
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
$<
imstat in=$source.$sb.qcln region=$statcet
$<
imstat in=$source.$sb.ucln region=$statcet

#goto terminate

########################################################
echo ' Make polarization intensity and vector images'
echo ' from the cleaned images of stokes i, q, and u'
########################################################
$<
POLVEC:
if ($source == '') then
goto POLMAP
endif
$<
set irms=0.01
set qurms=0.005
\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
