#!/bin/csh -f

source /star/etc/login  >& /dev/null
source /star/etc/cshrc  >& /dev/null
kappa >& /dev/null

if (`alias kappa` == "") then
  echo ""
  echo "Fatal: Initialization error of kappa environment"
  echo ""
  exit
endif

set kap = $KAPPA_DIR


# ----------------------------------------------------------------------
# Command line arguments

if ($#argv < 3 || $#argv > 4) then
cat << _EOT_

Name: sp_merge - Frequency-merge two cubes with a DC adjustment of each
                 spectrum determined from the overlap.
Version: 1.0

Use:  ps_merge file1 file2 outfile [nodisplay]

Merge two cubes in frequency by adjusting for the DC offset in the
overlap region of the spectra before using wcsmosaic to combine them.
The DC levels are determined from the median in the overlap region of 
each spectrum. If there is no overlap between the cubes the median over
the full frequency range of each cube is used instead.

_EOT_
  exit
endif

set file1 = `echo "$1" | sed s/".sdf"/""/`
if !(-e "${file1}.sdf") then
  echo "*** ERROR *** Can not find file ${file1}.sdf"
  echo " "
  exit
endif

set file2 = `echo "$2" | sed s/".sdf"/""/`
if !(-e "${file2}.sdf") then
  echo "*** ERROR *** Can not find file ${file2}.sdf"
  echo " "
  exit
endif

set outfile = `echo "$3" | sed s/".sdf"/""/`

set display = 1
if ($#argv == 4) then
  set display = 0
endif

# ----------------------------------------------------------------------

echo " "
echo "*** Merge ${file1} ${file2} ***"

# Use a manual DC leveling and wcsmosaic rather than wcsalign and makemos

# Find the ranges in each cube
${kap}/ndftrace $file1 >& /dev/null
set lbnd    = `${kap}/parget lbound ndftrace | head -1`
set ubnd    = `${kap}/parget ubound ndftrace | head -1`
set zlo1    =  $lbnd[3]
set zhi1    =  $ubnd[3]
set flbnd   = `${kap}/parget flbnd ndftrace | head -1`
set fubnd   = `${kap}/parget fubnd ndftrace | head -1`
set pixsc   = `${kap}/parget fpixscale ndftrace | head -1`
set flo1    = `echo "scale=5;(($flbnd[3])+($pixsc[3]))" | bc`
set fhi1    = `echo "scale=5;(($fubnd[3])-($pixsc[3]))" | bc`
echo "Spectral window ${file1}: ${flo1}:${fhi1}"

${kap}/ndftrace $file2 >& /dev/null
set lbnd    = `${kap}/parget lbound ndftrace | head -1`
set ubnd    = `${kap}/parget ubound ndftrace | head -1`
set zlo2    =  $lbnd[3]
set zhi2    =  $ubnd[3]
set flbnd   = `${kap}/parget flbnd ndftrace | head -1`
set fubnd   = `${kap}/parget fubnd ndftrace | head -1`
set pixsc   = `${kap}/parget fpixscale ndftrace | head -1`
set flo2    = `echo "scale=5;(($flbnd[3])+($pixsc[3]))" | bc`
set fhi2    = `echo "scale=5;(($fubnd[3])-($pixsc[3]))" | bc`
echo "Spectral window ${file2}: ${flo2}:${fhi2}"

# Sort out the ordering of the cube boundaries and set the
# overlap range frange = flo:fhi accordingly.

if (`echo "1000*(($fhi1)-($flo1))" | bc | cut -d\. -f1` > 0) then
  if (`echo "1000*(($flo2)-($flo1))" | bc | cut -d\. -f1` > 0) then
    set flo = $flo2
    set fhi = $fhi1
  else
    set flo = $flo1
    set fhi = $fhi2
  endif
else
  if (`echo "1000*(($fhi2)-($fhi1))" | bc | cut -d\. -f1` > 0) then
    set flo = $fhi2
    set fhi = $flo1
  else
    set flo = $fhi1
    set fhi = $flo2
  endif
endif

set frange = "${flo}:${fhi}"
if (`echo "1000*(($fhi)-($flo))" | bc | cut -d\. -f1` < 0) then
  echo ""
  echo "*** WARNING *** no overlap between the input cubes:"
  echo "using the median over the whole spectrum instead"
  echo ""
  set frange = ""
endif
echo "Overlap: ${frange}" 

# Remove previous files
\rm ${file1}_dc.sdf ${file2}_dc.sdf >& /dev/null

# Use collapse to find the median for each pixel
#
echo "Determining median in overlap region of ${file1}"
${kap}/collapse ${file1}'(,,'$frange')' spmerge_temp1 estimator=median \
   axis=3 variance=true wlim=0.0
echo "Determining median in overlap region of ${file2}"
${kap}/collapse ${file2}'(,,'$frange')' spmerge_temp2 estimator=median \
   axis=3 variance=true wlim=0.0

# DC level adjustment is half the difference between the medians
${kap}/maths "(ib-ia)/2.0" ia=spmerge_temp1 ib=spmerge_temp2 \
     out=spmerge_temp3

# Grow the plane with the adjustment levels back into a cube and add
# to file1
echo "Adding DC level to ${file1}"
${kap}/manic spmerge_temp3 spmerge_temp4 axes='[1,2,0]' \
     lbound=${zlo1} ubound=${zhi1}
${kap}/add ${file1} spmerge_temp4 ${file1}_dc
set in1 = ${file1}_dc

# Grow the plane with the adjustment levels back into a cube and subtract
# from file2
echo "Subtracting DC level from ${file2}"
${kap}/manic spmerge_temp3 spmerge_temp4 axes='[1,2,0]' \
     lbound=${zlo2} ubound=${zhi2}
${kap}/sub ${file2} spmerge_temp4 ${file2}_dc
set in2 = ${file2}_dc

\rm -f spmerge_temp[1-4].sdf >& /dev/null

# Create combined spectrum by wcsmosaicing the leveled cubes
echo "Merging ${in1} and ${in2}"
${kap}/wcsmosaic in="'"${in1}","${in2}"'" ref=\! \
          method=nearest lbnd=\! ubnd=\! out=${outfile} \
          variance=true genvar=false wlim=0

echo " "
echo "RESULTS:"
echo "   DC adjusted input files: ${in1} ${in2}"
echo "   Combined spectrum: ${outfile}"

# ----------------------------------------------------------------------
# Display result

if ($display == 1) then

  echo "Displaying result for spectrum (~1,~1,)..."

  # Find extent final spectrum 
  ${kap}/ndftrace ${outfile} >& /dev/null
  set lbnd   = `${kap}/parget flbnd ndftrace | head -1`
  set ubnd   = `${kap}/parget fubnd ndftrace | head -1`
  set pixsc  = `${kap}/parget fpixscale ndftrace | head -1`

  set fstep  = $pixsc[3]
  set flo    = `echo "scale=5;(($lbnd[3])-50*($fstep))" | bc`
  set fhi    = `echo "scale=5;(($ubnd[3])+50*($fstep))" | bc`
  if (`echo "1000*(($fhi1)-($flo1))" | bc | cut -d\. -f1` < 0) then
    set temp = $flo
    set flo = $fhi
    set fhi = $flo
  endif

  # Get minmax from the input files
  ${kap}/histat ${file1}'(~1,~1,)' >& /dev/null
  set max1  = `${kap}/parget maximum histat | head -1`
  set min1  = `${kap}/parget minimum histat | head -1`
  histat ${file2}'(~1,~1,)' >& /dev/null
  set max2  = `${kap}/parget maximum histat | head -1`
  set min2  = `${kap}/parget minimum histat | head -1`

  # Find common Y-scale
  set max = $max1
  if (`echo "1000*(($max2)-($max1))" | bc | cut -d\. -f1` > 0) then
    set max = $max2
  endif
  set min = $min2
  if (`echo "1000*(($min2)-($min1))" | bc | cut -d\. -f1` > 0) then
    set min = $min1
  endif

  set yr = `echo "scale=5;(($max)-($min))" | bc`
  set ybot = `echo "scale=5;(($min)-0.33*($yr))" | bc`
  set ytop = `echo "scale=5;(($max)+0.33*($yr))" | bc`

  ${kap}/linplot $file1'(~1,~1,)' device=xw mode=histogram \
     xleft=$flo xright=$fhi ybot=${ybot} ytop=${ytop} \
     style="'colour(Lines)=cyan'" 
  ${kap}/linplot $file2'(~1,~1,)' device=xw mode=histogram \
     noclear style="'colour(Lines)=red'" 
  ${kap}/linplot ${outfile}'(~1,~1,)' device=xw mode=histogram \
     noclear style="'colour(Lines)=white'"

endif

echo " "

exit