#! /bin/csh -f
#
# deconv.csh - deconvolution/selfcal script. 
#
set version=24sep15
#
# Usage: deconv.csh [outfile=] [doself=] [vfirst=] [vlast=] [refant=] \
#       [start=invert/deconv/selfcal/gainapply/maplsb/cleanup] 
#
# Suggested usage:
#       a. doself=y, start=invert (selfcal the data)
#       b. doself=y, start=gainapply (check convergence of gains; OPTIONAL)
#       c. Check for position offset between maps 1,3 or 1,4 using xcor.csh
#       d. doself=n, start=maplsb (check for continuum)
#       e. doself=n, start=invert, interact=y (final maps)
#
#-----------------------------------------------------------------------------

### REQUIRED PARAMETERS
set freq=89115
set vfirst=225
set vlast=248

### Filenames
set log=deconv.log

### General parameters
set doself=y
set interact=y         # Send plots to screen as they are written?
set mosaic=y            # Use INVERT's mosaicing (don't change for SONG)
set ext=sdi             # "cln" (clean), "sdi" (mossdi), or "mem" (mosmem)

### Invert parameters
set nmaps=              # setting this overrides vlast=
set sup=0               # sidelobe suppression region
set robust=0.5          # Briggs' robustness parameter
set exinvert=           # extra parameters for INVERT
set invfwhm=            # Gaussian taper
set offset=             # reference coord for mosaic

### Deconvolution parameters and region
set clnreg=             # region command to define clean region (not mosaic)
set senmsk=2            # sensitivity level to define clean region (mosaic)
set msklev=2.5          # sensitivity level to define final image
set res=5               # approx resolution of map (arcsec), used in imfit
set usemaj=             # major axis of restoring beam (optional)
set usemin=             # minor axis of restoring beam (optional)
set usepa=              # pa of restoring beam (optional)

### Selfcal parameters
set refant=3
set interval=20
set selfrun=1
set numruns=3
set varptg=n            # Do the visfiles have different pointing grids?

### Override above parameters if given on command line
foreach a ( $* )
   set check=`echo $a | awk -F= '{print NF}'`
   if ( "$check" == 2 ) set $a
end

if ($?allvis) then
    set allvis = `echo "$allvis" | sed 's/,/ /g'`
else
    set allvis = `ls -d *.$freq.cal`
endif
echo "Files: $allvis"

### Check inputs
if ($outfile == "") then
    echo "### Fatal error: outfile keyword required."
    exit 0
endif
if ($vfirst == "") echo "### Warning: vfirst keyword not specified."
if ($vlast == "") echo "### Warning: vlast keyword not specified."

### Defaults for INVERT and MOSSDI
if ($doself == 'y') then
   if !($?cell) set cell=2
   if !($?imsize) set imsize=65
   if !($?delv) set delv=1
   if !($?rmsfac) set rmsfac=3
   set outfile=$outfile.$selfrun
else
   if !($?cell) set cell=1
   if !($?imsize) set imsize=129
   if !($?delv) set delv=0.5
   if !($?rmsfac) set rmsfac=2.5
endif

### Local settings
alias date 'date "+%m/%d/%y %T %Z"'
alias mkpdf 'gs -dBATCH -dNOPAUSE -sPAPERSIZE=letter -dPDFFitPage -q -sDEVICE=pdfwrite -sOutputFile=\!:1.pdf \!:1.ps'

set pdflist=()
if !(-e $log) then
   echo "========== deconv.csh version $version ==========" > $log
endif

if ($?start) then
   goto $start
else
   goto continue
end
continue:
# Leave continue: here if you want to run selfcal loop.

#-----------------------------------------------------------------------------
#
# 1. Generate dirty map using all visdata:

invert:
if ($nmaps == "") set nmaps=`calc -i "1+($vlast-($vfirst))/$delv"`
set chanmaps=velocity,$nmaps,$vfirst,$delv,$delv
if ($offset != "") then 
    set exinvert="offset=$offset"
endif
if ($invfwhm != "") then
    set exinvert="$exinvert fwhm=$invfwhm"
endif
rm -rf $outfile.map $outfile.beam
if ($mosaic == "y") then
    set options=systemp,double,mosaic
else
    set options=systemp,double
endif
invert vis="$allvis" map=$outfile.map beam=$outfile.beam \
        imsize=$imsize cell=$cell line=$chanmaps robust=$robust \
        sup=$sup slop=0.3 options=$options $exinvert | tee log.tmp

# Flatten noise if mosaic:
if ($mosaic == "y") then
    rm -rf $outfile.sen $outfile.gain
    mossen in=$outfile.map sen=$outfile.sen gain=$outfile.gain
    set senmin=`histo in=$outfile.sen region='image(1)'|grep Min|awk '{print $3}'`
    set clip=`calc "$senmsk*$senmin"`
    rm -rf $outfile.senmsk
    maths exp=$outfile.sen out=$outfile.senmsk mask="<$outfile.sen>.lt.$clip"
    rm -rf $outfile.sennrm $outfile.sengn
    maths exp="$outfile.sen/($senmin)" out=$outfile.sennrm
    maths exp="$outfile.sen/($outfile.gain*$senmin)" out=$outfile.sengn
    rm -rf $outfile.mapnrm
    maths exp="$outfile.map/$outfile.sennrm" out=$outfile.mapnrm
    set usefile=$outfile.mapnrm
else
    set usefile=$outfile.map
endif

# OUTPUT STATS:
echo " " | tee -a $log
echo "*** $outfile.map generated by $user on `date`" | tee -a $log
echo "  Files: $allvis" >> $log
echo "  line=$chanmaps cell=$cell imsize=$imsize exinvert=$exinvert" >> $log
grep Theoretical log.tmp >> $log
set exec="sigest in=$usefile"
set drmsi=`$exec region="image(1,5)" | grep Est | awk '{print $4}'`
set nmin4 = `expr $nmaps - 4`
set drmsf=`$exec region="image($nmin4,$nmaps)" | grep Est | awk '{print $4}'`
set drms = `calc -f f7.5 "($drmsi+$drmsf)/2"`
echo "Estimated RMS noise, $usefile : $drms" | tee -a $log
set scale=`calc -f f7.3 "$drms/$senmin"`
echo "Scaling factor for rmsfac: $scale" | tee -a $log
puthd in=$outfile.map/rms value=$drms type=real
puthd in=$outfile.map/scale value=$scale type=real
echo " " | tee -a $log

# DISPLAY:
set pltfile=$usefile
set maxval=`histo in=$pltfile | grep Max | awk '{print $3}'`
cgdisp in=$pltfile,$pltfile type=c,p device=$pltfile.ps/ps labtyp=arcsec \
    nxy=4,3 slev=a,$drms levs1=2,3,4,5,6,8,10,12,14 range=0,$maxval,lin,1 \
    csize=0.5,0.8 options=full,3value,mirr 3format=f6.1 chan=2,2 >! cgdisp.tmp
mkpdf $pltfile
set pdflist=($pdflist $pltfile.pdf)

if ($doself != "y" && $interact == "y") then
    echo -n "Continue with deconvolution [y]? "
    set resp=$<
    if ($resp == 'n') goto end
endif

#-----------------------------------------------------------------------------
#
# 2. Deconvolve maps
#

deconv:
echo "*** Deconvolved $outfile.map with $ext on `date`" | tee -a $log

# GENERATE MODEL:
if ($ext == "mem") then
    echo "  rmsfac=$rmsfac senmsk=$senmsk region=mask($outfile.senmsk)" | tee -a $log
    rm -rf $outfile.$ext
    mosmem map=$outfile.map beam=$outfile.beam out=$outfile.$ext \
        rmsfac=$rmsfac region="mask($outfile.senmsk)" flux=0.1 measure=gull 
else if ($ext == "sdi") then
    echo "  rmsfac=$rmsfac senmsk=$senmsk region=mask($outfile.senmsk)" | tee -a $log
    set rms = `gethd in=$outfile.map/rms`
    echo "The rms is $rms"
    set cutoff=`calc -f f6.4 "$rmsfac*$rms"`
    echo "Using a CLEAN cutoff of $cutoff" | tee -a $log
    rm -rf $outfile.$ext
    echo "" >> mossdi.log
    echo "MOSSDI on $outfile.map begun on `date`" | tee -a mossdi.log
    mossdi map=$outfile.map beam=$outfile.beam out=$outfile.$ext \
        cutoff=$cutoff region="mask($outfile.senmsk)" niters=50000 \
        | tee -a mossdi.log
    echo "MOSSDI finished on `date`" | tee -a mossdi.log
else if ($ext == "cln") then
    set rms = `gethd in=$outfile.map/rms`
    echo "The rms is $rms"
    set cutoff=`calc -f f7.5 "$rmsfac*$rms"`
    echo "Using a CLEAN cutoff of $cutoff" | tee -a $log
    rm -rf $outfile.$ext
    clean map=$outfile.map beam=$outfile.beam out=$outfile.$ext \
        cutoff=$cutoff region=$clnreg niters=1000
endif

# DETERMINE RESTORING BEAM:
if ($usemaj != "" && $usepa != "") then
    set bmaj=$usemaj
    set bpa=$usepa
    if ($usemin != "") then
        set bmin=$usemin
    else
        set bmin=$usemaj
    endif
else if ($mosaic == "y") then 
    rm -rf $outfile.psf
    mospsf beam=$outfile.beam out=$outfile.psf
    imfit in=$outfile.psf object=gaussian \
        "region=arcsec,box(-$res,-$res,$res,$res)" > imfit.log
    set bmaj=`grep Major imfit.log | awk '{print $4}'`
    set bmin=`grep Minor imfit.log | awk '{print $4}'`
    set bpa=`grep angle imfit.log | awk '{print $4}'`
else
    imfit in=$outfile.beam object=gaussian fix=xy \
        "region=arcsec,box(-$res,-$res,$res,$res)" > imfit.log
    set bmaj=`grep 'Major axis' imfit.log | awk '{print $4}'`
    set bmin=`grep 'Minor axis' imfit.log | awk '{print $4}'`
    set bpa=`grep angle imfit.log | grep '+/-' | awk '{print $4}'`
endif
    
# GENERATE CONVOLVED MODEL AND RESIDUALS:
rm -rf $outfile.con
restor map=$outfile.map beam=$outfile.beam model=$outfile.$ext \
        out=$outfile.con mode=convol fwhm=$bmaj,$bmin pa=$bpa
rm -rf $outfile.res
restor map=$outfile.map beam=$outfile.beam model=$outfile.$ext \
        out=$outfile.res mode=residual fwhm=$bmaj,$bmin pa=$bpa 
copyhd in=$outfile.con out=$outfile.res items=bmaj,bmin,bpa 

# DISPLAY RESIDUALS AND CONVOLVED MODEL:
  if ($mosaic == "y") then
        set region="mask($outfile.senmsk)"
     rm -rf $outfile.resnrm
     maths exp="<$outfile.res>/<$outfile.sennrm>" out=$outfile.resnrm
     echo "Plotting residuals..."
     imhist in=$outfile.resnrm region=$region options=nbin,40,log \
        device=$outfile.res_imhist.ps/ps
     set usefile=$outfile.resnrm
  else
        set region=$clnreg
     echo "Plotting residuals..."
     cgdisp in=$outfile.res device=$outfile.res_imhist.ps/ps region=$region type=p csize=0,0.8 \
            nxy=4,3 options=full,3val range=0,0,lin,3 labtyp=arcsec
     imhist in=$outfile.res region=$region options=nbin,40,log \
        device=$outfile.res_imhist.ps/ps
     set usefile=$outfile.res
  endif
mkpdf $outfile.res_imhist
set pdflist=($pdflist $outfile.res_imhist.pdf)

# CALCULATE STATS:
histo in=$usefile region="$region" | tee histo.tmp
set rrms=`grep Rms  histo.tmp | awk '{print $4}'`
set rflux=`grep Flux histo.tmp | awk '{print $6}'`
set rintflux=`calc -f f8.2 "$rflux*$delv"`
set cflux=`histo in=$outfile.con | grep Flux | awk '{print $6}'`
set cintflux=`calc -f f8.2 "$cflux*$delv"`
echo "RMS  of residuals in region: $rrms Jy/bm" | tee -a $log
echo "Flux of residuals in region: $rintflux Jy km/s" | tee -a $log
echo "Flux of deconv. model in region: $cintflux Jy km/s" | tee -a $log

if ($doself != "y" && $interact == "y") then
    echo ""
    echo -n "Produce cleaned channel maps? [y] "
    set resp = $<
    if ($resp == "n") goto end
endif

# GENERATE THE FINAL CUBE:
rm -rf $outfile.cm
restor map=$outfile.map beam=$outfile.beam model=$outfile.$ext \
        out=$outfile.cm fwhm=$bmaj,$bmin pa=$bpa 
echo "Restoring beam: $bmaj x $bmin arcsec, PA $bpa deg" >> $log
copyhd in=$outfile.map out=$outfile.cm items=mask

if ($mosaic == "y") then
    rm -rf $outfile.cmnorm
    maths exp="$outfile.cm/$outfile.sennrm" out=$outfile.cmnorm
    set usefile=$outfile.cmnorm
else
    set usefile=$outfile.cm
endif

# Calculate statistics
set nmin4 = `expr $nmaps - 4`
set exec="sigest in=$usefile"
set crmsi=`$exec region="image(1,5)" | grep Est | awk '{print $4}'`
set crmsf=`$exec region="image($nmin4,$nmaps)" | grep Est | awk '{print $4}'`
set crms = `calc -f f7.5 "($crmsi+$crmsf)/2"`
echo "Estimated RMS noise, cleaned map: $crms" | tee -a $log
puthd in=$usefile/rms value=$crms type=real

# Produce gain-corrected image cube and noise cube
rm -rf $outfile.nse
maths exp="<$outfile.sengn>*$crms" out=$outfile.nse \
    mask="<$outfile.sengn>.lt.$msklev"
rm -rf $outfile.cmmsk
maths exp="<$outfile.cm>/<$outfile.gain>" \
    region="mask($outfile.nse)" out=$outfile.cmmsk
rm -rf $outfile.cmnse
maths exp="<$outfile.nse>" region="mask($outfile.nse)" \
    out=$outfile.cmnse
foreach type (cmmsk cmnse)
    rm -f $outfile.$type.fits
    fits in=$outfile.$type out=$outfile.$type.fits op=xyout
end

# Plot the channel maps
set maxval=`histo in=$usefile | grep Max | awk '{print $3}'`
cgdisp in=$usefile,$usefile type=c,p device=$usefile.ps/cps labtyp=arcsec \
    nxy=4,3 slev=a,$crms levs1=2,3,4,5,6,8,10,12,14 range=0,$maxval,lin,1 \
    csize=0.5,0.8 options=full,3value,mirr,beambl cols1=1 3format=f6.1 \
    chan=3,3 >! cgdisp.tmp
mkpdf $usefile
imspec in=$usefile options=style,step plot=flux device=$outfile.spec.ps/ps
mkpdf $outfile.spec
set pdflist=($pdflist $usefile.pdf $outfile.spec.pdf)


if ($doself != "y") then
    if ($interact == "y") then
        echo "Plotting average spectrum..."
    endif
    goto end
endif

#-----------------------------------------------------------------------------
#
# 3. Self-calibrate:
#
# look at gpplt plots: if phases are scattered around 0, don't apply;
# if secular trends are apparent and "Average Overall Uncertainty" is
# <~ 10 deg then apply them.  Use gpedit to zero certain antennas.

selfcal:
echo "### Selfcal run: $selfrun/$numruns begins on `date`" | tee -a $log
echo " interval=$interval  refant=$refant" | tee -a $log

if ($mosaic == "y" && $varptg == "n") then
    rm -rf $outfile:r.mod*
    demos map=$outfile.$ext vis=$allvis[1] out=$outfile:r.mod options=detaper
endif

# LOOP OVER VISIBILITY FILES:
set filenum=1
foreach vis ($allvis)
    if ($mosaic == "y") then
        if ($varptg == "y") then
           rm -rf $outfile:r.mod*
           demos map=$outfile.$ext vis=$vis out=$outfile:r.mod options=detaper
        endif
        selfcal vis=$vis interval=$interval model="$outfile:r.mod*" \
            refant=$refant options=mosaic,phase | tee scal.tmp
    else
        selfcal vis=$vis interval=$interval model="$outfile.$ext" \
            refant=$refant options=phase | tee scal.tmp
    endif
    echo "Processing file $vis..."
    tail -15 scal.tmp >> $log
    gpplt vis=$vis device=/ps yaxis=phase yrange=-180,180 nxy=3,2 \
        options=w
    mv pgplot.ps gains$selfrun.$vis.ps
    if ($interact == "y") gv gains$selfrun.$vis.ps &
    set filenum=`expr $filenum + 1`
    echo " " | tee -a $log
end

if ($doself == "y") then
    set selfrun=`expr $selfrun + 1`
    set outfile=$outfile:r.$selfrun
    if ($selfrun <= $numruns) goto invert
endif

goto end

### ----------------------- END OF SELFCAL LOOP

#-----------------------------------------------------------------------------
#
# 4. Final check: apply the gains and selfcal again.  Resulting gains
#    should be close to zero.

gainapply:
foreach vis ($allvis)
    set prefix=$vis:r
    rm -rf $prefix.cal2
    uvcat vis=$vis out=$prefix.cal2
end
set numruns=`expr $numruns + 1`
set selfrun=$numruns
set outfile=$outfile:r.$selfrun
set allvis = `ls -d *.$freq.cal2`
goto invert


#-----------------------------------------------------------------------------
#
# 6. Cleanup:

cleanup:
echo "outfile: $outfile"
rm -rf $outfile.[1-4].*
rm -rf $outfile.mod*
rm -rf $outfile.res* $outfile.sen*
rm -rf $outfile.con $outfile.nse $outfile.gain
rm -rf $outfile.psf $outfile.mapnrm
rm -f imfit.log
foreach vis ($allvis)
    rm -rf $vis:r.cal2
end
gzip -f *.ps

end:
rm -f *.tmp
if ($#pdflist != 0) then
    gs -q -dSAFER -dNOPAUSE -dBATCH -sOutputFile=$outfile.pdf \
        -sDEVICE=pdfwrite -c .setpdfwrite -f $pdflist
    echo Plotted images output to $outfile.pdf
    rm -f $pdflist
endif