#! /bin/csh -f
# handling standard 2 GHz polarization data: 
#                  Demonstrating to make polarization images of
#                  a target source from loading SMA data, calibration
#                  to imaging in Miriad.
# 
# Responsible: Josep Miquel Girart (CSIC-IEEC) & Jun-Hui Zhao (SAO)
# Version 2.0 
# testing data: 070728_02:39:09 from SMA archive
#
# demo illustration figures
# fig1.ps = image of polarization vector (linear lines)
#           overlaid on polarization intensity (greyscaled)
#           and continuum emission (contours) for g3141
#
##################################################################

#################################################################
#
# The correct quarter wave plate paramenters have been coded in 
# both Miriad-CVS and SMA beta Miriad  since Miriad-CVS 4.0.3
#
#################################################################
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
#
#jmg 2011-05-27 Major rearrange. The calibration steps follow a more 
#               detailed pathway. Added menus to select different options
#               for the different steps.
# 
#jhz 2011-06-02 add DATPATH
#jhz 2011-06-02 test this script and made a few minor corrections
#               delete "exit" before goto SMALOD in the MAIN MANU 
#               make hard copy fig1.ps for polarization vector
#               map of the target source.
#jhz 2011-06-03 change the script name sma2GHz1Rx_pol_src.csh
#               upgrade to 2.0
echo ''
echo 'Press return to continue'
$<

##########################################
#--> Define observational parameters:
#  date (yymmdd) 
#  time (hh:mm:ss)
#  side band (nn=0 for LSB, nn=1 for USB)**
#  reciever: rx1 for 345 Ghz 
#
# ** the calibration is done separately 
#    for the USB and LSB
##########################################
set dt    = 070728
set hour  = 02:39:09
set sb    = usb 
set nn    = 1
set rx    = rx1
set fname = "$dt""_""$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
#########################################
#--> Define reference antenna
#           average time interval
#########################################
set refant  = 2
set avetime = 5
#########################################


##########################################
#--> Define mappimg parameters:
#  polvis  
#  source
#  imsize  Image size in pixels
#  cell    Cellsize in arcsecond
#  rob     Robust weighting parameter
#  clreg   Cleaning region. it can be a ascii file with many boxes or
#          polygons defined (e.g. clreg=@filename)
#  clcut1  Cutoff for cleaning Stokes I (typically 2-3 * sigma)
#  clcut2  Cutoff for cleaning Stokes Q & U
#  disreg  Region for displaying the Stokes images  
#  statreg Region for statistical analysis (e.g. determination of sigma)
#  irms    rms noise of Stokes I
#  qurms   rms noise of Stokes Q and U
#
##########################################
set source  = $font
set polvis  = $source"_"$dt.$sb.tsys.bp.flx.gcal.pol
set imsize  = 256
set cell    = 0.15
set rob     = 0.5
set clreg   = 'arcsec,box(-3,-3,3,3)'
set clcut1  = 0.10
set clcut2  = 0.02
set disreg  = 'arcsec,box(-8,-8,8,8)'
set statreg = 'box(20,20,80,230)'
set irms    = 0.050
set qurms   = 0.008

##################################################################
# For more detailed and complete information on standard SMA 
# calibration, 
# please check: 
# http://www.cfa.harvard.edu/sma/miriad/miriad4sma/stepbystep.html 
################################################################## 
#   script based on the pol cal done of NGC 1333 IRAS 4A & G31.41+0.31:
#   Girart, Rao, Marrone 2006, Science, 313, 812                        
#   Girart, et al. 2009, Science, 324, 1408                             
echo;echo
echo "=====================================================================" 
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 (some are already included). "
echo "=====================================================================" 
echo;echo


MENU:
echo '';echo ''
echo '=========================='
echo '||        MAIN MENU     ||'
echo '=========================='
echo ''
echo '1. Load SMA data into miriad'
echo '2. Calibration'
echo '3. Mapping'
echo ''
echo ' Select one of the following options'
set sel=$<

switch ($sel)
case 1:
 echo 'Loading data into miriad' 
 goto SMALOD
breaksw
case 2:
 echo 'Standard Calibration'
 goto MENUcal
breaksw
case 3:
 echo 'Polarization map'
 goto MENUmap
breaksw
endsw
goto MENU


#####################################
# Convert raw data into miriad format
#####################################
#
SMALOD:
#
#users: please change the data path below
#
set DATPATH = /sci0/data
#
echo $fname
\rm -fr $fname
smalod in=${DATPATH}/$dt"_"$hour out=$dt rxif=1 \
   options=circular sideband=$nn refant=${refant} nscans=5,
echo 'Done. Press return to continue' $<
goto MENU



#####################################
#####################################
#
#     CALIBRATION PROCEDURES              
#
#####################################
#####################################

MENUcal:
echo;echo
echo '-------------------------------'
echo '=> SMA CAlibration procedure <='
echo '-------------------------------'
echo;echo
echo 'select (number) for the following options'
echo '   1. Inspect visibility data. Tsystem correction'
echo '   2. Flag bad data'
echo '   3. Bandpass and flux Calibration'
echo '   4. Gain Calibration for the source'
echo '   5. Polarization Calibration'
echo '   6. Go back to main menu'
echo
echo 'Select number: '; set sel=$<;echo 


MENUcal2:
echo; echo; echo 'Options:'
switch ($sel)
case 1:
 echo 'INSPECT   Create a file with a summary of the observations'
 echo '          Check for high Tsys and flag them'
 echo '          Plot visibility amplitude of all sources'
 echo 'TSYSFIX   Plot Tsys. Tsys corrections to the visibilies'
 echo; echo 'Select option (INSPECT, TSYSFIX):'
 set lab1=$<
 set out=tsys
breaksw
case 2:
 echo 'FLAG1     Flag data with incorrect pol tags'
 echo '          Flag data with Tsys and flag them'
 echo 'FLAG2     (optional in case data is corrupted)'
 echo '          Flag interatively bad data'
 echo '          If data is corrupted in a given interval, or for a given '
 echo '          antenna, baseline, etc uncomment uvflag and define properly'
 echo '          select and other parameters'
 echo 'APPLYF    Apply flags to the data'
 echo; echo 'Select option (FLAG1, FLAG2, APPLYF):'
 set lab1=$<
 set lab=tsys
 set out=tsys.flgd
breaksw
case 3:
 echo 'BPCAL     Bandpass calibration '
 echo '          Visualize solutions and data'
 echo 'APPLYBP   Apply the bandpass calibration'
 echo 'FLUX      Flux calibration'
 echo 'SPLIT     Split data'
 echo; echo 'Select option (BPCAL, APPLYBP, FLUX, SPLIT):'
 set lab1=$<
 set lab=tsys.flgd
 set out=tsys.flgd.bp
breaksw
case 4:
 echo 'GCAL      Gain calibration' 
 echo '          Visualize solutions and data'
 echo 'FLAG3     Flag data (only if necessary)'
 echo 'APPLYG    Apply gain calibration to the target'
 echo; echo 'Select option (GCAL, FLAG3, APPLYG):'
 set lab1=$<
 set lab=tsys.bp.flx
 set out=tsys.bp.flx.gcal
breaksw
case 5:
 echo 'GCAL2     First selfcalibrate pol.calibrator'
 echo '          Visualize solutions and data'
 echo 'POLAVE    Average data in 5 minutes interval for polcal'
 echo 'GETFLUX   Get the flux density for the pol calibrator'
 echo 'POLCAL    Derive the leakage solutions '
 echo 'APPLYPOL  Apply leakages to the target'
 echo; echo 'Select option (GCAL2, POLAVE, GETFLUX, POLCAL, APPLYPOL):'
 set lab1=$<
 set lab=tsys.bp.flx
 set out=tsys.bp.flx.gcal
 set out2=tsys.bp.flx.gcal.pol
breaksw
case 6:
 set lab1=MENU
breaksw
endsw
goto $lab1




#####################################
# Inspect Data
#####################################
INSPECT:
uvindex  vis=$fname >& Summary"_""$dt"".""$rx"
echo 'Generated file called: Summary_'$dt'.'$rx
echo 'Press return to continue'; $<
smacheck vis=$fname var=systemp flagval=flag range=0,1500
echo 'Press return to continue'; $<
smauvplt vis=$fname axis=time,ampl device=/xs nxy=1,1
echo 'Press return to continue'; $<
#####################################
# Tsystem correction:
#####################################
TSYSFIX:
\rm -fr    $fname.$out
smafix vis=$fname out=$fname.$out device=/xw nxy=1,1 \
       yrange=100,800 options=tsyscorr,dosour
echo 'Press return to continue'; $<
goto MENUcal




#####################################
# Flag data with incorrect pol tags
# Flag  
#####################################
FLAG1:
\rm -fr    $fname.flgd
uvflag vis=$fname.$lab select='-pol(ll,lr,rl,rr)' flagval=f
uvflag vis=$fname.$lab flagval=f edge=7 
echo 'Press return to continue'; $<
goto MENUcal2

#####################################
# Flag corrupted data: check your 
# visibility data to look for bad data
#####################################
FLAG2:
#uvflag vis=$fname.$lab flagval=f select='time(13:59,14:11),sou('$bcal')'
foreach src ($bcal $gcal)
 echo; echo 'Clik h on the screen for help'
 blflag vis=$fname.$lab axis=time,ampl device=/xs select='sou('$src')'
end
goto MENUcal2

#####################################
# Apply flags to data. create a new file
#####################################
APPLYF:
\rm -fr $fname.$out
uvcal  vis=$fname.$lab options=unflagged out=$fname.$out
echo 'Press return to continue'; $<
goto MENUcal




#####################################
# Bandpass calibration
#####################################
BPCAL:
set options=averrll
smamfcal select='pol(ll,rr),sou('$bcal')' refant=$refant \
	interval=3000 weight=2 \
        vis=$fname.$lab  options=$options edge=7 
echo 'Press return to continue. Check up the bandpass solutions'
$<
smagpplt vis=$fname.$lab select='pol(ll,rr),sou('$bcal')' \
         device=/xs yaxis=phase,amp \
         nxy=2,4 options=bandpass polyfit=5
echo 'Press return to continue'; $<
echo
echo 'Compare the raw and calibrated spectra (amplitude and phase).' 
echo 'First raw data'
echo
smauvspec vis=$fname.$lab select='pol(ll,rr),sou('$bcal')' \
          device=1/xs interval=9999 axis=freq,both nxy=4,6 options=nopass
echo 'Press return to continue. Now calibrated data'
$<
echo
smauvspec vis=$fname.$lab select='pol(ll,rr),sou('$bcal')' \
          device=2/xs interval=9999 axis=freq,both nxy=4,6
echo 'Press return to continue'; $<
goto MENUcal2

#####################################
# Apply bandpass calibration
#####################################
APPLYBP:
\rm -fr    $fname.$out
uvaver vis=$fname.$lab out=$fname.$out options=nocal
echo 'Press return to continue'; $<
goto MENUcal2

#################################################################
# Flux Calibration
#################################################################
FLUX:
echo $fcal
smaflux vis=$fname.$out select='pol(ll,rr),so('$fcal')' mirhome=$MIR
echo 'Press return to continue'; $<
goto MENUcal2

#################################################################
# Split visibility on source based.  
#################################################################
SPLIT:
foreach so ($fcal $gcal $pcal $bcal $prcal $font)
\rm -fr $so'_'$dt.$sb.tsys.bp.flx
uvaver vis=$fname.$out select='sou('$so')' \
       out=$so'_'$dt.$sb.tsys.bp.flx line=ch,24,7,114,128
end
echo 'Press return to continue'; $<
goto MENUcal




###################################################
# Gain Calibation for target. 
###################################################
GCAL:
set vis=$gcal'_'$dt.$sb.$lab 
smamfcal vis=$vis refant=${refant} \
	 options=nopassol select='pol(ll,rr)' interval=5 
echo 'Plot amp & phase solutions. Press return to continue'; $<
smagpplt vis=$vis device=1/xs yaxis=amp,pha nxy=2,4
echo 'Plot calibrated phase. Press return to continue'; $<
smauvplt vis=$vis device=1/xs axis=tim,pha nxy=2,3
echo 'Plot calibrated amp. Press return to continue'; $<
smauvplt vis=$vis device=1/xs axis=tim,pha nxy=2,3
echo 'Press return to continue'; $<
goto MENUcal2

#####################################
# Flag corrupted data: check your 
# visibility data to look for bad data
#####################################
FLAG3:
set vis=$gcal'_'$dt.$sb.$lab 
echo; echo 'Clik h on the screen for help'
blflag vis=$vis axis=time,ampl device=/xs options=noba
blflag vis=$vis axis=time,phas device=/xs options=noba
echo 'Press return to continue'; $<
goto MENUcal2

###################################################
# Applying Gain calibration to the target
###################################################
APPLYG:
gpcopy  vis=$gcal"_"$dt.$sb.$lab out=$font"_"$dt.$sb.$lab mode=copy
\rm -rf $font"_"$dt.$sb.$out
uvaver  vis=$font"_"$dt.$sb.$lab out=$font"_"$dt.$sb.$out interval=${avetime}
echo 'Press return to continue'; $<
goto MENUcal




###################################################
# Remove the phase error of pol calibrator
###################################################
GCAL2:
set vis=$pcal'_'$dt.$sb.$lab 
selfcal vis=$vis refant=${refant} \
	options=pha select='pol(ll,rr)' \
	interval=3
echo 'Plot phase solutions. Press return to continue'; $<
smagpplt vis=$vis device=/xs yaxis=phas nxy=2,3
echo 'Plot calibrated phase. Press return to continue'; $<
smauvplt vis=$vis device=/xs axis=time,phas nxy=2,3
echo 'Press return to continue'; $<
goto MENUcal2

###################################################
# Averaging visibility to a long time interval 
###################################################
POLAVE:
\rm -fr $pcal'_'$dt.$sb.$out
uvaver vis=$pcal'_'$dt.$sb.$lab out=$pcal'_'$dt.$sb.$out \
	interval=${avetime}
echo 'Press return to continue'; $<
goto MENUcal2

###################################################
# get the flux density for the pol calibrator
###################################################
GETFLUX:
uvflux vis=$pcal"_"$dt.$sb.$out stokes=i
echo 'Press return to continue'; $<
goto MENUcal2

###################################################
# Polarization calibration
###################################################
POLCAL:
set vis=$pcal'_'$dt.$sb.$out 
echo 'Delete bad points'
blflag vis=$vis axis=real,imag options=nobase \
	device=/xs
echo 'Press return to continue'; $<
gpcal vis=$vis refant=${refant} \
	options=noxy,circular,qusolve \
	interval=${avetime} flux=14.47
echo ''
echo 'Pol calibration done. First visualize uncalibrated data. '
echo 'Press return to continue'; $<
smauvplt vis=$vis device=1/xs axis=real,imag options=nocal,nopol,noba
echo 'Visualize calibrated data. Press return to continue'; $<
smauvplt vis=$vis device=2/xs axis=real,imag options=noba
echo 'Press return to continue'; $<
goto MENUcal2

###################################################
# Applying Pol calibration to the target
###################################################
APPLYPOL:
gpcopy vis=$pcal'_'$dt.$sb.$out  out=$font'_'$dt.$sb.$out options=nocal
\rm -rf $font'_'$dt.$sb.$out2
uvaver vis=$font'_'$dt.$sb.$out  out=$font'_'$dt.$sb.$out2 interval=${avetime}
echo 'Press return to continue'; $<
goto MENUcal








#####################################
#####################################
#
#  POLARIZATION IMAGING PROCEDURES        
#
#####################################
#####################################

MENUmap:
echo;echo;echo
echo '-------------------------------'
echo '=> SMA Polarization imaging  <='
echo '-------------------------------'
echo;echo
echo; echo; echo 'Options:'
echo 'DMAP      Generate dirty map for all the Stokes'
echo '          Display the dirty maps'
echo 'CLEAN     Make cleaned images all Stokes parameters'
echo 'DISPLAY   Display the cleaned maps for all Stokes parameters'
echo 'STAT      Make statistics from the Stokes images'
echo 'POLVEC    Compute polariation intensity and vector images'
echo '          Display polarization vectors and intensity.'
echo; echo 'Select option (DMAP, CLEAN, DISPLAY, STAT, POLVEC)'
set lab1=$<
goto $lab1


DMAP:
########################################
echo ' Make dirty images for four stokes'
echo ' Press return to continue'; $<
#######################################
\rm -fr $source.$sb.i $source.$sb.q $source.$sb.u $source.$sb.v $source.$sb.b
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,systemp  \
	stokes=i,q,u,v robust=$rob
echo 'Press return to continue'; $<
set dev=dirty_map.ps/cps
set dev=/xs
cgdisp in=$source.$sb.i,$source.$sb.i type=p,c \
	device=$dev \
	labtyp=arcsec nxy=1,1 slev=p,10 \
	levs1=1,2,3,4,5,6,7,8,9 \
	options=wedge,mirror region=$disreg
echo 'Press return to continue'; $<
goto MENUmap


CLEAN:
#######################################
echo ' Clean dirty images for four stokes'
echo ' Press return to continue'; $<
#######################################
\rm -fr $source.$sb.icmp
clean map=$source.$sb.i beam=$source.$sb.b out=$source.$sb.icmp \
	region=$clreg gain=0.1  niters=500 cutoff=$clcut1
foreach stk (q u v)
\rm -fr $source.$sb.$stk'cmp'
clean map=$source.$sb.$stk beam=$source.$sb.b out=$source.$sb.$stk'cmp' \
        region=$clreg gain=0.05 niters=100 cutoff=$clcut2
end

#######################################
echo ' Restore cleaned pixel images for' 
echo ' four the stokes i, q, u, and v'
echo ' Press return to continue'; $<
#######################################
foreach stk (i q u v)
\rm -fr $source.$sb.$stk'cln'
restor map=$source.$sb.$stk beam=$source.$sb.b out=$source.$sb.$stk'cln' \
	model=$source.$sb.$stk'cmp'
end
echo 'Press return to continue'; $<
goto MENUmap


DISPLAY:
#######################################
echo 'Display Stokes I image '
#######################################
set dev=$font'_'i.ps/cps 
set dev=/xs
cgdisp in=$source.$sb.icln,$source.$sb.icln type=p,c \
       device=$dev \
       labtyp=arcsec nxy=1,1 slev=p,10 levs1=1,2,3,4,5,6,7,8,9 \
       options=wedge,mirror,blacklab cols1=7 \
       region=$disreg range=0.14,14,lin,2
echo 'Press return to continue'; $<

#######################################
echo 'Display Stokes Q, U and V images'
#######################################
foreach stk (q u v)
set dev==$font'_'$stk.ps/cps 
set dev=/xs
echo 'Display Stokes '$stk
echo 'Press return to continue'; $<
cgdisp in=$source.$sb.$stk'cln',$source.$sb.icln type=c,p device=$dev \
       labtyp=arcsec nxy=1,1 slev=a,$qurms \
       levs1=-10,-8,-6,-4,-2,2,4,6,8,10 \
       options=wedge,mirror,blacklab cols1=7 \
       region=$disreg range=-0.08,14,lin,-2
end
echo 'Press return to continue'; $<
goto MENUmap


STAT:
########################################################
echo 'Make statistics from the images'
########################################################
foreach stk (i q u v)
echo 'Stokes  '$stk'. Press return to continue'; $<
imstat in=$source.$sb.$stk'cln' region=$statreg 
end
goto MENUmap


POLVEC:
########################################################
echo ' Make polarization intensity and vector images'
echo ' from the cleaned images of stokes i, q, and u'
########################################################
\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
echo 'Press return to continue'; $<
set dev=$font'_'pol.ps/cps
set dev=/xs
cgdisp in=$source.$sb.icln,$source.$sb.pi,$source.$sb.pm,$source.$sb.pa \
	type=c,pix,amp,ang device=$dev \
        labtyp=arcsec options=mirror,blacklab,rot90 \
        vecfac=2.0,3.0 cols1=0.5 \
        range=0,0,lin,1 slev=p,10 levs1=-0.5,0.5,1,2,3,4,5,6,7,8,9  \
        lines=5,8,8 region=$disreg
$<
cgdisp in=$source.$sb.icln,$source.$sb.pi,$source.$sb.pm,$source.$sb.pa \
        type=c,pix,amp,ang device=fig1.ps/cps \
        labtyp=arcsec options=mirror,blacklab,rot90 \
        vecfac=2.0,3.0 cols1=0.5 \
        range=0,0,lin,1 slev=p,10 levs1=-0.5,0.5,1,2,3,4,5,6,7,8,9  \
        lines=5,8,8 region=$disreg
echo 'Press return to continue'; $<
#–––––––––– 
exit
