#!/bin/csh -f
#
# edgeimage.csh - CARMA EDGE imaging script (Tony Wong, UIUC)
#
# All parameters can be set on the command line:
#   edgeimage.csh files= source= vsys= delv= vwidth=
#
# History
#  2014nov20: original
#  2014dec05: HTML log output
#  2014dec07: bug fixes, choose min noise from both ends of spectrum
#  2014dec08: reduce maxsen to 2.5 and msklev to 2
#  2014dec18: do not output sengain as fits cube
#  2015jan21: output mJy as well as K rms to HTML file
#  2015feb25: bug fixes
#  2015mar02: allow robust to be unset
#  2015apr02: allow continuum subtraction
#  2015apr05: adjust algorithm for calculating restricted velocity range
#  2015apr28: further flexibility in choosing vrange; dvvis=5
#  2015jul14: allow more flexible correlator settings (pre-MSIP)
#  2015sep03: output cmnorm to fits, and trim blank space at edges
#  2016jun08: output cmnormgain, copy beam info to noise images; tpeak in K
#--------------------------------------------------------------------------

# Basic params
set files=
set source=
set showplot=n
set purge=y
set arr=e           # Array label for HTML link
set dropch=1        # Edge channels to avoid
set contsub=n       # Perform continuum subtraction on vis spectra?
set dvvis=5         # channel spacing for visfile in km/s
set savevis=y       # Export visibilities in a tarfile?

# Mapping and deconvolution params
set vsys=           # systemic velocity [REQ'D]
set delv=           # channel spacing for image cube in km/s [REQ'D]
set vwidth=         # total velocity range to image [REQ'D]
set vfirst=         # velocity of first channel [optional, auto selected]
set nmaps=          # number of channels [optional, auto selected]
set imsize=257      # for INVERT
set cell=1          # for INVERT
set robust=0.5      # for INVERT
set sup=            # for INVERT
set offset=         # for INVERT
set exinvert=       # extra parameters for INVERT
set rmsfac=2.5      # rmsfac for mossdi
set clnbox=         # deconvolution region, senmsk if unset
set maxsen=2.5      # sensitivity contour level for deconvolution
set msklev=2        # sensitivity contour level for signal detection
set bmaj=           # beam major axis to override default, arcsec
set bmin=           # beam minor axis to override default, arcsec
set bpa=            # beam position angle to override default, deg
set decon='mossdi2' # deconvolution task, alternative: mosmem
set exmem=          # extra keyword for mosmem
set sdmap=          # sd map for joint deconvolution
set sdbeam=         # sd beam for joint deconvolution

# Moment params
set fwhm=15         # FWHM of smoothed map for masking method
set mom0lev=1       # First contour level for mom0, in K km/s
set mom1del=20      # Spacing between velocity contours, in km/s
set clipsmo=2.5     # Sigma clip level to make mask for smoothed cube
set clip1pk=3.5     # Sigma clip level for vpeak image

### LOCAL SETTINGS
alias date 'date "+%m/%d/%y %T %Z"'
set version=2016jun08
set log=edgeimage.log

### OVERRIDE ABOVE PARAMETERS IF GIVEN ON COMMAND LINE
foreach par ( $* )
   set check=`echo $par | awk -F= '{print NF}'`
   if ( "$check" >= 2 ) set $par
end

### CHECK THAT CUBE IS FULLY SPECIFIED
if ($vsys == "" || $delv == "" || $vwidth == "") then
    echo "### FATAL ERROR: vsys, delv, and vwidth must be specified"
    exit 0
endif

#============================================================================

# Output some interesting information
echo ""
echo "edgeimage.csh $version executed by $user on `date`" | tee -a $log
echo "" | tee -a $log

#--------------------------------------------------------------------------
# 0. Generate concatenated visibility file
#--------------------------------------------------------------------------

# Determine the available velocity range
set filelis=`echo "$files" | sed 's/,/ /g'`
set restfreq=`prthd in=$filelis[1] | grep GHz | head -1 | awk '{print $5}'`
if ($restfreq =~ 110*) then
    set line=13co
    echo "This is a 13CO dataset" | tee -a $log
else
    set line=co
    echo "This is a CO dataset" | tee -a $log
endif

# Determine the available velocity range
set allveli=()
set allvelf=()
foreach visfile ($filelis)
    uvspec vis=$visfile axis=vel,amp log=tmp.uvspec1 interval=10000 options=nobase,avall
    sort -n tmp.uvspec1 > tmp.uvspec2
    set drphead=`calc -i "1+$dropch"`
    set veli=`head -$drphead tmp.uvspec2 | tail -1 | awk '{printf "%.1f", $1}'`
    set velf=`tail -$drphead tmp.uvspec2 | head -1 | awk '{printf "%.1f", $1}'`
    echo "$visfile has velocity range $veli to $velf" | tee -a $log
    set allveli=($allveli $veli)
    set allvelf=($allvelf $velf)
end
set veli=`echo $allveli | tr ' ' '\n' | sort -k1,1n | tail -1`
set velf=`echo $allvelf | tr ' ' '\n' | sort -k1,1n | head -1`
echo "Common velocity range is $veli to $velf" | tee -a $log
set vvisi=`echo $veli | awk -v dv=$dvvis '{printf "%.0f\n", ($1+dv)/dv}' | awk -v dv=$dvvis '{print $1*dv}'`
set vvisf=`echo $velf | awk -v dv=$dvvis '{printf "%.0f\n", ($1-dv)/dv}' | awk -v dv=$dvvis '{print $1*dv}'`
echo "Visibility channels will be $vvisi to $vvisf spaced by $dvvis km/s" | tee -a $log
set nchvis=`calc -i "1+($vvisf-$vvisi)/$dvvis"`
set tpc=`calc -i "($nchvis*0.1)+0.5"`
set npc=`calc -i "($nchvis*0.9)+0.5"`

# Generate averaged file and subtract continuum if needed
rm -rf $source.vis
uvaver vis="$files" out=$source.vis line=vel,$nchvis,$vvisi,$dvvis,$dvvis
if ($contsub == 'y') then
    echo "Fitting continuum to channels 1,$tpc,$npc,$nchvis" | tee -a $log
    rm -rf $source.uvlin $source.chan0
    uvlin vis=$source.vis out=$source.uvlin chans=1,$tpc,$npc,$nchvis
    uvlin vis=$source.vis out=$source.chan0 chans=1,$tpc,$npc,$nchvis \
        mode=chan0
endif

# Make a list of original track files
set rawlis = ()
set i = 1
foreach file ($filelis)
    set raw = `grep miriad $file/history | head -1 | awk -F= '{print $2}'`
    set rawlis = ($rawlis $raw)
    set i = `expr $i + 1`
end


#--------------------------------------------------------------------------
# 1. Make dirty map:
#--------------------------------------------------------------------------

# Set up the default velocity range
if ($nmaps == '') set nmaps=`calc -i "1+$vwidth/$delv"`
set vtmp = `calc -i "$vsys/$delv-($nmaps-1)/2"`
if ($vfirst == '') set vfirst = `expr $vtmp \* $delv`
set vlast=`calc -i "$vfirst+$delv*($nmaps-1)"`

# Check for any hard limits based on spectral coverage
if ($vfirst < $vvisi) then
    set nclip=`echo $vfirst $vvisi $delv | awk '{printf "%.0f\n",0.51+($2-$1)/$3}'`
    set vfirst=`calc -i "$vfirst+$nclip*$delv"`
    set nmaps=`calc -i "$nmaps-$nclip"`
    echo "Dropping $nclip channels to bring vfirst to $vfirst" | tee -a $log
endif
if ($vlast > $vvisf) then
    set nclip=`echo $vlast $vvisf $delv | awk '{printf "%.0f\n",0.51+($1-$2)/$3}'`
    set vlast=`calc -i "$vlast-$nclip*$delv"`
    set nmaps=`calc -i "$nmaps-$nclip"`
    echo "Dropping $nclip channels to bring vlast to $vlast" | tee -a $log
endif

# Make the map
rm -rf $source.map $source.beam
echo "Included files: $files" | tee -a $log
echo "The source velocity appears to be $vsys" | tee -a $log
echo "Making $nmaps maps from $vfirst to $vlast and spaced by $delv km/s" \
    | tee -a $log
echo "" | tee -a $log

if ($sup != '') then
    set exinvert="sup=$sup $exinvert"
endif
if ($robust != '') then
    set exinvert="robust=$robust $exinvert"
endif
if ($offset != '') then
    set exinvert="offset=$offset $exinvert"
endif
if ($exinvert != '') then
    echo "Extra invert parameters: $exinvert" | tee -a $log
endif

if ($contsub == 'y') then
    set infile=$source.uvlin
else
    set infile=$source.vis
endif
invert vis=$infile map=$source.map beam=$source.beam cell=$cell \
    imsize=$imsize line=velocity,$nmaps,$vfirst,$delv,$delv \
    options=systemp,double,mosaic slop=1 $exinvert \
    | tee tmp.invert

# Get the number of pointings
grep Number tmp.invert >> $log
set nptng=`grep Number tmp.invert | awk '{print $4/3}'`
puthd in=$source.map/nptng value=$nptng type=integer

# Get the (flattened) noise level.
rm -rf $source.sen $source.gain
mossen in=$source.map sen=$source.sen gain=$source.gain
set minval=`histo in=$source.sen region='image(1)'|grep Min|awk '{print $3}'`
echo "Normalizing $source.sen by $minval" | tee -a $log
rm -rf $source.sennorm $source.sengain
maths exp="$source.sen/($minval)" out=$source.sennorm
maths exp="$source.sen/($source.gain*$minval)" out=$source.sengain

# Normalize the maps
rm -rf $source.mapnorm
maths exp="$source.map/$source.sennorm" out=$source.mapnorm
rm -rf $source.senmsk
maths exp=$source.sen out=$source.senmsk mask="<$source.sengain>.lt.$maxsen"
set nmin1 = `expr $nmaps - 1`
set exec="sigest in=$source.mapnorm"
set drmsi=`$exec region="image(1,2)" | grep Est | awk '{print $4}'`
set drmsf=`$exec region="image($nmin1,$nmaps)" | grep Est | awk '{print $4}'`
#set drms = `calc -f f7.5 "($drmsi+$drmsf)/2"`
set drms = `echo $drmsi $drmsf | awk '{if ($1<$2) {print $1} else {print $2}}'`
echo "Estimated RMS noise, dirty map: $drms" | tee -a $log
set scale=`calc -f f7.3 "$drms/$minval"`
echo "Scaling factor for rmsfac: $scale" | tee -a $log
puthd in=$source.mapnorm/rms value=$drms type=real
puthd in=$source.mapnorm/scale value=$scale type=real
echo " " | tee -a $log

# Make the plot.
cgdisp in=$source.mapnorm device=$source.dmaps.ps/ps \
    region='arc,box(-90,-90,90,90)' \
    options=full,3value labtyp=arcsec nxy=5,4 type=contour \
    slev=a,$drms levs1=2,3,4,5,6,8,10,12,14 csize=0.5,0.7
if ($showplot != 'n') then
    gv $source.dmaps.ps &
endif 

#--------------------------------------------------------------------------
# 2. Clean the map:
#--------------------------------------------------------------------------

deconv:
# Deconvolve using MOSSDI
set scale=`gethd in=$source.mapnorm/scale`
set cutoff=`calc "$scale*$rmsfac"`
echo "${decon}: rmsfac=$rmsfac, cutoff=$cutoff" | tee -a $log

rm -rf $source.psf
mospsf beam=$source.beam out=$source.psf
if ($bmaj == '') then
    imfit in=$source.psf object=beam region="a,box(-10,-10,10,10)" \
        | tee tmp.imfit
    set bmaj=`grep 'Major axis' tmp.imfit | sed 's/.*://' | awk '{printf "%.2f",$1}'`
    set bmin=`grep 'Minor axis' tmp.imfit | sed 's/.*://' | awk '{printf "%.2f",$1}'`
    set bpa=`grep 'Position angle' tmp.imfit | sed 's/.*://' | awk '{printf "%.1f",$1}'`
else
    if ($bmin == '') set bmin=$bmaj
    if ($bpa == '') set bpa=0
endif

rm -rf $source.cln

if ($decon == 'mossdi2') then
    if ($clnbox == '') then
        mossdi2 map=$source.map beam=$source.beam out=$source.cln \
            cutoff=$cutoff niters=20000 region="mask($source.senmsk)"
    else
        mossdi2 map=$source.map beam=$source.beam out=$source.cln \
            cutoff=$cutoff niters=20000 region=$clnbox
    endif
endif

if ($decon == 'mosmem') then
    rm -rf $sdmap.rg
    regrid in=$sdmap out=$sdmap.rg tin=$source.map
    if ($clnbox == '') then
        mosmem map=$source.map,$sdmap.rg beam=$source.beam,$sdbeam \
            out=$source.cln rmsfac=$scale \
            niters=20000 region="mask($source.senmsk)" $exmem
    else
        mosmem map=$source.map,$sdmap.rg beam=$source.beam,$sdbeam \
            out=$source.cln rmsfac=$scale \
            niters=20000 region=$clnbox $exmem
    endif
endif

# Produced cleaned maps and noise-flattened image cube (cm, cmnorm, cmnormgain)
rm -rf $source.cm
restor map=$source.map beam=$source.beam model=$source.cln out=$source.cm \
    fwhm=$bmaj,$bmin pa=$bpa
copyhd in=$source.map out=$source.cm items=mask
rm -rf $source.cmnorm $source.cmnorm.flux
maths exp="$source.cm/$source.sennorm" out=$source.cmnorm region="mask($source.cm)"
maths exp="($source.gain*$minval)/$source.sen" out=$source.cmnormgain \
    region="mask($source.cm)"
copyhd in=$source.cmnorm out=$source.cmnormgain items=bmaj,bmin,bpa
set nmin4 = `expr $nmaps - 4`
set exec="sigest in=$source.cmnorm"
set crmsi=`$exec region="image(1,4)" | 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"`
#set crms = `echo $crmsi $crmsf | awk '{if ($1<$2) {print $1} else {print $2}}'`
echo "Estimated RMS noise, cleaned map: $crms" | tee -a $log
puthd in=$source.cmnorm/rms value=$crms type=real

# Produce gain-corrected image cube (cmmsk) and noise cube (cmnse)
set minval=`histo in=$source.sen region='image(1)'|grep Min|awk '{print $3}'`
rm -rf $source.nse $source.nsesub $source.nse1
maths exp="<$source.sengain>*$crms" out=$source.nse \
    mask="<$source.sengain>.lt.$msklev"
imsub in=$source.nse region='image(1)' out=$source.nsesub
imblr in=$source.nsesub out=$source.nse1 value=1
rm -rf $source.cmmsk
maths exp="<$source.cm>/<$source.gain>" region="mask($source.nse)" \
    out=$source.cmmsk
rm -rf $source.cmnse
maths exp="<$source.nse>" region="mask($source.nse)" out=$source.cmnse
copyhd in=$source.cmmsk out=$source.cmnse items=bmaj,bmin,bpa
set ra=`gethd in=$source.cm/crval1 format=hms`
set dc=`gethd in=$source.cm/crval2 format=dms`
set dvhel=`cotra radec=$ra,$dc uvw=10.3,15.3,7.7 | grep dVlsr | awk '{print (-1)*$2}'`
set wdate=`gdate -u '+%Y-%m-%d'`
foreach type (cmnorm cmnormgain cmmsk cmnse)
    puthd in=$source.$type/dvhel value=$dvhel type=real
    puthd in=$source.$type/date value=$wdate type=ascii
    rm -f $source.$type.fits
    fits in=$source.$type out=$source.$type.fits op=xyout
end

# Produced noise-flattened residual cube (resnorm)
rm -rf $source.res
restor map=$source.map beam=$source.beam model=$source.cln out=$source.res \
    fwhm=$bmaj,$bmin pa=$bpa mode=residual
rm -rf $source.resnorm
maths exp="$source.res/$source.sennorm" out=$source.resnorm
copyhd in=$source.cm out=$source.resnorm items=bmaj,bmin,bpa
set rflux=`histo in=$source.resnorm | grep Flux | awk '{print $6}'`
set rrms=`sigest in=$source.resnorm | grep Est | awk -v c=$crms '{print $4/c}'`
set rintflux=`calc -f f8.2 "$rflux*$delv"`
set cflux=`histo in=$source.cln | grep Flux | awk '{print $6}'`
set cintflux=`calc -f f8.2 "$cflux*$delv"`
echo "Used restoring beam of $bmaj x $bmin arcsec, pa $bpa" | tee -a $log
echo "Flux in residuals map: $rintflux Jy km/s" | tee -a $log
echo "Flux in deconv. model: $cintflux Jy km/s" | tee -a $log
echo "RMS of residual in sigma units: $rrms" | tee -a $log
echo " " | tee -a $log

# Make the plot.
set maxval=`histo in=$source.cmnorm | grep Max | awk '{print $3}'`
cgdisp in=$source.cmnorm,$source.cmnorm device=$source.cmaps.ps/ps \
    region='arc,box(-90,-90,90,90)' options=full,3value,mirr \
    type=c,p labtyp=arcsec nxy=5,4 lines=1,1 slev=a,$crms \
    levs1=2,3,4,5,6,8,10,12,14 csize=0.5,0.7 range=0,$maxval,lin,1
imspec in=$source.cmmsk plot=flux device=$source.imspec.ps/ps \
    region='arc,box(-20,-20,20,20)' options=style,step
if ($showplot != 'n') then
    gv $source.cmaps.ps &
    gv $source.imspec.ps &
endif

#--------------------------------------------------------------------------
# 3. Make moment maps:
#--------------------------------------------------------------------------

domom:

# Calculate Jy to K conversion
if !($?restfreq) set restfreq=115.2712
set jyperk=`calc "2*1.3807e7/((30/$restfreq)*206265)**2"`
#set bmaj=`gethd in=$source.cmnorm/bmaj format=arcsec`
#set bmin=`gethd in=$source.cmnorm/bmin format=arcsec`
set bmarea=`calc "1.133*$bmaj*$bmin"`

# Peak brightness map (convert to K)
rm -rf $source.pknorm
moment in=$source.cmnorm out=$source.pknorm mom=-2
rm -rf $source.tpeak
maths exp="<$source.pknorm>/($bmarea*$jyperk)" out=$source.tpeak
puthd in=$source.tpeak/bunit value='K' type=ascii
set crms=`gethd in=$source.cmnorm/rms`
set crms2=`calc "$crms/($bmarea*$jyperk)"`
echo "The new rms is $crms2"
puthd in=$source.tpeak/rms value=$crms2
rm -rf $source.sengain1
imsub in=$source.sengain out=$source.sengain1 region='image(1)'
rm -rf $source.fov
regrid in=$source.sengain1 out=$source.fov axes=1,2 tin=$source.tpeak
rm -rf $source.imgen
imgen in=$source.fov out=$source.imgen object=disk factor=0 \
    object=disk spar=1,0,0,70,70,0
rm -rf $source.tpeakmsk
maths exp=$source.tpeak out=$source.tpeakmsk \
    mask="<$source.imgen>.gt.0"
set maxplt=`histo in=$source.tpeakmsk | grep Max | awk '{print $3*1.}'`
cgdisp in=$source.tpeak,$source.tpeak,$source.fov,$source.imgen \
    type=p,c,c,c slev=a,$crms2,a,1,a,1 levs1=4,5,6,8,10 levs2=$msklev levs3=1 \
    region='arc,box(-60,-60,60,60)' device=$source.tpeak.ps/cps \
    options=full,solneg2,black,wedge beamtyp=b,l labtyp=arcsec nxy=1,1 \
    cols1=4 cols2=1 cols3=6 csize=1,0.7 range=$crms2,$maxplt,lin,-8 lines=2,2,2,2
if ($showplot != 'n') then
    gv $source.tpeak.ps &
endif
foreach type (tpeak fov)
    rm -f $source.$type.fits
    fits in=$source.$type out=$source.$type.fits op=xyout
end

# Velocity at peak brightness map
rm -rf $source.vt $source.vpeak
at_moment in=$source.cmnorm out=$source.vt mom=-3
maths exp="<$source.vt>" mask="<$source.pknorm>/$crms.gt.$clip1pk" \
    out=$source.vpeak
minmax in=$source.vpeak
set vlast=`calc "$vfirst+($nmaps-1)*$delv"`
set nlevs=`calc -i "($vlast-$vfirst)/$mom1del"`
awk -v v0=$vfirst -v dv=$mom1del -v nl=$nlevs \
    'BEGIN { for (i = 1; i < nl; i++) {printf "%d,",v0+i*dv} printf "%d\n",v0+nl*dv}' \
    > vlevs.txt
cgdisp in=$source.vpeak,$source.vpeak type=p,c nxy=1,1 \
    device=$source.vpeak.ps/cps region='arc,box(-60,-60,60,60)' \
    range=$vfirst,$vlast,lin,3 slev=a,1 levs1=@vlevs.txt \
    labtyp=arcsec options=black,wedge,full cols1=1 csize=1 lines=2,1
if ($showplot != 'n') then
    gv $source.vpeak.ps &
endif
rm -f $source.vpeak.fits
fits in=$source.vpeak out=$source.vpeak.fits op=xyout

# Generate the smoothed mask
rm -rf $source.snr
maths exp="$source.cm/$source.sen" out=$source.snr
rm -rf $source.smo
convol map=$source.snr fwhm=$fwhm pa=0 out=$source.smo options=final
echo "Smoothed $source.snr to fwhm=$fwhm" | tee -a $log
set nmin1 = `expr $nmaps - 1`
set exec="sigest in=$source.smo"
set srmsi=`$exec region="image(1,2)" | grep Est | awk '{print $4}'`
set srmsf=`$exec region="image($nmin1,$nmaps)" | grep Est | awk '{print $4}'`
set srms = `echo $srmsi $srmsf | awk '{if ($1<$2) {print $1} else {print $2}}'`
#set srms = `calc -f f7.5 "($srmsi+$srmsf)/2"`
set clip=`calc -f f10.5 "$srms*$clipsmo"`
echo "Masking $source.smo at $clipsmo sigma value: $clip" | tee -a $log
echo "" | tee -a $log

# Apply the mask
rm -rf $source.cm.mskd $source.cmmsk.mskd
maths exp="<$source.cm>" out=$source.cm.mskd mask="<$source.smo>.gt.$clip"
maths exp="<$source.cm.mskd>/(<$source.gain>*$bmarea*$jyperk)" \
    region="mask($source.nse)" out=$source.cmmsk.mskd
puthd in=$source.cmmsk.mskd/bunit value='K' type=ascii
rm -rf $source.cmmsk.unit
maths exp="<$source.cmmsk.mskd>*0+1" out=$source.cmmsk.unit
delhd in=$source.cmmsk.unit/mask

# Generate the masked moment map
rm -rf $source.mmom0
moment in=$source.cmmsk.mskd out=$source.mmom0 mom=0
rm -rf $source.mmom0.nch
moment in=$source.cmmsk.unit out=$source.mmom0.nch mom=0
rm -rf $source.knse1
maths exp="<$source.nse1>/($bmarea*$jyperk)" out=$source.knse1 \
    region="mask($source.nsesub)"
puthd in=$source.knse1/bunit value='K' type=ascii
# We only multiply by sqrt($delv) since .nch has an extra factor of $delv
rm -rf $source.mmom0.nse
maths exp="sqrt(<$source.mmom0.nch>*abs($delv))*<$source.knse1>" \
    out=$source.mmom0.nse

# Plot the mom-0 map
set maxval=`histo in=$source.mmom0 | grep Max | awk '{print $3}'`
set sencon=`histo in=$source.knse1 | grep Min | awk -v c=$msklev '{print 0.95*$3*c}'`
mv $source.mmom0/mask $source.mmom0/mask1
cgdisp in=$source.mmom0,$source.mmom0,$source.knse1 type=p,c,c nxy=1,1 \
    device=$source.mmom0.ps/cps range=0,$maxval,lin,-8 slev=a,$mom0lev,a,1 \
    levs1=1,4,9,16,25,36,49,81,100 levs2=$sencon labtyp=arcsec cols1=4 \
    options=black,wedge,full,solneg2 cols2=1 csize=1 lines=2,1 beamtyp=b,l
if ($showplot != 'n') then
    gv $source.mmom0.ps &
endif
mv $source.mmom0/mask1 $source.mmom0/mask
rm -f $source.mmom0.fits
fits in=$source.mmom0 out=$source.mmom0.fits op=xyout

mv mossdi2_iteration.log $source'_iteration.log'
echo "---------------------------------------------------------------------" \
    | tee -a $log
echo "" | tee -a $log
if ($purge == 'n') goto end

#--------------------------------------------------------------------------
# 5. Output to HTML log
#--------------------------------------------------------------------------

echo "-----HTML------"

if ($line == 'co') then
    rm -f imagelog.html
    echo "<tr>" >> imagelog.html
    echo "<td>$source:r</td>" >> imagelog.html
    echo "<td>$vsys</td>" >> imagelog.html
    echo -n "<td>$rawlis[1]:r" >> imagelog.html
    set srctim=(0)
    set srctim[1]=`uvindex vis=$filelis[1] | grep 'Total obs' | awk '{print $5}'`
    set tottim = $srctim[1]
    echo "Source time is $tottim"
    if ($#rawlis > 1) then
        set i = 2
        while ($i <= $#rawlis)
            echo -n "<br>$rawlis[$i]:r" >> imagelog.html
            set addtim=`uvindex vis=$filelis[$i] | grep 'Total obs' | awk '{print $5}'`
            echo "Adding additional source time of $addtim hr"
            set tottim = `calc -f f5.2 "$tottim+$addtim"`
            set i = `expr $i + 1`
        end
    endif
    echo "</td>" >> imagelog.html
    echo "Total time is $tottim"
    echo "<td>$tottim</td>" >> imagelog.html
    set nptng=`gethd in=$source.map/nptng`
    echo "<td>$nptng</td>" >> imagelog.html
    echo "<td>$bmaj x $bmin</td>" >> imagelog.html
    set jyrms=`gethd in=$source.cmnorm/rms`
    set mjyrms=`echo $jyrms | awk '{printf "%.1f",$1*1e3}'`
    echo "<td align="\"right\"">$mjyrms</td>" >> imagelog.html
    set trms=`gethd in=$source.tpeak/rms`
    set trmsint=`echo $trms | awk '{printf "%.1f",$1*1000.}'`
    set tmax=`histo in=$source.tpeakmsk | grep Max | awk '{printf "%.0f",$3*1000.}'`
    set snr=`histo in=$source.tpeakmsk | grep Max | awk -v s=$trms '{printf "%.2f",$3/s}'`
    echo "<td align="\"right\"">$trmsint</td>" >> imagelog.html
    echo "<td align="\"right\"">$tmax</td>" >> imagelog.html
    echo "<td align="\"right\"">$snr</td>" >> imagelog.html
    set flux=`histo in=$source.mmom0 | grep Sum | awk -v c=$cell '{printf "%4.3e",$6*c*c}'`
    echo "<td align="\"right\"">$flux</td>" >> imagelog.html
    if (-e $source.cmaps.ps) then
        ps2pdf $source.cmaps.ps
        mv $source.cmaps.pdf ../gal_cm
        echo "<td><a target="\"_NEW\"" href="\"http://mmwave.astro.illinois.edu/carma/edge/${arr}_gal_cm/$source.cmaps.pdf\"">link</a></td>" >> imagelog.html
    else
        echo "<td>ERROR</td>" >> imagelog.html
    endif
    foreach imtype (tpeak vpeak mmom0)
        if (-e $source.$imtype.ps) then
            epspdf $source.$imtype.ps
            #echo "$source:r" | enscript -B -f Courier-Bold@32 -o- | ps2pdf - | \
            #    pdftk $source.$imtype.pdf stamp - output output.pdf
            if ($imtype == 'tpeak') set dir = "gal_tpk"
            if ($imtype == 'vpeak') set dir = "gal_vpk"
            if ($imtype == 'mmom0') set dir = "gal_mom"
            #if (-e output.pdf) mv output.pdf $source.$imtype.pdf
            mv $source.$imtype.pdf ../$dir
            echo "<td><a target="\"_NEW\"" href="\"http://mmwave.astro.illinois.edu/carma/edge/${arr}_$dir/$source.$imtype.pdf\"">link</a></td>" >> imagelog.html
        else
            echo "<td>ERROR</td>" >> imagelog.html
        endif
    end
    echo "<td>`date`</td>" >> imagelog.html
    echo "</tr>" >> imagelog.html
endif

#--------------------------------------------------------------------------
# 6. Clean up files
#--------------------------------------------------------------------------

echo "-----dorm------"

dorm:
# Clean up:
rm -f tmp.* vlevs.txt
rm -rf $source.{cm.mskd,cmmsk.mskd,cmmsk.unit,imgen,mapnorm,nse,nse1,knse1,nsesub,pknorm,res,sengain,sengain1,senmsk,sennorm,smo,snr,vt}

# compress files:
gzip -f *.ps
if ($contsub == y && -e $source.uvlin && $savevis == y) then
    tar cvf $source.uvlin.tar $source.uvlin
endif
if ($savevis == y) then
    tar cvf $source.vis.tar $source.vis
endif
rm -rf $source.{vis,uvlin}

end: