#!/bin/csh -f
#
# edgecalib.csh - CARMA EDGE calibration script (Tony Wong, UIUC)
#
# All parameters can be set on the command line:
#   edgecalib.csh file= [refant=] [restfreq=]
#
# History
#  2014oct24: original
#  2014oct25: use mfboot mode=vector for MWC349.
#  2014oct27: fix bug in determining velocity limits
#  2014oct31: improved plotting and cleaning; support for 13co; bug fixes
#  2014nov01: plotting tweaks
#  2014nov02: peak brightness map in mK
#  2014nov06: output fluxes for all calibrators
#  2014nov13: set a default flux for 3C273
#  2014dec03: unused fluxcals are considered gaincals
#  2014dec04: allow fluxflux to override mwcflux and flux3c; html out
#  2014dec10: prevent varplt from trying to send plots to screen
#  2015jan18: bug fixes for apfile mode
#  2015feb23: don't add fluxcal to gaincal list if it's already there;
#             use uvredo to recompute Doppler corrections; split CO and 13CO
#  2015mar02: revert to triple mode for MWC349
#  2015jul14: allow more flexible correlator settings
#  2015jul17: output fluxflux to reduclog
#--------------------------------------------------------------------------

set version=2015jul17
set log=edgecalib.log
set showplots=y

### INFORMATION ABOUT THIS TRACK
set file=           # REQUIRED, several may be specified
set sources=        # more than one allowed [choose automatically]
set gaincals=       # more than one allowed [choose automatically]
set passcals=       # more than one allowed [choose automatically]
set fluxcal=        # only one allowed [choose automatically]
set winselect=9,10,11,12,13,14,15,16      # Window selection; others discarded
set wide=3          # (re-numbered) spectral window for phase calibration
set f12=115.271204
set f13=110.201353
set restfreq=$f13,$f13,$f13,$f12,$f12,$f12,$f12,$f12 # (re-numbered) rest frequencies
set win12=4,5,6,7,8 # 12CO windows
set win13=1,2,3     # 13CO windows
set apfile=         # antenna position file, if correction required
set fluxflux=       # Flux of primary, if not a planet
set mwcflux=1.2     # Adopted flux of MWC349
set flux3c=12       # Adopted flux of 3C273
set arr=e           # Array label for HTML link

### USER PREFERENCES
set dolinecal=y     # Apply linelength calibration?
set hanning=n       # Hanning smooth the data?
set tsyslim=1000    # Flag data with Tsys larger than this
set bpedge=0        # Edge channels to drop during mfcal
set polyfit=0       # Order of polynomial fit to bandpass if >0
set edge=2          # Edge channels to flag on source
set refant=8        # Ref ant for selfcal (first available if unset)
set pcalint=6       # selfcal interval (mins) for phasecal
set delv=40         # channel spacing in km/s
set imsize=129      # image dimension for INVERT
set cell=2          # pixel size for INVERT
set rmsfac=3        # clean cutoff in units of map RMS
set senlim=2.5      # sensitivity contour level for senmsk
set xyreg='arc,box(-60,-60,60,60)'  # region for plotting (entire field if unset)
set exinvert=       # extra parameters for INVERT
set do_linecal=y    # Apply the linelength calibration?
set smooth_int=     # Time interval for smoothing linelengths [unset]
set zero_ants=      # Antennas to zero linelengths [unset]
set do_csflag=n     # flag shadowed data? [n]
set do_tsysflag=y   # flag high Tsys values? [y]
set do_deconv=y     # deconvolve maps? [y]
set purge=n         # Delete scratch files and compress plots? [n]

### CUSTOM FLAGGING
# Example: 
#     flagsrc=all,3C84
#     flagsel='ant(1)'+'time(10:00,11:00)'
# will flag ant(1) for all data and times 10:00-11:00 for 3C84.
set flagsrc=        # optional source(s) for extra flagging
set flagsel=        # optional selection(s) for extra flagging
set flaglin=all     # optional linetype specification, default is all channels

### Optional flux calibration parameters (planet bootstrapping)
# interval to average ampgains over before applying to planet
#set fluxavint=
# optional time range where cal has similar elev/Tsys to planet
#set gainsel=

### LOCAL SETTINGS
alias date 'date "+%m/%d/%y %T %Z"'
alias xspause 'set tmpr=$<'

### 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 INPUTS

if ($file == "") then
    echo "### Fatal Error: file= must be specified"
    exit 1
endif

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

#--------------------------------------------------------------------------
# 1. Diagnostics:
#--------------------------------------------------------------------------

# Run an initial listobs including NOISE integrations:
listobs vis=$file log=listobs1

# Fix the baselines if requested
if ($apfile != "") then
	set nfiles=`echo $file | awk -F, '{print NF}'`
    if ($nfiles > 1) then
        set outfile = `echo $file | sed 's/$/_c/' | sed 's/,/_c /g'`
    else
        set outfile = `echo $file | sed 's/$/_c/'`
    endif
    rm -rf $outfile
    uvedit vis=$file apfile=$apfile
    echo "Applied antenna position file $apfile"
else
    set outfile=$file
endif

# Recompute Doppler terms.  Window selection if requested
if ($winselect != '') then
    rm -rf uvdata.wins
    uvcat vis="$outfile" out=uvdata.wins select="win($winselect)" \
        options=nowide
    rm -rf uvredo.mir
    uvredo vis=uvdata.wins out=uvredo.mir select='-source(NOISE)' \
        options=velocity
    set wins=`echo "$winselect" | sed 's/,/ /g'`
    set nwide=$#wins
    rm -rf uvdata.mir
    uvwide vis=uvredo.mir out=uvdata.mir nwide=$nwide
    echo "Selected spectral windows $winselect"
else
    rm -rf uvdata.mir uvredo.mir
    uvcat vis="$outfile" out=uvdata.mir
    uvredo vis=uvdata.mir out=uvredo.mir select='-source(NOISE)' \
        options=velocity
    rm -rf uvdata.mir
    uvwide vis=uvredo.mir out=uvdata.mir nwide=16
endif

# Fix rest frequencies
if ($restfreq != "") then
    rm -rf uvputhd.mir
    uvputhd vis=uvdata.mir out=uvputhd.mir hdvar=restfreq varval=$restfreq
    rm -rf uvdata.mir
    mv uvputhd.mir uvdata.mir
    echo "Entered rest freqs of $restfreq"
endif

# Hanning smoothing if requested
if ($hanning == 'y') then
    uvcal vis=uvdata.mir options=hanning out=uvdata.hann
    rm -rf uvdata.mir
    mv uvdata.hann uvdata.mir
    echo "Applied hanning smoothing to all data"
endif

# Generate the chronology of observations
listobs vis=uvdata.mir log=listobs2
cat listobs2
if !($?obsdate) set obsdate=`awk '/Chronology/ {print $5}' listobs2`

# Choose a new reference antenna if default is unavailable
set reffound = 0
set beglin=`awk '/Sys Temps/ {print NR}' listobs2`
@ beglin ++
set ants=`sed -n "${beglin}p" listobs2 | sed 's/\(.*[deg,mode]\)//'`
foreach ant ($ants)
    if ($refant == $ant) then
        set reffound = 1
    endif
end
if ($reffound == 0) then
    echo "*** WARNING: Requested reference antenna $refant not found"
    echo "*** WARNING: Using antenna $ants[1] instead"
    set refant = $ants[1]
endif

# Get the source list and their intents
set beglin2=`awk '/Observed Sources/ {print NR+2}' listobs2`
set endlin2=`awk '/Frequency Set-up/ {print NR-2}' listobs2`
set allobs=`sed -n "$beglin2,${endlin2}p" listobs2 | awk '{print $1}'`
set purpose=`sed -n "$beglin2,${endlin2}p" listobs2 | awk '{print $2}'`
set i=1
set defsources=
set defgaincals=
set defpasscals=
set deffluxcal=
echo ""
while ($i <= $#allobs)
    if ($purpose[$i] =~ *B*) set defpasscals = ($defpasscals $allobs[$i])
    if ($purpose[$i] =~ *G*) set defgaincals = ($defgaincals $allobs[$i])
    if ($purpose[$i] =~ *S*) set defsources = ($defsources $allobs[$i])
    if ($purpose[$i] =~ *F*) set deffluxcal = ($deffluxcal $allobs[$i])
    if ($purpose[$i] =~ *A*) set defsources = ($defsources $allobs[$i])
    if ($purpose[$i] =~ *O*) set defgaincals = ($defgaincals $allobs[$i])
    echo "Source $allobs[$i] has purpose $purpose[$i]"
    set i = `expr $i + 1`
end

# If there are 2 or more passcals, remove any planets:
if ($#defpasscals > 1) then
    foreach passcal ($defpasscals)
        if ($passcal =~ {MARS,JUPITER,SATURN,URANUS,NEPTUNE,MWC349}) then
            if ($#defpasscals > 1) then
                set defpasscals=`echo $defpasscals | sed "s/$passcal//"`
            endif
        endif
    end
endif

# If there are 2 or more gaincals, remove any planets:
if ($#defgaincals > 1) then
    foreach gaincal ($defgaincals)
        if ($gaincal =~ {MARS,JUPITER,SATURN,URANUS,NEPTUNE,MWC349}) then
            if ($#defgaincals > 1) then
                set defgaincals=`echo $defgaincals | sed "s/$gaincal//"`
            endif
        endif
    end
endif

# If there are 2 or more fluxcals, choose a planet:
if ($#deffluxcal > 1) then
    foreach obj ($deffluxcal)
        if ($obj !~ {MARS,JUPITER,SATURN,URANUS,NEPTUNE,MWC349}) then
            if ($#deffluxcal > 1) then
                set deffluxcal=`echo $deffluxcal | sed "s/$obj//"`
                set gaincalstmp=`echo $defgaincals | sed 's/ /,/g'`
                if ($obj !~ {$gaincalstmp}) then
                    set defgaincals=($defgaincals $obj)
                endif
            endif
        endif
    end
    if ($#deffluxcal > 1) set deffluxcal=$deffluxcal[1]
endif

# State the default classifications
echo ""
echo "Default assignments, based on source intent:"
echo "sources : $defsources"
echo "gaincals: $defgaincals"
echo "passcals: $defpasscals"
echo "fluxcal: $deffluxcal"

# Fill in the source names if unspecified
if ($sources == '') set sources=`echo $defsources | sed 's/ /,/g'`
if ($gaincals == '') set gaincals=`echo $defgaincals | sed 's/ /,/g'`
if ($passcals == '') set passcals=`echo $defpasscals | sed 's/ /,/g'`
# If there is no passcal, choose the first gaincal:
if ($passcals == '') set passcals=$defgaincals[1]
if ($fluxcal == '') set fluxcal=$deffluxcal

# Compile lists of sources
set sourcelis=`echo "$sources" | sed 's/,/ /g'`
set gaincalis=`echo "$gaincals" | sed 's/,/ /g'`
set passcalis=`echo "$passcals" | sed 's/,/ /g'`
set allcals = ()
set allobjs = ()
foreach obj ($gaincalis $passcalis $fluxcal)
    set val=`echo $allcals | grep $obj`
    if ("$val" == "") set allcals = ($obj $allcals)
end
foreach obj ($sourcelis $gaincalis $passcalis $fluxcal)
    set val=`echo $allobjs | grep $obj`
    if ("$val" == "") set allobjs = ($obj $allobjs)
end

# Open log and output some interesting information
if !(-e $log) then
    mv $log $log.bak
endif
echo ""
echo "edgecalib.csh $version executed by $user on `date`" | tee $log
echo "" | tee -a $log
echo "Project: $file  Obs Date: $obsdate" | tee -a $log
echo "" | tee -a $log
echo "Sources: $sources  Gain cals: $gaincals" | tee -a $log
echo "Flux cal: $fluxcal   Passband cal: $passcals" | tee -a $log
if ($fluxflux != "") echo "Flux of $fluxcal set to $fluxflux" | tee -a $log
echo "" | tee -a $log
echo "===============================================================" \
    >> $log

# Get the elapsed project time, note this includes NOISE integrations
echo "" | tee -a $log
echo "------ diagnostics run on `date`" >> $log
set startlin=`awk '/hhmmss/ {print NR+1}' listobs1`
set endlin=`wc listobs1 | awk '{print $1-1}'`
set starttime=`sed -n "${startlin}p" listobs1 | awk '{print $2}'`
set stoptime=`sed -n "${endlin}p" listobs1 | awk '{print $2}'`
set starthours=`echo $starttime | awk '{split ($1,h,""); print h[1] h[2]}'`
set startminutes=`echo $starttime | awk '{split ($1,h,""); print h[3] h[4]}'`
set stophours=`echo $stoptime | awk '{split ($1,h,""); print h[1] h[2]}'`
set stopminutes=`echo $stoptime | awk '{split ($1,h,""); print h[3] h[4]}'`
set addminutes=`sed -n "${endlin}p" listobs1 | awk '{print $4}'`
set starttime = `calc -f f6.3 "$starthours+$startminutes/60.0"`
set stoptime  = `calc -f f6.3 "$stophours+($stopminutes+$addminutes)/60.0"`
set stoptime  = `echo $starttime $stoptime | awk '{if ($1 > $2) print ($2 + 24); else print $2}'`
set totaltime = `calc -f f6.1 "$stoptime-$starttime"`
echo "Total project time $totaltime hrs from start to finish" | tee -a $log

# Report frequency setup
if ($restfreq != "") then
    echo "Entered rest frequency of $restfreq" >> $log
endif
if ($winselect != '') then
    echo "Selected spectral windows $winselect" >> $log
endif
if ($hanning == 'y') then
    echo "Applied hanning smoothing to all data" >> $log
endif
uvlist vis=uvdata.mir select="-auto" options=spectra log=spectra.log
cat spectra.log
prthd in=uvdata.mir options=full log=tmp.prthd 
set beglin=`awk '/Bandwidth/ {print NR}' tmp.prthd`
@ beglin ++
set endlin=`awk '/J2000/ {print NR}' tmp.prthd`
set endlin=`expr $endlin - 3`
sed -n "$beglin,${endlin}p" tmp.prthd > tmp.corr2
set bandwsrc = `cat tmp.corr2 | awk '{printf "%5d", 1000*sqrt($3*$3)}'`
echo "Correlator setup:" | tee -a $log
cat tmp.corr2 | tee -a $log
echo "Using window $wide for calibration" | tee -a $log

# Determine mean rmspath, tau230
varplt vis=uvdata.mir yaxis=tau230 device=tau230.ps/ps log=tau.log
if (-e tau.log) then
  set tauavg=`cat tau.log | tail +5 | awk '{s += $3} END {printf "%.2f", s/NR}'`
  echo "Mean value of tau230 is $tauavg" | tee -a $log
  rm -f tau.log
endif
varplt vis=uvdata.mir yaxis=rmspath device=rmspath.ps/ps log=rmp.log
if (-e rmp.log) then
  set rmpavg=`cat rmp.log | tail +5 | awk '{s += $3} END {printf "%.0f", s/NR}'`
  echo "Mean value of rmspath is $rmpavg" | tee -a $log
  rm -f rmp.log
endif
smavarplt vis=uvdata.mir yaxis=systemp nxy=5,3 device=tsys1.ps/cps 
if ($showplots == 'y') gv tsys1.ps &
smavarplt vis=uvdata.mir yaxis=systemp nxy=5,3 device=tsys2.ps/cps \
    yrange=0,1000
if ($showplots == 'y') gv tsys2.ps &

# --- ADDITIONAL DIAGNOSTICS IN INTERACTIVE MODE

if ($showplots == 'y') then
    # LOOK AT CALIBRATORS, UV COVERAGE
    echo ""
    echo "UV-coverage on sources:"
    uvplt vis=uvdata.mir select="source($sources),-auto" axis=uc,vc \
        device=/xs options=nobase,equal nxy=1
    echo "Calibrator amps and phases:"
    smauvplt vis=uvdata.mir select="source($gaincals),-auto" device=/xs \
        line=wide,1,$wide axis=time,amp nxy=5,3 title=BL
    xspause
    smauvplt vis=uvdata.mir select="source($gaincals),-auto" device=/xs \
        line=wide,1,$wide axis=time,phase nxy=5,3 yrange=-180,180 title=BL
    xspause
    # LOOK AT SOURCE SPECTRA
    echo ""
    echo "Uncalibrated source spectra:"
    uvspec vis=uvdata.mir select="source($sources),-auto" interval=1000 \
        nxy=3,3 device=/xs 
    xspause
endif

#--------------------------------------------------------------------------
# 2. Linelength calibration:
#--------------------------------------------------------------------------

dolinecal:

# Log entry:
echo "" | tee -a $log
echo "------ linecal run on `date`" >> $log

# Get and examine linelengths:
linecal vis=uvdata.mir
gpplt vis=uvdata.mir yaxis=phase yrange=-180,180 options=wrap,dots \
    device=linecal.ps/ps nxy=3,5
if ($showplots == 'y') gv linecal.ps &

# Time smoothing of linelengths (optional):
if ($smooth_int != "") then
    gpaver vis=uvdata.mir interval=$smooth_int
    if ($showplots == 'y') then
        gpplt vis=uvdata.mir device=/xs yaxis=ph nxy=3,5 yrange=-180,180 \
            options=wrap,dots
    endif
    echo "Smoothed linelengths over $smooth_int minute interval" | tee -a $log
    echo ""
endif

# Nullify bad linelengths (optional):
if ($zero_ants != "") then
    gpedit vis=uvdata.mir select="ant($zero_ants)" options=replace
    echo "Zeroed phases for antennas $zero_ants" | tee -a $log
    echo ""
    if ($showplots == 'y') then
        gpplt vis=uvdata.mir device=/xs yaxis=ph nxy=3,5 yrange=-180,180 \
            options=wrap,dots
    endif
endif

# Don't apply linelengths at all (optional):
if ($do_linecal == 'n') then
    echo "NO LINELENGTH CALIBRATION APPLIED." | tee -a $log
    delhd in=uvdata.mir/gains
endif

#--------------------------------------------------------------------------
# 3. Data flagging:
#--------------------------------------------------------------------------

doflag:

# Log entry:
echo "" | tee -a $log
echo "------ uvflag run on `date`" >> $log

# Flagging of shadowed data
if ($do_csflag != "n") then 
  echo ""
  echo "Flagging shadowed data..."
  csflag vis=uvdata.mir carma=true
  echo "Shadowed data flagged with csflag." >> $log
endif

# Flagging of high Tsys data
if ($do_tsysflag != "n") then 
  echo ""
  echo "Flagging high Tsys data..."
  uvflag vis=uvdata.mir tsys=$tsyslim flagval=f
  echo "High Tsys data (>$tsyslim) flagged with uvflag." >> $log
endif

# Custom flagging specified in script
set flagsrc=`echo "$flagsrc" | sed 's/,/ /g'`
set flagsel=`echo "$flagsel" | sed 's/+/ /g'`
set flaglin=`echo "$flaglin" | sed 's/+/ /g'`
set i = $#flagsrc
if ($i == 0) goto edgeflag
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
if ($#flagsrc > $#flaglin) then
    set j=$#flaglin
    while ($j < $i)
	    set flaglin = ($flaglin $flaglin[$j])
	    set j=`expr $j + 1`
    end
endif
echo "Flag sources: $flagsrc"
echo "Flag selection: $flagsel"
echo "Line selection: $flaglin"
while ($i > 0)
    if ($flaglin[$i] == 'all') then
        set extra=
	else
	    set extra="line=$flaglin[$i]"
	endif
    if ($flagsrc[$i] == 'all') then
	    echo ""
	    echo "Flagging data for all srcs, select=$flagsel[$i] $extra" \
	        | tee -a $log
        uvflag vis=uvdata.mir select="$flagsel[$i]" flagval=f $extra
    else if ($flagsrc[$i] != "") then
	    echo ""
	    echo "Flagging data for source($flagsrc[$i]),$flagsel[$i] $extra" \
	        | tee -a $log
	    uvflag vis=uvdata.mir select="source($flagsrc[$i]),$flagsel[$i]" \
	        flagval=f $extra
    endif
    set i = `expr $i - 1`
end

edgeflag:
# Flag edge channels on sources
if ($edge != "") then
    uvflag vis=uvdata.mir select="source($sources)" flagval=f edge=$edge
    echo "Flagged $edge edge channels on $sources" >> $log
endif
echo ""

# Remake the visibility sets if desired:
# if ($flagredo == 'y') then
#     foreach flgsrc ($source $fluxcal phasecals)
#         uvcat vis=$flgsrc out=$flgsrc.tmp options=unflagged
#         rm -rf $flgsrc
#         mv $flgsrc.tmp $flgsrc
#     end
# endif

# Apply linelength gains to unflagged data
rm -rf lineflg.mir
uvcat vis=uvdata.mir out=lineflg.mir select="-auto" options=unflagged

#--------------------------------------------------------------------------
# 4. Calibrate the passband:
#--------------------------------------------------------------------------

dopass:

# Inspect phase and amplitude across passband:
if ($showplots == 'y') then
    echo "Plotting the amps and phases of the pb calibrator vs. chan..."
    uvspec vis=lineflg.mir interval=600 device=/xs nxy=5,3 axis=chan,amp \
        select="source($passcals)"
    xspause
    uvspec vis=lineflg.mir interval=600 device=/xs nxy=5,3 axis=chan,pha \
        select="source($passcals)" yrange=-180,180
    xspause
endif

# Log entry:
echo "" | tee -a $log
echo "------ passband cal using $passcals on `date`" >> $log
echo "Used antenna $refant as reference for selfcal." >> $log
echo "Edge channels to ignore in each window: $bpedge" >> $log

# Calibrate the passband:
if ($polyfit != '0' && $polyfit != '') then
    smamfcal vis=lineflg.mir refant=$refant interval=1 edge=$bpedge \
        options=opolyfit polyfit=$polyfit select="source($passcals)"
else
    mfcal vis=lineflg.mir refant=$refant interval=1 edge=$bpedge \
        select="source($passcals)"
endif
smagpplt vis=lineflg.mir options=bandpass,nofit device=bpamp.ps/cps xaxis=chan \
    yrange=0,2 nxy=3,5
if ($showplots == 'y') gv bpamp.ps &
smagpplt vis=lineflg.mir options=bandpass,nofit,wrap device=bpph.ps/cps xaxis=chan \
    yaxis=phase yrange=-180,180 nxy=3,5
if ($showplots == 'y') gv bpph.ps &

# Apply the passband gains:
rm -rf pbcald.mir
uvcat vis=lineflg.mir out=pbcald.mir options=nocal
rm -rf pbcald.nw
uvwide vis=pbcald.mir out=pbcald.nw

#--------------------------------------------------------------------------
# 5. Phase calibration:
#--------------------------------------------------------------------------

docal:

# Log Entry:
echo "" | tee -a $log
echo "------ gain calibration run on `date`" >> $log
echo "Reference antenna: $refant" >> $log
echo "Solution interval: $pcalint min" >> $log

# Solve For Gains:
set gaincalis=`echo "$gaincals" | sed 's/,/ /g'`
foreach gaincal ($gaincalis)
    rm -rf $gaincal.sgains
    mselfcal vis=pbcald.nw refant=$refant line=wide,1,$wide interval=$pcalint \
        options=amp select="source($gaincal)" out=$gaincal.sgains
end
foreach gaincal ($gaincalis)
    if ($gaincal == $gaincalis[1]) then
        gpcopy vis=$gaincalis[1].sgains out=pbcald.nw mode=copy
    else
        gpcopy vis=$gaincal.sgains out=pbcald.nw mode=merge
    endif
end
puthd in=pbcald.nw/senmodel value='GSV' type=ascii
puthd in=pbcald.nw/interval value=0.1
gpplt vis=pbcald.nw device=gainph.ps/ps yrange=-180,180 nxy=5,3 \
    yaxis=phase options=wrap
if ($showplots == 'y') gv gainph.ps &
gpplt vis=pbcald.nw device=gainamp.ps/ps yrange=0,4 nxy=5,3
if ($showplots == 'y') gv gainamp.ps &

# Report gains in log:
echo "Mean amplitude gains, window ${wide}:" | tee -a $log
gplist vis=pbcald.nw | grep Means | sed 's/.*:[ ]*//g' | tee -a $log
echo "... these multiply the online values:" >> $log
uvlist vis=pbcald.nw options=var,full | grep jyperka >> $log

# Plot phases on calibrator:
uvplt vis=pbcald.nw axis=time,amp nxy=5,3 line=wide,1,$wide \
    device=calamps.ps/ps options=2p,source \
    size=1.8,5 select="source($gaincals)"
uvplt vis=pbcald.nw axis=time,phase nxy=5,3 line=wide,1,$wide \
    device=calphases.ps/ps options=2p,source yrange=-180,180 \
    size=1.8,5 select="source($gaincals)"

## Optional: map calibrator and check its size:
if ($showplots == 'y') then
    gv calamps.ps &
    gv calphases.ps &
    uvamp vis=pbcald.nw device=/xs line=wide,1,$wide bin=45,1,klam
#     echo ""
#     echo -n "Map the calibrator [n]? "
#     set resp = $<
#     if ($resp == "y") then
#         echo "Inverting the $gaincalis[1] data..." | tee -a $log
#         rm -rf $gaincalis[1].map $gaincalis[1].beam
#         invert vis=pbcald.nw map=$gaincalis[1].map beam=$gaincalis[1].beam \
#             imsize=256 cell=1 sup=0 line=wide,1,$wide options=systemp \
#             select="source($gaincalis[1])"
#         imfit in=$gaincalis[1].map region='arcsec,box(-15,-15,15,15)' \
#             object=gaussian > imfit1.tmp
#         set bmaj=`grep Major imfit1.tmp | awk '{print $4}'`
#         set bmin=`grep Minor imfit1.tmp | awk '{print $4}'`
#         echo "Gaussian fit to map : $bmaj arcsec x $bmin arcsec" | tee -a $log
#         imfit in=$gaincalis[1].beam region='arcsec,box(-15,-15,15,15)' \
#             object=gaussian > imfit2.tmp
#         set bmaj=`grep Major imfit2.tmp | awk '{print $4}'`
#         set bmin=`grep Minor imfit2.tmp | awk '{print $4}'`
#         echo "Gaussian fit to beam: $bmaj arcsec x $bmin arcsec" | tee -a $log
#         rm -rf $gaincalis[1].clean
#         clean map=$gaincalis[1].map beam=$gaincalis[1].beam out=$gaincalis[1].clean 
#         rm -rf $gaincalis[1].cm
#         restor map=$gaincalis[1].map beam=$gaincalis[1].beam \
#             model=$gaincalis[1].clean out=$gaincalis[1].cm
#         cgdisp in=$gaincalis[1].cm,$gaincalis[1].cm type=p,c nxy=1 device=/xs  \
#             labtyp=arcsec region='arcsec,box(-60,-60,60,60)' range=0,0,lin,3 \
#             slev=p,10 levs1=1,2,3,4,5,6,7,8,9 options=full,beambl,wedge,mirr
#         xspause
#     endif
endif

#--------------------------------------------------------------------------
# 6. Absolute flux calibration
#--------------------------------------------------------------------------

doflux:

if ($fluxcal == "") then
    echo "No flux calibrator specified."
    set bootfac=
    goto domaps
endif

# --- Gain fiddling: select a time range where calibrator is 
# --- close to planet's elevation and average gains in time
# if ($gainsel != "") then
#     echo "Selecting only gains in time range $gainsel" | tee -a $log
#     gpedit vis=$fcalfile select="-time($gainsel)" options=flag
# endif
# if ($fluxavint != "") then
#     echo "Averaging gains over $fluxavint min" | tee -a $log
#     gpaver vis=$fcalfile interval=$fluxavint options=scalar
# endif

# Bootstrap The Fluxes:
echo "" | tee -a $log
echo "------ flux calibration run on `date`" >> $log
echo "The flux calibrator is $fluxcal" | tee -a $log
if ($fluxcal =~ {MARS,JUPITER,SATURN,URANUS,NEPTUNE}) then
    mfboot vis=pbcald.nw select="source($fluxcal)" device=boot.ps/cps \
        line=wide,1,$wide clip=0.2 mode=scalar | tee tmp.mfboot
    if ($showplots == 'y') gv boot.ps &
else if ($fluxcal =~ {MWC349}) then
    if ($fluxflux == "") set fluxflux=$mwcflux
    echo "Adopting a flux of $fluxflux" | tee -a $log
    mfboot vis=pbcald.nw select="source($fluxcal)" \
        line=wide,1,$wide flux=$fluxflux mode=triple | tee tmp.mfboot
else if ($fluxcal =~ {3C273}) then
    if ($fluxflux == "") set fluxflux=$flux3c
    echo "Adopting a flux of $fluxflux" | tee -a $log
    mfboot vis=pbcald.nw select="source($fluxcal)" \
        line=wide,1,$wide flux=$fluxflux mode=triple | tee tmp.mfboot
else
    calflux source=$fluxcal freq=110 delfreq=30 date=0 | tee calflux.log
    if ($fluxflux == "") then
        set fluxflux=`grep Jy calflux.log | awk -F: '{print $3}' | awk '{print $1}'`
        echo "Adopting the most recent flux of $fluxflux" | tee -a $log
    endif
    mfboot vis=pbcald.nw select="source($fluxcal)" \
        line=wide,1,$wide flux=$fluxflux mode=triple | tee tmp.mfboot
endif
grep Scaling tmp.mfboot >> $log
set bootfac=`grep Scaling tmp.mfboot | awk '{print $5}'`
uvflux vis=pbcald.nw select="source($fluxcal,$passcals,$gaincals)" \
    line=wide,1,$wide | tee -a $log

# Plot The Selfcal Gains On Fluxcal Alongside The Adopted Gains:
rm -rf $fluxcal.sgains
mselfcal vis=pbcald.nw options=amp,apriori,noscale interval=$pcalint \
    line=wide,1,$wide out=$fluxcal.sgains select="source($fluxcal)"
gpcopy vis=pbcald.nw out=$fluxcal.sgains mode=merge
gpplt vis=$fluxcal.sgains yaxis=amp device=$fluxcal.sgains.ps/ps \
    nxy=5,3 yrange=0,3
if ($showplots == 'y') gv $fluxcal.sgains.ps &

#--------------------------------------------------------------------------
# 7. Produce dirty maps.
#--------------------------------------------------------------------------

domaps:

# Log Entry
echo " " | tee -a $log
echo "Producing dirty cubes..."
echo "------ imaged sources on `date`" >> $log

foreach line (co 13co)
# DETERMINE THE VELOCITY RANGE (DEFAULTS FOR EDGE)
    if ($line == 'co') then
        set linewin=$win12
        set winvec=`echo "$win12" | sed 's/,/ /g'`
        set red=$winvec[1]
        set blue=$winvec[$#winvec]
        uvspec vis=uvdata.mir select="window($blue)" axis=vel,amp log=tmp.uvspec1 \
            interval=10000 options=nobase,avall
        set chani=`wc -l tmp.uvspec1 | awk '{printf "%d", $1*7/8}'`
        set veli=`sed -n "${chani}p" tmp.uvspec1 | awk '{printf "%d", $1/20}' | awk '{print $1*20}'`
        uvspec vis=uvdata.mir select="window($red)" axis=vel,amp log=tmp.uvspec2 \
            interval=10000 options=nobase,avall
        set chanf=`wc -l tmp.uvspec2 | awk '{printf "%d", $1*1/8}'`
        set velf=`sed -n "${chanf}p" tmp.uvspec2 | awk '{printf "%d", $1/20}' | awk '{print $1*20}'`
        echo "The velocity range is $veli to $velf"
        set nmaps=`calc -i "1+($velf-$veli)/$delv"`
        set velf=`calc -i "$veli+$delv*($nmaps-1)"`
        echo "Making $nmaps maps from velocity $veli to $velf and spaced by $delv km/s" 
        set obsfreq=`grep 'starting freq' spectra.log | head -1 | awk '{print $9}'`
    else
        set linewin=$win13
        set winvec=`echo "$win13" | sed 's/,/ /g'`
        set red=$winvec[1]
        set blue=$winvec[$#winvec]
    endif

# MAP SOURCES WITH NATURAL WEIGHTING
    foreach source ($sourcelis)
        # -- generate visibility file
        rm -rf $source.$line.cal
        uvaver vis=pbcald.nw out=$source.$line.cal select="source($source),window($linewin)"
        # -- image the visibilities
        if ($line == 'co') then
            rm -rf $source.$line.map $source.$line.beam
            invert vis=$source.$line.cal map=$source.$line.map beam=$source.$line.beam \
                imsize=$imsize cell=$cell line=velocity,$nmaps,$veli,$delv,$delv \
                options=systemp,double,mosaic slop=1 sup=0 $exinvert | tee tmp.invert
            # -- report theoretical noise
            echo "${source}: Produced $nmaps maps from $veli to $velf spaced by $delv km/s" \
                >> $log
            echo -n "  " >> $log
            grep Theoretical tmp.invert >> $log
            echo -n "  " >> $log
            grep Number tmp.invert >> $log
            set nptng=`grep Number tmp.invert | awk '{print $4/3}'`
            # -- report actual noise
            rm -rf $source.$line.sen $source.$line.gain
            mossen in=$source.$line.map sen=$source.$line.sen gain=$source.$line.gain
            rm -rf $source.$line.norm $source.$line.sengain
            set minval=`histo in=$source.$line.sen region='image(1)'|grep Min|awk '{print $3}'`
            echo "Normalizing $source.$line.sen by $minval"
            maths exp="$source.$line.map*$minval/$source.$line.sen" out=$source.$line.norm
            maths exp="$source.$line.sen/($source.$line.gain*$minval)" out=$source.$line.sengain
            set nmin1 = `expr $nmaps - 1`
            set exec="histo in=$source.$line.norm"
            set box="box(-60,-60,60,60)"
            set crmsi=`$exec region="a,$box(1,2)" |grep Rms |awk '{printf "%f",$4}'`
            set crmsf=`$exec region="a,$box($nmin1,$nmaps)"|grep Rms|awk '{printf "%f",$4}'`
            set crms = `calc -f f7.5 "($crmsi+$crmsf)/2"`
            echo "  RMS estimates at field center: $crmsi $crmsf  Mean: $crms" | tee -a $log
            # -- Generate mask for cleaning
            rm -rf $source.$line.senmsk
            maths exp=$source.$line.sengain out=$source.$line.senmsk \
                mask="<$source.$line.sengain>.lt.$senlim"
            # -- Determine scaling factor for rmsfac
            set scale=`calc -f f7.3 "$crms/$minval"`
            echo "  Scaling factor for rmsfac: $scale" | tee -a $log
            puthd in=$source.$line.map/rms value=$crms type=real
            puthd in=$source.$line.map/scale value=$scale type=real
            puthd in=$source.$line.map/nptng value=$nptng type=integer
            # -- Plot normalized channel maps
            cgdisp in=$source.$line.norm,$source.$line.sengain slev=a,$crms,a,1 \
                levs1=2,4,6,8,10,12,14 levs2=$senlim region=$xyreg \
                device=$source.$line.dmaps.ps/ps options=full,3value,mirr,solneg2 \
                labtyp=arcsec nxy=5,4 type=contour csize=0.5,0.7
            if ($showplots == 'y') gv $source.$line.dmaps.ps &
        endif
    end # of loop over sources

end # of loop over line

if ($do_deconv != y) then
    ls -lt *.ps
    if ($purge == y) then
        goto dorm
    else
        goto end
    endif
else
    goto doclean
endif


#--------------------------------------------------------------------------
# 8. Deconvolve maps, if desired.
#--------------------------------------------------------------------------

doclean:

# Log Entry
echo " " | tee -a $log
echo "------ deconvolved maps using mossdi2 on `date`" >> $log
echo "Region defined by $senlim x min. noise; rmsfac: $rmsfac" >> $log

set line=co

# CLEAN maps with MOSSDI2:
foreach source ($sourcelis)
    set crms=`gethd in=$source.$line.map/rms`
    set scale=`gethd in=$source.$line.map/scale`
    set cutoff=`calc "$scale*$rmsfac"`
    rm -rf $source.$line.sdi
    mossdi2 map=$source.$line.map beam=$source.$line.beam out=$source.$line.sdi \
        cutoff=$cutoff region="mask($source.$line.senmsk)" niters=10000
    rm -rf $source.$line.psf
    mospsf beam=$source.$line.beam out=$source.$line.psf
    imfit in=$source.$line.psf object=beam region="a,box(-10,-10,10,10)" \
        | tee tmp.imfit
    set fitmaj=`grep 'Major axis' tmp.imfit | sed 's/.*://' | awk '{print $1}'`
    set fitmin=`grep 'Minor axis' tmp.imfit | sed 's/.*://' | awk '{print $1}'`
    set fitpa=`grep 'Position angle' tmp.imfit | sed 's/.*://' | awk '{print $1}'`
    rm -rf $source.$line.cm
    restor map=$source.$line.map beam=$source.$line.beam model=$source.$line.sdi \
        out=$source.$line.cm fwhm=$fitmaj,$fitmin pa=$fitpa
    rm -rf $source.$line.res
    restor map=$source.$line.map beam=$source.$line.beam model=$source.$line.sdi \
        out=$source.$line.res fwhm=$fitmaj,$fitmin pa=$fitpa mode=residual
    copyhd in=$source.$line.cm out=$source.$line.res items=bmaj,bmin,bpa
    set minval=`histo in=$source.$line.sen region='image(1)'|grep Min|awk '{print $3}'`
    rm -rf $source.$line.cmnorm
    maths exp="$source.$line.cm*$minval/$source.$line.sen" out=$source.$line.cmnorm
    rm -rf $source.$line.resnorm
    maths exp="$source.$line.res*$minval/$source.$line.sen" out=$source.$line.resnorm
    set maxplt=`histo in=$source.$line.cmnorm | grep Max | awk '{print $3*1.1}'`
    cgdisp in=$source.$line.cmnorm,$source.$line.cmnorm,$source.$line.sengain \
        slev=a,$crms,a,1 levs1=2,4,6,8,10,12,14 levs2=$senlim region=$xyreg \
        device=$source.$line.cmaps.ps/ps options=full,3value,mirr,solneg2 \
        labtyp=arcsec nxy=5,4 type=p,c,c csize=0.5,0.7 range=0,$maxplt,lin,1
    if ($showplots == "y") then
        gv $source.$line.cmaps.ps &
        echo "Plotting residuals for ${source}..."
        cgdisp in=$source.$line.resnorm device=/xs region=$xyreg type=p csize=0,0.8 \
            nxy=5,4 options=full,3val labtyp=arcsec range=0,0,lin,3
        xspause
    endif
    set rflux=`histo in=$source.$line.resnorm region="mask($source.$line.senmsk)" | grep Flux | awk '{print $6}'`
    set rintflux=`calc -f f8.2 "$rflux*$delv"`
    set cflux=`histo in=$source.$line.sdi | grep Flux | awk '{print $6}'`
    set cintflux=`calc -f f8.2 "$cflux*$delv"`
    echo "${source}: Used restoring beam of $fitmaj x $fitmin arcsec, pa $fitpa" \
        | tee -a $log
    echo "  Flux of norm. residuals in region: $rintflux Jy km/s" | tee -a $log
    echo "  Flux of deconv. model in region: $cintflux Jy km/s" | tee -a $log
    # --- peak brightness map (convert to mK)
    rm -rf $source.$line.pknorm
    moment in=$source.$line.cmnorm out=$source.$line.pknorm mom=-2
    set jyperk=`calc "2*1.3807e7/((30/$obsfreq)*206265)**2"`
    set bmaj=`gethd in=$source.$line.pknorm/bmaj format=arcsec`
    set bmin=`gethd in=$source.$line.pknorm/bmin format=arcsec`
    set bmarea=`calc "1.133*$bmaj*$bmin"`
    rm -rf $source.$line.tpeak
    maths exp="1000*<$source.$line.pknorm>/($bmarea*$jyperk)" out=$source.$line.tpeak
    puthd in=$source.$line.tpeak/bunit value='mK' type=ascii
    set crms2=`calc "1000*$crms/($bmarea*$jyperk)"`
    echo "The new rms is $crms2"
    puthd in=$source.$line.tpeak/rms value=$crms2
    rm -rf $source.$line.sengain1
    imsub in=$source.$line.sengain out=$source.$line.sengain1 region='image(1)'
    rm -rf $source.$line.tpeakmsk
    maths exp=$source.$line.tpeak out=$source.$line.tpeakmsk \
        mask="<$source.$line.sengain1>.lt.2"
    set maxplt=`histo in=$source.$line.tpeakmsk | grep Max | awk '{print $3*1.}'`
    cgdisp in=$source.$line.tpeak,$source.$line.tpeak,$source.$line.sengain1 \
        type=p,c,c slev=a,$crms2,a,1 levs1=4,5,6 levs2=$senlim cols1=4 cols2=1 \
        region=$xyreg device=${source}_$obsdate.$line.tpeak.ps/cps \
        options=full,solneg2,black,wedge,beambl labtyp=arcsec nxy=1,1 \
        csize=0.8,0.7 range=$crms2,$maxplt,lin,-8 lines=2,2,2
end

echo ""
ls -lt *.ps

#--------------------------------------------------------------------------
# 9. Output for HTML log.
#--------------------------------------------------------------------------

set line=co
rm -f reduclog.html
foreach source ($sourcelis)
    echo "<tr>" >> reduclog.html
    echo "<td>$obsdate</td>" >> reduclog.html
    echo "<td>$file:r</td>" >> reduclog.html
    echo "<td>$tauavg</td>" >> reduclog.html
    echo "<td>$rmpavg</td>" >> reduclog.html
    echo "<td>$fluxcal</td>" >> reduclog.html
    echo "<td>$fluxflux</td>" >> reduclog.html
    echo "<td>$bootfac</td>" >> reduclog.html
    echo "<td>$source</td>" >> reduclog.html
    uvindex vis=$source.$line.cal > tmp.uvindex
    set srctime=`grep 'Total obs' tmp.uvindex | awk '{print $5}'`
    echo "<td>$srctime</td>" >> reduclog.html
    set nptng=`gethd in=$source.$line.map/nptng`
    echo "<td>$nptng</td>" >> reduclog.html
    set trms=`gethd in=$source.$line.tpeak/rms`
    set trmsint=`echo $trms | awk '{printf "%.0f",$1}'`
    set tmax=`histo in=$source.$line.tpeakmsk | grep Max | awk '{printf "%.0f",$3}'`
    set snr=`histo in=$source.$line.tpeakmsk | grep Max | awk -v s=$trms '{printf "%.1f",$3/s}'`
    echo "<td align="\"right\"">$trmsint</td>" >> reduclog.html
    echo "<td align="\"right\"">$tmax</td>" >> reduclog.html
    echo "<td align="\"right\"">$snr</td>" >> reduclog.html
    if (-e ${source}_$obsdate.$line.tpeak.ps) then
        echo "<td><a target="\"_top\"" href="\"http://mmwave.astro.illinois.edu/carma/edge/${arr}_tracks_tpk/${source}_$obsdate.$line.tpeak.pdf\"">link</a></td>" >> reduclog.html
    else
        echo "<td>ERROR</td>" >> reduclog.html
    endif
    echo "<td>`date`</td>" >> reduclog.html
    echo "</tr>" >> reduclog.html
end

if ($purge == y) then
    goto dorm
else
    goto end
endif
    
#--------------------------------------------------------------------------
# 10. When satisfied, clean up files.
#--------------------------------------------------------------------------

dorm:
set echo

# initial fixes:
rm -rf uvdata.{mir,wins} uvedit.mir lineflg.mir pbcald.{mir,nw}
rm -rf uvputhd.mir uvredo.mir
#rm -rf *_c

# gain calibration / mapping calibrator:
# foreach gaincal ($gaincalis)
#     rm -rf $gaincal.sgains
# end

# mapping/deconvolution:
foreach source ($sourcelis)
    rm -rf $source.co.res $source.co.resnorm $source.co.norm 
    rm -rf $source.co.gain $source.co.sen* $source.co.psf
end

# Archive the peak temperature maps
foreach source ($sourcelis)
    if (-e ${source}_$obsdate.co.tpeak.ps) then
        if !(-d ../tpkplts) mkdir ../tpkplts
        echo "Archiving ${source}_$obsdate.co.tpeak.ps to ../tpkplts"
        epspdf ${source}_$obsdate.co.tpeak.ps
        echo "${source}_$obsdate" | enscript -B -f Courier-Bold@32 -o- | \
            ps2pdf - | pdftk ${source}_$obsdate.co.tpeak.pdf stamp - output output.pdf
        if (-e output.pdf) mv output.pdf ${source}_$obsdate.co.tpeak.pdf
        mv ${source}_$obsdate.co.tpeak.pdf ../tpkplts
    endif
end        

# compress files:
gzip -f *.ps
unset echo

end:
echo "Output to $log"
rm -f tmp.* *.tmp