#!/bin/csh -f
#-------------------------------------------------------------------
# General-purpose reduction script for ATCA data
# This script is not designed for multi-pointing or full Stokes data, 
#    nor does it allow flagging based on spectral channels (e.g. in 
#    cases of RFI).
# Requires: 
#    atlod   >23oct04
#    atfix   >27feb06 (tsys and gain-elev by default)
#    gpedit  >02jan05 (amp selection for gains)
#    uvflag  >26aug04 (edge=-1)
#    uvsplit >28jan05 (options=clobber)
#-------------------------------------------------------------------
set version=15sep15
set today=`date "+%y%m%d"`
set log=atcalib.$today.log

# Setup parameters for given dataset, edit this part

# Execution parameters
#set showplots=y
set interact=y
set purge=n
set xdev=/xw

# Datafiles: give wildcards in quotes
set datadir=
set files=

# Atlod options
set atlodopt=

# Atfix options (for 3mm only)
set dotsys=n         # y to interpolate paddle-derived Tsys
set doelev=n         # y to apply gain-elevation correction
# array configuration name or positions (3-element 3mm data)
set array=
# baseline correction in ns, e.g. set dantpos=@$MIRCAT/dantpos.030910
set dantpos=

# Sources (if selfcal, leave sources blank and give as gaincal)
set sources=
set gaincals=
set testcals=
set bpcal=
set fluxcal=

# Frequencies in MHz (if calibrating phase on just 1 of 2 freqs, make it freq2)
set freq1=
set freq2=
# To use chan 1 freq instead of centre freq to label windows, set this to y:
set ch1split=n

# Rest freqs in GHz (leave blank for continuum)
set restfq1=
set restfq2=

# Usually leave these parameters unset.  By default, if a rest frequency 
# is specified the corresponding spectral window is mapped as a datacube.
# Frequencies to map as velocity datacubes:
set velmapfq=
# Frequencies to map as continuum images:
set conmapfq=

# Flagging parameters
set fluxtime=         # optional time range for flux calibrator
set edge1=-1          # edge chans to flag for f1 (default: nch/8 on each end)
set edge2=-1          # edge chans to flag for f2 (default: nch/8 on each end)
set shadow=           # shadowing limit in metres
set flagsrc=          # optional source(s) for extra flagging
set flagsel=          # optional selection(s) for extra flagging

# Calibration strategy
set transfer=n        # Calibrate freq1 using freq2, assuming constant offset
set phaseonly=n       # y to calibrate only phase gains
# Optional chans to average before calibration, if not all (e.g. maser):
# e.g. set calchan='line=chan,10,40,1,1'
set calchan=

# General calibration parameters
set refant=3          # reference antenna for phase
set bpbin1=           # no. of channels to bin in F1 during bandpass cal
set bpbin2=           # no. of channels to bin in F2 during bandpass cal
set bpint=1           # time interval (mins) for bandpass calibration
set gainint=5,2       # time interval (mins) for gain calibration
set gainclip=         # clip amp gains above this level

# Optional flux calibration parameters (quasar bootstrapping)
set calflux=          # gain calibrator's flux, usually unset
set fluxflux=         # flux calibrator's flux, if a qsr not known by mfcal
# time/antenna selection for gaincal for gpboot
# e.g. set bootsel='select=time(10:00,10:30)'
set bootsel=

# Optional flux calibration parameters (planet bootstrapping)
# interval to average ampgains over before applying to planet
set fluxavint=
# optional time range where cal and planet have similar elev/Tsys
set gainsel=          
# antenna selection for planet for plboot
# e.g. set plselect='select=-ant(2)(4)'
set plselect=

# Mapping parameters
set chbin=2           # binning factor for velocity channels (default: 2)
set imline=           # linetype for invert to override vel linetype
set impol='xx,yy'     # Polarisation(s) to image
set cellsz=           # pixel size in arcsec (default: 1/3 of resolution)
set imagsz=           # image size in pixels (default: to HPBW)
set no6ant=n          # discard CA06 when mapping?
set contsub=n         # subtract continuum before imaging line?
set uvlinch=          # ranges of line-free channels (e.g. 1,50,80,128)
set fitorder=1        # polynomial order fit for continuum subtraction
set rmsfac=3          # CLEAN cutoff in units of map RMS
set maxsen=2.5        # sensitivity contour level for deconvolution
set msklev=2          # sensitivity contour level for signal detection

#######################################################################

alias mkpdf 'gs -dBATCH -dNOPAUSE -sPAPERSIZE=letter -dPDFFitPage -q -sDEVICE=pdfwrite -sOutputFile=\!:1.pdf \!:1.ps'

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

# Get lists of datafiles, sources and calibrators
set files = `echo "$files"     | sed 's/,/ /g'`
set sources = `echo "$sources" | sed 's/,/ /g'`
set gaincals=`echo "$gaincals" | sed 's/,/ /g'`
set testcals=`echo "$testcals" | sed 's/,/ /g'`
set velmapfq=`echo "$velmapfq" | sed 's/,/ /g'`
set conmapfq=`echo "$conmapfq" | sed 's/,/ /g'`
set flagsrc=`echo "$flagsrc" | sed 's/,/ /g'`
set flagsel=`echo "$flagsel" | sed 's/+/ /g'`

# Figure out some parameters
set allcals = ()
foreach obj ($gaincals $testcals $bpcal $fluxcal)
    set val=`echo $allcals | grep $obj`
    if ("$val" == "") set allcals = ($obj $allcals)
end
echo "Calibrators: $allcals"
if ("$velmapfq" == "" && "$conmapfq" == "") then
    if ($restfq1 != "") then
        set velmapfq = ($freq1)
    else
        set conmapfq = ($freq1)
    endif
    if ($restfq2 != "") then
        set velmapfq = ($velmapfq $freq2)
    else
        set conmapfq = ($conmapfq $freq2)
    endif
endif

# Output some interesting information
echo ""
echo "atcalib.csh $version executed by $user on `date`" | tee $log
echo "" | tee -a $log
set pdflist=()

if ($interact == 'n') then
    goto load
endif

# MAIN MENU
echo ""
echo "ATCA Miriad processing script version $version"
menu:
echo ""
echo "Where would you like to start?"
echo "  0. EXIT" 
echo "  1. atlod & atfix"
echo "  2. split by source and frequency"
echo "  3. non-interactive data flagging (uvflag)"
echo "  4. interactive data flagging (blflag)"
echo "  5. bandpass calibration"
echo "  6. determine gain offsets between freqs"
echo "  7. gains vs. time calibration"
echo "  8. bootstrap flux from planet or primary calibrator"
echo "  9. copy final gains vs. time to source"
echo " 10. map the sources (spectral line)"
echo " 11. map the sources (continuum)"
echo " 12. purge non-essential files (when done)"
echo ""
echo -n "? "
set start = $<
switch ($start)
    case 0:
    goto end
    case 1:
    goto load
    case 2: 
    goto split
    case 3: 
    goto uvflagger
    case 4: 
    goto blflagger
    case 5: 
    goto bpcal
    case 6: 
    goto goffset
    case 7: 
    goto gcal
    case 8: 
    goto bootflux
    case 9: 
    goto gapply
    case 10: 
    goto mapline
    case 11: 
    goto mapcont
    case 12: 
    goto purge
endsw
goto menu

load:
#######################################################################
# 1. Read in data and apply Tsys and gain-elevation correction
#######################################################################

# Read in data with atlod
set i = 1
foreach file ($files)
    echo "Working on file $datadir/$file"
    if !(-e atfile.$i) then
        rm -rf atfile.$i
        atlod in=$datadir/$file out=atfile.$i options=$atlodopt
    else
        echo "NOTE: Using existing file atfile.$i"
    endif
    if !(-e atfile.$i) then
        echo "Running ATLOD a second time with options=noif..."
        atlod in=$datadir/$file out=atfile.$i options=noif,$atlodopt
    endif
    set i = `expr $i + 1`
end

# atfix for 3mm observations
if ($dotsys == 'y' || $doelev == 'y' || $dantpos != '' || $array != '') then
    set fixopts=
    if ($dantpos != "") set fixopts="$fixopts dantpos=$dantpos"
    if ($array != "") set fixopts="$fixopts array=$array"
    if ($doelev == "n") set fixopts="$fixopts options=nogainel"
    if ($dotsys == "n") set fixopts="$fixopts tsyscal=any"
    foreach file (atfile.*)
        set seq=$file:e
        echo "Running atfix with $fixopts"
        rm -rf atfix.$seq
        atfix vis=$file out=atfix.$seq $fixopts
    end
endif

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


split:
#######################################################################
# 2. Split data by source and frequency
#######################################################################

# Split data by sources:
if (-e atfix.1) then
    set dofile=atfix
else
    set dofile=atfile
endif

# Accept only data in proper time ranges:
foreach obj ($sources $allcals)
    if (-e $obj) rm -rf $obj
    if (($obj == $fluxcal) && ($fluxtime != "")) then
        echo "Restricting $fluxcal to time $fluxtime"
        uvcat vis="$dofile.*" out=$obj options=unflagged \
            select="time($fluxtime),source($obj)"
    else if (`echo $obj | grep mos` != "") then
        set base=`echo $obj | sed 's/mos//g'`
        echo "Extracting data for mosaic source $obj"
        uvcat vis="$dofile.*" out=$obj options=unflagged \
            select="source($base*)"     
    else
        echo "Extracting data for source $obj"
        uvcat vis="$dofile.*" out=$obj options=unflagged \
            select="source($obj)" 
    endif
end

# Plot UV coverage:
if ($#sources != 0) then
    uvplt vis="$sources" axis=uc,vc device=uvcov.ps/cps \
        options=nobase,equal select='pol(xx)' inc=17
else if ($#gaincals != 0) then
    uvplt vis="$gaincals[1]" axis=uc,vc device=uvcov.ps/cps \
        options=nobase,equal select='pol(xx)' inc=17
endif
mkpdf uvcov
echo "$files" | enscript -r -s 20 -B -f Courier-Bold@12 -o- | \
    ps2pdf - | pdftk uvcov.pdf stamp - output output.pdf
if (-e output.pdf) mv output.pdf uvcov.pdf
set pdflist=($pdflist uvcov.pdf)

# Split data by frequencies:
foreach obj ($sources $allcals)
    if ($ch1split == 'y') then
        set ch1fqs = `uvindex vis=$obj | grep GHz | awk '{print $2}'`
        foreach ch1 ($ch1fqs)
            set fq=`calc -i "$ch1*1000"`
            set f1=`calc -f f8.4 "$ch1-0.001"`
            set f2=`calc -f f8.4 "$ch1+0.001"`
            rm -rf $obj.$fq
            uvcat vis=$obj select="freq($f1,$f2)" out=$obj.$fq
        end
    else if (`echo $obj | grep mos` != "") then
        uvsplit vis=$obj options=clobber,mosaic
        set base=`echo $obj | awk -F_ '{print $1}'`
        rm -rf $obj.[0-9]*[0-9]
        foreach file ($base.*)
            echo "Renaming $file to $obj.$file:e"
            mv $file $obj.$file:e
        end
    else
        uvsplit vis=$obj options=clobber
    endif
    if ($restfq1 != "") puthd in=$obj.$freq1/restfreq value=$restfq1
    if ($restfq2 != "") puthd in=$obj.$freq2/restfreq value=$restfq2
end

if ($interact != 'n') then
    echo ""
    echo -n "Continue with initial data flagging [y]? "
    set resp = $<
    if ($resp == 'n' || $resp == 'N') goto end
endif


uvflagger:
#######################################################################
# 3. Automated flagging of bad data
#######################################################################

# Flagging of edge channels and shadowed data:
foreach obj ($sources $allcals)
    if ($edge1 != "") then
        echo "Flagging edge channels $edge1 on $obj.$freq1"
        uvflag vis=$obj.$freq1 edge=$edge1 flagval=f > /dev/null
    endif
    if ($edge2 != "" && $freq2 != "") then
        echo "Flagging edge channels $edge2 on $obj.$freq2"
        uvflag vis=$obj.$freq2 edge=$edge2 flagval=f > /dev/null
    endif
    if ($shadow != "") then
        foreach freq ($freq1 $freq2)
            echo "Flagging data shadowed up to $shadow m on $obj.$freq"
            uvflag vis=$obj.$freq select="shad($shadow)" flagval=f > /dev/null
        end
    endif
end

# Flagging from previous runs of blflag:
foreach obj ($sources $allcals)
    foreach freq ($freq1 $freq2)
        if (-e $obj.$freq.flgselect) then
            echo ""
            echo "uvflag vis=$obj.$freq select=@$obj.$freq.flgselect"
            uvflag vis=$obj.$freq select=@$obj.$freq.flgselect flagval=f 
        endif
    end
end

# Flagging specified in preamble:
set i = $#flagsrc
echo ""
echo "There are $i flag sources"
if ($#flagsrc > $#flagsel) then
    set j=$#flagsel
    while ($j < $i)
        set flagsel = ($flagsel $flagsel[$j])
        set j=`expr $j + 1`
    end
endif
echo "Flag sources: $flagsrc"
echo "Flag selection: $flagsel"
while ($i > 0)
    if ($flagsrc[$i] != "") then
        echo ""
        echo "Flagging data for $flagsrc[$i] select=$flagsel[$i]"
        foreach freq ($freq1 $freq2)
            uvflag vis=$flagsrc[$i].$freq select="$flagsel[$i]" flagval=f
        end
    endif
    set i = `expr $i - 1`
end

# Plot amp vs. time for calibrators:
# foreach obj ($allcals)
#     foreach freq ($freq1 $freq2)
#         uvplt vis="$obj.$freq" axis=time,ampl device=$obj.$freq.tamp.ps/cps \
#             options=nobase stokes=xx,yy size=1,4 > /dev/null
#         if ($showplots == "y") then
#             gv $obj.$freq.tamp.ps &
#         endif
#     end
# end

if ($interact == 'n') then
    goto bpcal
else
    echo ""
    echo -n "Continue with interactive data flagging [n]? "
    set resp = $<
    if ($resp != 'y' && $resp != 'Y') goto bpcal
endif


blflagger:
#######################################################################
# 4. Interactive flagging of bad data
#######################################################################

# Use blflag to flag points in time
foreach obj ($sources $allcals)
    foreach freq ($freq1 $freq2)
        echo ""
        echo "Displaying data for $obj.$freq"
        blflag vis=$obj.$freq stokes=xx,yy options=selgen device=/xs
        if (-e blflag.select) mv blflag.select $obj.$freq.flgselect
    end
end

if ($bpcal == "") goto endbpcal
if ($interact != 'n') then
    echo ""
    echo -n "Continue with bandpass calibration [y]? "
    set resp = $<
    if ($resp == 'n' || $resp == 'N') goto menu
endif


bpcal:
#######################################################################
# 5. Calibrate the bandpass
#######################################################################

if ($bpcal == "") then
    echo 'ERROR: No bandpass calibrator selected\!'
    goto endbpcal
endif

# Examine bandpass calibrator
foreach freq ($freq1 $freq2)
    uvplt vis=$bpcal.$freq axis=time,amp device=$bpcal.$freq.tamp.ps/cps \
        options=nocal,nopass stokes=xx,yy nxy=4,3 size=1.5,6
    uvplt vis=$bpcal.$freq axis=time,phase device=$bpcal.$freq.tph.ps/cps \
        options=nocal,nopass stokes=xx,yy nxy=4,3 size=1.5,6 yrange=-180,180
    uvspec vis=$bpcal.$freq axis=chan,phase device=$bpcal.$freq.pbp.ps/cps \
        interval=1000 nxy=4,3 stokes=xx,yy yrange=-180,180 \
        options=nocal,nopass
    uvspec vis=$bpcal.$freq axis=chan,amp device=$bpcal.$freq.pba.ps/cps \
        interval=1000 nxy=4,3 stokes=xx,yy options=nocal,nopass,ampsca
    foreach file (tamp tph pba pbp)
        mkpdf $bpcal.$freq.$file
        echo $bpcal.$freq.$file.pdf | enscript -r -s 400 -B -i 5.5i -f Courier-Bold@18 \
            -o- | ps2pdf - | pdftk $bpcal.$freq.$file.pdf stamp - output output.pdf
        if (-e output.pdf) mv output.pdf $bpcal.$freq.$file.pdf
        set pdflist=($pdflist $bpcal.$freq.$file.pdf)
    end
end

if ($interact != 'n') then
    echo ""
    echo -n "Proceed with bandpass calibration [y]? "
    set resp = $<
    if ($resp == 'n' || $resp == 'N') goto endbpcal
endif

# Calibrate the bandpass
foreach freq ($freq1 $freq2)
    if ($freq == $freq1 && $bpbin1 != "") then
        uvlist vis=$bpcal.$freq1 options=spectra log=uvlist.tmp
        set totchan=`grep Number uvlist.tmp | awk '{print $5}'`
        set binchan=`calc -i "($totchan-1)/$bpbin1"`
        set bpextra="line=chan,$binchan,1,$bpbin1,$bpbin1"
    else if ($freq != $freq1 && $bpbin2 != "") then
        uvlist vis=$bpcal.$freq2 options=spectra log=uvlist.tmp
        set totchan=`grep Number uvlist.tmp | awk '{print $5}'`
        set binchan=`calc -i "($totchan-1)/$bpbin2"`
        set bpextra="line=chan,$binchan,1,$bpbin2,$bpbin2"
    else
        set bpextra=
    endif
    if ($fluxcal != "") then
        if ($bpcal == $fluxcal && $fluxflux != "") then
            mfcal vis=$bpcal.$freq interval=$bpint refant=$refant $bpextra \
                flux=$fluxflux
        else
            mfcal vis=$bpcal.$freq interval=$bpint refant=$refant $bpextra
        endif
    else
        mfcal vis=$bpcal.$freq interval=$bpint refant=$refant $bpextra
    endif
    gpplt vis=$bpcal.$freq yaxis=phase device=bpgains.$freq.pbp.ps/ps \
        nxy=3,2 options=bandpass,wrap yrange=-180,180
    gpplt vis=$bpcal.$freq yaxis=amp device=bpgains.$freq.pba.ps/ps \
        nxy=3,2 options=bandpass yrange=0,2
    foreach file (pba pbp)
        mkpdf bpgains.$freq.$file
        echo bpgains.$freq.$file | enscript -r -s 400 -B -i 7i -f Courier-Bold@18 \
            -o- | ps2pdf - | pdftk bpgains.$freq.$file.pdf stamp - output output.pdf
        if (-e output.pdf) mv output.pdf bpgains.$freq.$file.pdf
        set pdflist=($pdflist bpgains.$freq.$file.pdf)
    end
end

set bpquery=n
if ($interact != 'n') then
    echo ""
    echo -n "Apply to all sources and calibrators [y]? "
    set resp = $<
    if ($resp == 'n' || $resp == 'N') set bpquery=y
endif

# If desired, apply the bandpass table:
foreach obj ($sources $allcals)
    if ($obj != $bpcal && $bpquery == 'y') then
        echo -n "Apply to source $obj [n]? "
        set bpresp = $<
    else
        set bpresp = 'y'
    endif
    foreach freq ($freq1 $freq2)
        if ($obj != $bpcal && $bpresp == 'y') then
            echo "Copying to $obj.$freq"
            gpcopy vis=$bpcal.$freq out=$obj.$freq options=nocal
        endif
    end
end

endbpcal:
echo ""
if ($interact == 'n') then
    if ($transfer != "y") goto gcal
else
    if ($transfer == "y") then
        echo -n "Estimate gain offset between 2 freqs [y]? "
        set resp = $<
        if ($resp == "n" || $resp == 'N') goto menu
    else
        goto endgoff
    endif
endif


goffset:
#######################################################################
# 6. Determine the gain offsets between the 2 frequencies
#######################################################################

if ($transfer == "n") then
    echo 'ERROR: you must set transfer=y to perform this step\!'
    goto end
endif
if ($bpcal == "") then
    echo 'ERROR: a bandpass calibrator is required\!'
    goto end
endif

# Make a second copy of the data set:
if !(-e RELCAL) mkdir RELCAL
cd RELCAL
foreach freq ($freq1 $freq2)
    rm -rf $bpcal.$freq.relcal
    cp -r ../$bpcal.$freq $bpcal.$freq.relcal
end

# Calibrate $freq2 data and copy gain table to $freq1:
mfcal vis=$bpcal.$freq2.relcal interval=0.1 refant=$refant flux=1
gpcopy vis=$bpcal.$freq2.relcal out=$bpcal.$freq1.relcal options=nopass
uvplt vis=$bpcal.$freq1.relcal axis=time,phase device=relcal.ph.ps/cps \
    yrange=-180,180 nxy=3,2 size=2,6 select='pol(xx,yy)'
uvplt vis=$bpcal.$freq1.relcal axis=time,amp device=relcal.amp.ps/cps \
    yrange=0,2 nxy=3,2 size=2,6 select='pol(xx,yy)'
# if ($showplots == "y") then
#     gv relcal.ph.ps &
#     gv relcal.amp.ps &
# endif

if ($interact != 'n') then
    echo ""
    echo -n "Proceed to determine offsets [y]? "
    set resp = $<
    if ($resp == 'n' || $resp == 'N') goto end
endif

# Re-calibrate data on $freq1:
rm -rf $bpcal.$freq1.relcal2
uvcat vis=$bpcal.$freq1.relcal out=$bpcal.$freq1.relcal2

# Determine the number of polns
set npol = `gethd in=$bpcal.$freq1.relcal2/npol`
echo "I find $npol polarisation(s) in this dataset"
if ($npol == 1) then
    set dopol = (x)
else
    set dopol = (x y)
endif

# Calculate the median amp and phase gains:
foreach pol ($dopol)
    mfcal vis=$bpcal.$freq1.relcal2 interval=0.1 refant=$refant flux=1 \
        options=nopass select="pol($pol$pol)"
    if ($interact == 'y') then
        set gpdev=$xdev
    else
        set gpdev=/null
    endif
    gpplt vis=$bpcal.$freq1.relcal2 yaxis=phase device=$gpdev \
        nxy=3,2 log=gpplt.ph.tmp
    gpplt vis=$bpcal.$freq1.relcal2 yaxis=amp device=$gpdev \
        nxy=3,2 log=gpplt.amp.tmp
    foreach type (amp ph)
        grep -v '#' gpplt.$type.tmp | awk '{print $3}' > a1$pol.$type.out
        grep -v '#' gpplt.$type.tmp | awk '{print $4}' > a2$pol.$type.out
        grep -v '#' gpplt.$type.tmp | awk '{print $5}' > a3$pol.$type.out
        grep -v '#' gpplt.$type.tmp | awk '{print $6}' > a4$pol.$type.out
        grep -v '#' gpplt.$type.tmp | awk '{print $7}' > a5$pol.$type.out
        grep -v '#' gpplt.$type.tmp | awk '{print $8}' > a6$pol.$type.out
        foreach ant (1 2 3 4 5 6)
            cat a$ant$pol.$type.out | sort -n > sort.tmp
            set nlines=`wc -l sort.tmp | awk '{print $1}'`
            set midline=`expr \( $nlines + 1 \) / 2`
            set median = `head -n $midline sort.tmp | tail -1`
            echo "The median $type gain for CA0$ant $pol is $median"
            echo $median > a$ant$pol.$type.med
        end
    end
end

# Apply back to the bandpass cal to check:
foreach ant (1 2 3 4 5 6)
    foreach pol ($dopol)
        set gamp = `cat a"$ant"$pol.amp.med`
        set gph = `cat a"$ant"$pol.ph.med`
        if ($gamp != '0.000') then
            echo "Ant $ant Pol $pol - Replacing gains by $gamp $gph"
            gpedit vis=$bpcal.$freq1.relcal select="ant($ant)" feeds=$pol \
                gain=$gamp,$gph options=multiply
        endif
    end
end
uvplt vis=$bpcal.$freq1.relcal axis=time,phase device=post.ph.ps/vcps \
    nxy=1,3 size=2,6 select='pol(xx,yy)' yrange=-180,180
uvplt vis=$bpcal.$freq1.relcal axis=time,amp device=post.amp.ps/vcps \
    nxy=1,3 size=2,6 select='pol(xx,yy)' yrange=0,2
# if ($showplots == "y") then
#     gv post.ph.ps &
#     gv post.amp.ps &
# endif
rm -f *.tmp
cd ..

endgoff:
if ($interact != 'n') then
    echo ""
    echo -n "Continue with gains vs. time calibration [y]? "
    set resp = $<
    if ($resp == 'n' || $resp == 'N') goto menu
endif


gcal:
#######################################################################
# 7. Calibrate gains vs. time
#######################################################################

# Select the frequencies:
if ($transfer == "n") then
    set dofreq = ($freq1 $freq2)
else
    set dofreq = ($freq2)
endif

foreach obj ($gaincals)
    foreach freq ($dofreq)
        # ----- Apply the bandpass
        rm -rf $obj.$freq.cal
        uvcat vis=$obj.$freq out=$obj.$freq.cal options=unflagged,nocal
        # ----- Plot amps and phases before calibration
        uvplt vis=$obj.$freq.cal axis=time,amp device=$obj.$freq.tamp.ps/cps \
            options=nobase stokes=xx,yy size=1,4
        mkpdf $obj.$freq.tamp
        set pdflist=($pdflist $obj.$freq.tamp.pdf)
        if ($interact == 'y') then
            uvplt vis=$obj.$freq.cal axis=time,phase device=$xdev $calchan \
                select='pol(xx,yy)' options=source size=2,3 yrange=-180,180
        endif
        # ----- Do the calibration
        if ($calflux != "") then
            mfcal vis=$obj.$freq.cal interval=$gainint refant=$refant \
                $calchan options=nopassol flux=$calflux
        else
            mfcal vis=$obj.$freq.cal interval=$gainint refant=$refant \
                $calchan options=nopassol
        endif
        # ----- Plot amps and phases after calibration
        if ($interact == 'y') then
            uvplt vis=$obj.$freq.cal axis=time,phase device=$xdev $calchan \
                select='pol(xx,yy)' options=source size=2,3 yrange=-180,180
            uvplt vis=$obj.$freq.cal axis=time,amp device=$xdev $calchan \
                select='pol(xx,yy)' options=source size=2,3
        endif
    end
end

# Consolidate the gain tables into first gaincal and plot:
foreach freq ($dofreq)
    foreach obj ($gaincals)
        if ($obj != $gaincals[1]) then
            gpcopy vis=$obj.$freq.cal out=$gaincals[1].$freq.cal \
                mode=merge options=nopass
        endif
    end
    if ($gainclip != "") then
        echo "Clipping gains with amplitudes above $gainclip..."
        gpedit vis=$gaincals[1].$freq.cal select="-amp(0,$gainclip)" \
            options=flag
        gpplt vis=$gaincals[1].$freq.cal yaxis=amp nxy=3,2 \
            device=gains.$freq.amp.ps/ps yrange=0,$gainclip 
    else
        gpplt vis=$gaincals[1].$freq.cal yaxis=amp nxy=3,2 \
            device=gains.$freq.amp.ps/ps
    endif
    gpplt vis=$gaincals[1].$freq.cal yaxis=phase nxy=3,2 \
        device=gains.$freq.ph.ps/ps options=wrap yrange=-180,180
    foreach file (amp ph)
        mkpdf gains.$freq.$file
        echo gains.$freq.$file | enscript -r -s 400 -B -i 7i -f Courier-Bold@18 \
            -o- | ps2pdf - | pdftk gains.$freq.$file.pdf stamp - output output.pdf
        if (-e output.pdf) mv output.pdf gains.$freq.$file.pdf
        set pdflist=($pdflist gains.$freq.$file.pdf)
    end
#     if ($showplots == "y") then
#         gv gains.$freq.amp.ps &
#         gv gains.$freq.ph.ps &
#     endif
end

# If transferring to $freq1:
if ($transfer == "y") then
    set npol = `gethd in=$gaincals[1].$freq2.cal/npol`
    echo "I find $npol polarisation(s) in this dataset"
    if ($npol == 1) then
        set dopol = (x)
    else
        set dopol = (x y)
    endif
    rm -rf $gaincals[1].$freq1.cal
    uvcat vis=$gaincals[1].$freq1 out=$gaincals[1].$freq1.cal options=unflagged
    gpcopy vis=$gaincals[1].$freq2.cal out=$gaincals[1].$freq1.cal \
        options=nopass
    foreach ant (1 2 3 4 5 6)
        foreach pol ($dopol)
            if (-e RELCAL/a$ant$pol.amp.med) then
                set gamp = `cat RELCAL/a"$ant"$pol.amp.med`
                set gph = `cat RELCAL/a"$ant"$pol.ph.med`
                echo "Ant $ant Pol $pol - Replacing gains by $gamp $gph"
                gpedit vis=$gaincals[1].$freq1.cal select="ant($ant)" \
                    feeds=$pol gain=$gamp,$gph options=multiply
            endif
        end
    end
    gpplt vis=$gaincals[1].$freq1.cal yaxis=phase nxy=3,2 \
        device=gains.$freq1.ph.ps/ps options=wrap yrange=-180,180
    gpplt vis=$gaincals[1].$freq1.cal yaxis=amp nxy=3,2 \
        device=gains.$freq1.amp.ps/ps
#     if ($showplots == "y") then
#         gv gains.$freq1.ph.ps &
#         gv gains.$freq1.amp.ps &
#     endif
    # ----- Verify on main phase calibrator or selfcal source
    if ($interact == 'y') then
        uvplt vis=$gaincals[1].$freq1.cal axis=time,phase device=$xdev \
            $calchan select='pol(xx,yy)' 
        uvplt vis=$gaincals[1].$freq1.cal axis=time,amp device=$xdev \
            $calchan select='pol(xx,yy)'
    endif
endif

if ($interact == 'n') then
    if ($fluxcal == "") goto gapply
else
    if ($fluxcal == "") goto endboot
    echo ""
    echo -n "Continue with absolute flux calibration [y]? "
    set resp = $<
    if ($resp == 'n' || $resp == 'N') goto menu
endif


bootflux:
#######################################################################
# 8. Determine calibrator fluxes and rescale amp gains
#######################################################################

if ($fluxcal == "") then
    echo "For this step you must specify a flux calibrator\!"
    goto menu
endif

set planets='mercury venus mars jupiter saturn uranus neptune'
if (`echo $planets | grep $fluxcal` == "") then
    set planetcal = 'n'
else
    set planetcal = 'y'
endif

foreach freq ($freq1 $freq2)
    if ($fluxcal != $bpcal) then
# Examine flux calibrator
        uvplt vis=$fluxcal.$freq axis=time,amp device=$fluxcal.$freq.tamp.ps/cps \
            nxy=4,3 size=1.5,6 select='pol(xx,yy)' options=nocal
        mkpdf $fluxcal.$freq.tamp
        echo "$fluxcal, pb cal only" | enscript -r -s 400 -B -i 6i -f Courier-Bold@18 \
            -o- | ps2pdf - | pdftk $fluxcal.$freq.tamp.pdf stamp - output output.pdf
        if (-e output.pdf) mv output.pdf $fluxcal.$freq.tamp.pdf
        set pdflist=($pdflist $fluxcal.$freq.tamp.pdf)
#         if ($showplots == "y") then
#             gv $fluxcal.$freq.tamp.ps &
#         endif
# If flux cal is a quasar, first calibrate the data:
        if ($planetcal == 'n' && $fluxflux == "") then
            mfcal vis=$fluxcal.$freq interval=$gainint refant=$refant \
                options=nopass
        else if ($planetcal == 'n' && $fluxflux != "") then
            mfcal vis=$fluxcal.$freq interval=$gainint refant=$refant \
                options=nopass flux=$fluxflux
# If flux cal is a planet, vector avg in time and copy ampgains from gain cal:
        else if ($planetcal == 'y') then
            rm -rf $fluxcal.$freq.ave
            uvaver vis=$fluxcal.$freq out=$fluxcal.$freq.ave \
                interval=$gainint options=nocal
            gpcopy vis=$gaincals[1].$freq.cal out=$fluxcal.$freq.ave \
                options=nopass
            # --- Gain fiddling: select a time range where calibrator is 
            # --- close to planet's elevation and average gains in time
            if ($fluxavint != "") then
                if ($gainsel != "") then
                    gpedit vis=$fluxcal.$freq.ave select="-time($gainsel)" \
                        options=flag
                    echo "Selecting only gains in time range $gainsel"
                endif
                echo "Averaging gains over $fluxavint min"
                gpaver vis=$fluxcal.$freq.ave interval=$fluxavint \
                    options=scalar
                puthd in=$fluxcal.$freq.ave/interval value=1.
                gpplt vis=$fluxcal.$freq.ave yaxis=amp nxy=3,2 \
                    device=$fluxcal.$freq.gains.ps/ps 
                mkpdf $fluxcal.$freq.gains
                set pdflist=($pdflist $fluxcal.$freq.gains.pdf)
#                 if ($showplots == "y") then
#                     gv $fluxcal.$freq.gains.ps &
#                 endif
            endif
        endif
    endif
end

if ($interact != 'n') then
    echo ""
    echo -n "Proceed with flux calibration [y]? "
    set resp = $<
    if ($resp == 'n' || $resp == 'N') goto end
endif

# Select the frequencies to work on:
set bootfreq = ($freq1 $freq2)

# Bootstrap the gains:
foreach freq ($bootfreq)

    set doflux = ($fluxcal.$freq.ave $gaincals[1].$freq.cal)
    if ($bpcal != "" && $bpcal != $fluxcal && $bpcal != $gaincals[1]) then
        rm -rf $bpcal.$freq.cal
        uvcat vis=$bpcal.$freq out=$bpcal.$freq.cal
        set doflux = ($doflux $bpcal.$freq.cal)
    endif

    if ($#bootfreq == 1) then
        set doboot = ($doflux $gaincals[1].$freq1.cal)
    else
        set doboot = ($doflux)
    endif
    echo "$doboot"

    if ($planetcal == y) then
        plboot vis="$doboot" $plselect device=plboot.$freq.ps/cps \
            | tee -a $log
        echo " " | tee -a $log
#         mkpdf plboot.$freq
#         set pdflist=($pdflist plboot.$freq.pdf)
#         if ($showplots == "y") then
#             gv plboot.$freq.ps &
#         endif
        uvflux vis="$doflux" $plselect
    else
        gpboot vis=$gaincals[1].$freq.cal cal=$fluxcal.$freq $bootsel \
            | tee -a $log
        uvflux vis=$gaincals[1].$freq.cal,$fluxcal.$freq
    endif
    
end

endboot:
if ($interact != 'n') then
    echo ""
    echo -n "Copy the gain tables to the sources [y]? "
    set resp = $<
    if ($resp == 'n' || $resp == 'N') goto menu
endif


gapply:
#######################################################################
# 9. Copy gains vs. time to sources
#######################################################################

if ($#sources == 0) then
    echo 'ERROR: selfcal project, cannot copy gains\!'
    if ($interact != 'n') then
        echo -n "Proceed to map $gaincals[1] [y]? "
        set resp = $<
        if ($resp == 'n' || $resp == 'N') goto menu
    endif
    goto mapline
endif

# Copy calibration data to source(s):
foreach source ($sources $testcals)
    foreach freq ($freq1 $freq2)
        gpcopy vis=$gaincals[1].$freq.cal out=$source.$freq options=nopass
        if ($phaseonly == 'y') then
            gpedit vis=$source.$freq options=phase
        endif
        rm -rf $source.$freq.cal
        uvcat vis=$source.$freq out=$source.$freq.cal options=unflagged
    end
end

# Plot spectra
foreach source ($sources)
    foreach freq ($freq1 $freq2)
        if (`echo $velmapfq | grep $freq` == "") then
            uvspec vis=$source.$freq.cal device=$source.$freq.spec.ps/cps \
                axis=freq,amp interval=10000 nxy=3,4 select="pol($impol)"
        else
            uvspec vis=$source.$freq.cal device=$source.$freq.spec.ps/cps \
                axis=vel,amp interval=10000 nxy=3,4 select="pol($impol)"
        endif    
        mkpdf $source.$freq.spec
        echo $source.$freq.spec | enscript -r -s 420 -B -i 6i -f Courier-Bold@18 \
            -o- | ps2pdf - | pdftk $source.$freq.spec.pdf stamp - output output.pdf
        if (-e output.pdf) mv output.pdf $source.$freq.spec.pdf
        set pdflist=($pdflist $source.$freq.spec.pdf)
#         if ($showplots == "y") then
#             gv $source.$freq.spec.ps &
#         endif
    end
end
foreach source ($testcals)
    foreach freq ($freq1 $freq2)
        uvspec vis=$source.$freq.cal device=$source.$freq.spec.ps/cps \
            axis=freq,phase interval=10000 nxy=3,4 select="pol($impol)" \
            yrange=-180,180
        mkpdf $source.$freq.spec
        uvplt vis=$source.$freq.cal device=$source.$freq.xph.ps/cps \
            axis=time,phase size=1.5,7 nxy=3,4 select="pol($impol)" \
            yrange=-180,180
        mkpdf $source.$freq.xph
        set pdflist=($pdflist $source.$freq.spec.pdf $source.$freq.xph.pdf)
#         if ($showplots == "y") then
#             gv $source.$freq.spec.ps &
#             gv $source.$freq.xph.ps &
#         endif
    end
end

if ($interact != 'n') then
    echo ""
    echo -n "Proceed to map the sources [y]? "
    set resp = $<
    if ($resp == 'n' || $resp == 'N') goto end
endif


mapline:
#######################################################################
# 10. Map the sources (spectral line)
#######################################################################

if ("$velmapfq" == "") goto mapcont
if ($#sources == 0) then
    echo "No source specified, mapping $gaincals[1]"
    set sources = ($gaincals[1])
endif

# Make the velocity datacubes:
foreach source ($sources)
  foreach mapfq ($velmapfq)
    # ------- Set up datacube
    uvspec vis=$source.$mapfq.cal axis=chan,amp log=uvspec.tmp \
        interval=10000 options=nobase,avall
    if ($imline == "") then
        set chani=`head -1 uvspec.tmp | awk '{printf "%d", $1}'`
        set chanf=`tail -1 uvspec.tmp | awk '{printf "%d", $1}'`
        echo "The first and last channels are $chani and $chanf"
        set nchan=`calc -i "(($chanf-$chani+1)/$chbin)-0.5"`
        echo "Binning by $chbin leads to $nchan channels"
    else
        set nchan=`echo $imline | awk -F, '{print $2}'`
    endif
    
#     uvlist vis=$source.$mapfq.cal options=spectra log=uvlist.tmp
#     set vi=`grep 'Start velocity' uvlist.tmp | head -1 | awk '{print int($4)}'`
#     set vf=`grep 'End velocity' uvlist.tmp | head -1 | awk '{print int($4)}'`
#     set vinc=`grep 'city increment' uvlist.tmp | head -1 | awk '{print $3}'`
#     if ($delv == "") then
#         set dv = `calc -f f6.2 "abs($vinc*2)"`
#     else
#         set dv = $delv
#     endif
#     if ($vi > $vf) then
#         set tmp=$vf
#         set vf=$vi
#         set vi=`expr $tmp + 1`
#     else
#         set vi=`expr $vi + 1`
#     endif
#     if (($nchan == "") && ($vstart == "")) then
#         set nch=`calc -i "0.75*($vf-($vi))/$dv"`
#         set v0=`calc -i "$vi+($vf-($vi))/8"`
#     else if ($vstart == "") then
#         set nch=$nchan
#         set v0=`calc -i "($vi+($vf)-$nchan*$dv)/2"`
#     else if ($nchan == "") then
#         set nchan=`calc -i "($vf-($vstart)/$dv)"`
#         set v0=$vstart
#     else
#         set nch=$nchan
#         set v0=$vstart
#     endif
#         
    # ------- subtract continuum if desired
    if ($contsub == "y") then
        rm -rf $source.$mapfq.uvlin
        if ($uvlinch != "") then
            echo "Using line-free channels $uvlinch"
            uvlin vis=$source.$mapfq.cal out=$source.$mapfq.uvlin mode=line \
                order=$fitorder chans=$uvlinch
        else
            uvlin vis=$source.$mapfq.cal out=$source.$mapfq.uvlin mode=line \
                order=$fitorder
        endif
        uvspec vis=$source.$mapfq.uvlin device=$source.$mapfq.spec2.ps/cps \
            axis=vel,amp interval=10000 nxy=2,3 select="pol($impol)"
#         if ($showplots == "y") then
#             gv $source.$mapfq.spec2.ps &
#         endif
        set infile=uvlin
    else
        set infile=cal
    endif
    
    # ------- Make the datacube
    echo "Working on $source.$mapfq.$infile" | tee -a $log
    rm -rf $source.$mapfq.map $source.$mapfq.beam
    if ($no6ant == "y") then
        set imsel="pol($impol),-ant(6)"
    else
        set imsel="pol($impol)"
    endif
    set imagepar = ""
    if ($cellsz != "") set imagepar = "cell=$cellsz"
    if ($imagsz != "") set imagepar = "$imagepar imsize=$imagsz"
    if ($imline == "") then
        set imline = "chan,$nchan,$chani,$chbin,$chbin"
    endif
    set options="double,systemp"
    if (`echo $source | grep mos` != "") then
        set options="$options,mosaic"
    endif
    echo "INVERT: $imagepar line=$imline options=$options" | tee -a $log    
    invert vis=$source.$mapfq.$infile map=$source.$mapfq.map \
        beam=$source.$mapfq.beam $imagepar sup=0 slop=0.5 \
        line=$imline select=$imsel options=$options

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

    # Normalize the maps
    rm -rf $source.$mapfq.mapnrm
    maths exp="$source.$mapfq.map/$source.$mapfq.sennorm" out=$source.$mapfq.mapnrm
#     rm -rf $source.$mapfq.senmsk
#     maths exp=$source.$mapfq.sen out=$source.$mapfq.senmsk \
#         mask="<$source.$mapfq.sengn>.lt.$maxsen"
    set drms=`sigest in=$source.$mapfq.mapnrm  | grep Est | awk '{print $4}'`
    set nmin4 = `expr $nchan - 4`
    set exec="sigest in=$source.$mapfq.mapnrm"
    set drmsi=`$exec region="image(1,5)" | grep Est | awk '{print $4}'`
    set drmsf=`$exec region="image($nmin4,$nchan)" | grep Est | awk '{print $4}'`
    set drms = `calc -f f7.5 "($drmsi+$drmsf)/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.$mapfq.mapnrm/rms value=$drms type=real
    puthd in=$source.$mapfq.mapnrm/scale value=$scale type=real
    echo " " | tee -a $log

    set slev = `calc "$drms*2"`
    set maxval=`histo in=$source.$mapfq.mapnrm | grep Max | awk '{print $3}'`
    cgdisp in=$source.$mapfq.mapnrm,$source.$mapfq.mapnrm type=c,p \
        device=$source.$mapfq.mapnrm.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=3,3 >! cgdisp.tmp
    mkpdf $source.$mapfq.mapnrm
    set pdflist=($pdflist $source.$mapfq.mapnrm.pdf)
#     if ($showplots == "y") then
#         gv $source.$mapfq.mapnrm.ps &
#     endif
  end
end

if ($interact != 'n') then
    echo ""
    echo -n "CLEAN the datacubes [y]? "
    set resp = $<
    if ($resp == 'n' || $resp == 'N') goto endline
endif

# Clean the maps:
foreach source ($sources)
  foreach mapfq ($velmapfq)
    rm -rf $source.$mapfq.cln
    histo in=$source.$mapfq.map region='quarter(1,3)' >! histo.tmp
    set rms = `gethd in=$source.$mapfq.mapnrm/rms`
    echo "The rms is $rms"
    set cutoff=`calc -f f6.4 "$rmsfac*$rms"`
    echo "Using a CLEAN cutoff of $cutoff"  | tee -a $log
    if (`echo $source | grep mos` != "") then
        mossdi map=$source.$mapfq.map beam=$source.$mapfq.beam cutoff=$cutoff \
            out=$source.$mapfq.cln niters=10000
    else
        clean map=$source.$mapfq.map beam=$source.$mapfq.beam cutoff=$cutoff \
            out=$source.$mapfq.cln niters=1000 region=quarter
    endif
    rm -rf $source.$mapfq.cm
    restor map=$source.$mapfq.map beam=$source.$mapfq.beam \
        model=$source.$mapfq.cln out=$source.$mapfq.cm
    set bmaj=`gethd in=$source.$mapfq.cm/bmaj format=arcsec | awk '{printf "%.2f",$1}'`
    set bmin=`gethd in=$source.$mapfq.cm/bmin format=arcsec | awk '{printf "%.2f",$1}'`
    set bpa=`gethd in=$source.$mapfq.cm/bpa | awk '{printf "%.1f",$1}'`
    echo "Generated $source.$mapfq.cm" | tee -a $log
    echo "Used restoring beam of $bmaj x $bmin arcsec, pa $bpa" | tee -a $log
    copyhd in=$source.$mapfq.map out=$source.$mapfq.cm items=mask
    rm -rf $source.$mapfq.cmnorm
    maths exp="$source.$mapfq.cm/$source.$mapfq.sennorm" out=$source.$mapfq.cmnorm
    set nmin4 = `expr $nchan - 4`
    set exec="sigest in=$source.$mapfq.cmnorm"
    set crmsi=`$exec region="image(1,5)" | grep Est | awk '{print $4}'`
    set crmsf=`$exec region="image($nmin4,$nchan)" | 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=$source.$mapfq.cmnorm/rms value=$crms type=real

    # Produce gain-corrected image cube and noise cube
    set minval=`histo in=$source.$mapfq.sen region='image(1)'|grep Min|awk '{print $3}'`
    rm -rf $source.$mapfq.nse
    maths exp="<$source.$mapfq.sengn>*$crms" out=$source.$mapfq.nse \
        mask="<$source.$mapfq.sengn>.lt.$msklev"
#     rm -rf $source.$mapfq.nse $source.$mapfq.nsesub $source.$mapfq.nse1
#     imsub in=$source.$mapfq.nse region='image(1)' out=$source.$mapfq.nsesub
#     imblr in=$source.$mapfq.nsesub out=$source.$mapfq.nse1 value=1
    rm -rf $source.$mapfq.cmmsk
    maths exp="<$source.$mapfq.cm>/<$source.$mapfq.gain>" \
        region="mask($source.$mapfq.nse)" out=$source.$mapfq.cmmsk
    rm -rf $source.$mapfq.cmnse
    maths exp="<$source.$mapfq.nse>" region="mask($source.$mapfq.nse)" \
        out=$source.$mapfq.cmnse
    foreach type (cmmsk cmnse)
        rm -f $source.$mapfq.$type.fits
        fits in=$source.$mapfq.$type out=$source.$mapfq.$type.fits op=xyout
    end

    # Produced noise-flattened residual cube
    rm -rf $source.$mapfq.res
    restor map=$source.$mapfq.map beam=$source.$mapfq.beam model=$source.$mapfq.cln \
        out=$source.$mapfq.res mode=residual
    rm -rf $source.$mapfq.resnorm
    maths exp="$source.$mapfq.res/$source.$mapfq.sennorm" out=$source.$mapfq.resnorm
    copyhd in=$source.$mapfq.cm out=$source.$mapfq.resnorm items=bmaj,bmin,bpa
    set rflux=`histo in=$source.$mapfq.resnorm | grep Flux | awk '{print $6}'`
    set rrms=`sigest in=$source.$mapfq.resnorm | grep Est | awk -v c=$crms '{print $4/c}'`
    set dv=`gethd in=$source.$mapfq.map/cdelt3`
    set rintflux=`calc -f f8.2 "$rflux*abs($dv)"`
    set cflux=`histo in=$source.$mapfq.cln | grep Flux | awk '{print $6}'`
    set cintflux=`calc -f f8.2 "$cflux*abs($dv)"`
    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
    
    # Plot the channel maps
    set maxval=`histo in=$source.$mapfq.cmnorm | grep Max | awk '{print $3}'`
    cgdisp in=$source.$mapfq.cmnorm,$source.$mapfq.cmnorm type=c,p \
        device=$source.$mapfq.cmnorm.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 $source.$mapfq.cmnorm
    set pdflist=($pdflist $source.$mapfq.cmnorm.pdf)
#     if ($showplots == "y") then
#         gv $source.$mapfq.cmnorm.ps &
#     endif
  end
end

endline:
if ("$conmapfq" == "" && "$testcals" == "") goto endcont
if ($interact != 'n') then
    echo ""
    echo -n "Continue with continuum maps [y]? "
    set resp = $<
    if ($resp == 'n' || $resp == 'N') goto endcont
endif


mapcont:
#######################################################################
# 11. Map the sources (continuum)
#######################################################################

if ($#sources == 0) then
    echo "No source specified, mapping $gaincals[1]"
    set sources = ($gaincals[1])
endif

# Make the continuum maps:
foreach source ($sources $testcals)
  if (`echo $testcals | grep $source` != "") then
    set dofreq = ($freq1 $freq2)
  else
    set dofreq = ($conmapfq)
  endif
  foreach mapfq ($dofreq)
    # ------- Set up datacube
    uvlist vis=$source.$mapfq options=spectra log=uvlist.tmp
    set totchan=`grep Number uvlist.tmp | awk '{print $5}'`
    # ------- Image all channels
    rm -rf $source.$mapfq.map $source.$mapfq.beam
    if ($no6ant == "y") then
        set imsel="pol($impol),-ant(6)"
    else
        set imsel="pol($impol)"
    endif
    set options="mfs,double,systemp"
    if (`echo $source | grep mos` != "") then
        set options="$options,mosaic"
    endif
    invert vis=$source.$mapfq.cal map=$source.$mapfq.map \
        beam=$source.$mapfq.beam cell=$cellsz imsize=$imagsz sup=0 slop=0.5 \
        line=chan,1,1,$totchan select=$imsel options=$options
    cgdisp in=$source.$mapfq.map type=c device=$source.$mapfq.map.ps/cps \
        labtyp=hms,dms nxy=1,1 slev=p,10 levs1=1,2,3,4,5,6,7,8,9 \
        csize=0.8,0.8 options=full,mirr cols1=1
#     cgdisp in=$source.$mapfq.beam type=c device=$source.$mapfq.beam.ps/ps \
#         labtyp=arcsec nxy=1,1 slev=p,10 levs1=1,2,3,4,5,6,7,8,9 \
#         csize=0.8,0.8 options=full,mirr
    mkpdf $source.$mapfq.map
    set pdflist=($pdflist $source.$mapfq.map.pdf)
#     if ($showplots == "y") then
#         gv $source.$mapfq.map.ps &
#     endif
  end
end

if ($interact != 'n') then
    echo ""
    echo -n "CLEAN the continuum maps [y]? "
    set resp = $<
    if ($resp == 'n' || $resp == 'N') goto end
endif

# Clean the maps:
foreach source ($sources $testcals)
  if (`echo $testcals | grep $source` != "") then
    set dofreq = ($freq1 $freq2)
  else
    set dofreq = ($conmapfq)
  endif
  foreach mapfq ($dofreq)
    rm -rf $source.$mapfq.cln
    set rms = `gethd in=$source.$mapfq.map/rms`
    if ($rms == "") then
        set rms=`histo in=$source.$mapfq.map | grep Rms | awk '{print $4}'`
    endif
    set cutoff=`calc -f f6.4 "$rmsfac*$rms"`
    echo "Using a CLEAN cutoff of $cutoff"  | tee -a $log
    if (`echo $source | grep mos` != "") then
        mossdi map=$source.$mapfq.map beam=$source.$mapfq.beam cutoff=$cutoff \
            out=$source.$mapfq.cln niters=10000
    else
        clean map=$source.$mapfq.map beam=$source.$mapfq.beam cutoff=$cutoff \
            out=$source.$mapfq.cln niters=1000 
    endif
    rm -rf $source.$mapfq.cm
    restor map=$source.$mapfq.map beam=$source.$mapfq.beam \
        model=$source.$mapfq.cln out=$source.$mapfq.cm
    cgdisp in=$source.$mapfq.cm type=c device=$source.$mapfq.cm.ps/cps \
        labtyp=hms,dms nxy=1,1 slev=p,10 levs1=1,2,3,4,5,6,7,8,9 \
        csize=0.8,0.8 options=full,mirr,beambl cols1=1
    mkpdf $source.$mapfq.cm
    set pdflist=($pdflist $source.$mapfq.cm.pdf)
#     if ($showplots == "y") then
#         gv $source.$mapfq.cm.ps &
#     endif
  end
end

endcont:
if ($interact != 'n') then 
    echo ""
    echo -n "Purge unwanted files? [n]? "
    set resp = $<
    if ($resp == 'y' || $resp == 'Y') goto purge
else
    if ($purge == 'y') goto purge
endif
goto end


purge:
#######################################################################
# 12. Purge unwanted files
#######################################################################

# Move plots to subdirectory
if !(-e PLOTS) mkdir PLOTS
mv *.ps PLOTS
gzip -f PLOTS/*.ps &
if (-e RELCAL) gzip -f RELCAL/*.ps &

foreach obj ($sources $gaincals $testcals $fluxcal $bpcal)
    rm -rf $obj
end

foreach source ($sources $testcals)
    foreach freq ($freq1 $freq2)
        rm -rf $source.$freq.beam $source.$freq.cln
    end
end

foreach freq ($freq1 $freq2)
    rm -rf $bpcal.$freq.cal 
end

rm -rf atfix.[1-9]
rm -f *.med *.out uvlist.log

end:
echo "Output to $log"
rm -f tmp.* *.tmp
gs -q -dSAFER -dNOPAUSE -dBATCH -sOutputFile=allplots.$today.pdf -sDEVICE=pdfwrite \
    -c .setpdfwrite -f $pdflist
echo Plotted images output to allplots.$today.pdf
rm -f $pdflist