#!/bin/bash
banner[1]='_______________________________________________________________________________'
banner[2]='                                                                               '
banner[3]='          HIRESCOLDENS = Compute High-Resolution Column Density Images         '
banner[4]='              Alexander Men`shchikov, SAp IRFU CEA Saclay, France              '
banner[5]='                                Version 1.140127                               '
banner[6]='_______________________________________________________________________________'
#
# Development initiated on 2014/12/24 by A.M. (alexander.menshchikov@cea.fr)
#___________________________________________________________________________________________________________________________________

CONDITIONS ()
{
  echo '#  _____________________________________________________________________________'
  echo '# |                                                                             '
  echo '# |         HIRESCOLDENS = Compute High-Resolution Column Density Images        '
  echo '# |             Alexander Men''shchikov, SAp IRFU CEA Saclay, France             '
  echo '# |_____________________________________________________________________________'
  echo '# |                                                                             '
  echo '# |    CONDITIONS OF USE (*1)                                                   '
  echo '# |                                                                             '
  echo '# | In what follows, "this software" refers to the Bash scripts GETSOURCES,     '
  echo '# | CONVOLVE, HIRESCOLDENS, IOSPEED, PREPAREOBS, RESAMPLE, & the associated     '
  echo '# | FORTRAN utilities CLEANBG, ELLIPSES, EXPANDA, FFTCONV, FINALCAT, FITFLUXES, '
  echo '# | FMEASURE, IMGSTAT, MODFITS, OPERATE, READHEAD, SEPARATE, SFINDER, SMEASURE, '
  echo '# | SPLITCUBE, and the library TOOLS, and "the author" refers to Alexander      '
  echo '# | Men'\''shchikov of the Service d'\''Astrophysique, IRFU, CEA Saclay, France.      '
  echo '# |                                                                             '
  echo '# | It is assumed that anyone using this code ("the user") has read, understood,'
  echo '# | and agreed to the following conditions of use:                              '
  echo '# |                                                                             '
  echo '# | 1. Distribution of this software shall remain the purview of the author. A  '
  echo '# |    user is free to share this code with co-workers and students, but if it  '
  echo '# |    is requested by a colleague not working directly with the user, the user '
  echo '# |    is asked to redirect such requests to: alexander.menshchikov@cea.fr      '
  echo '# |                                                                             '
  echo '# | 2. This software shall be used exclusively for education, research, non-    '
  echo '# |    profit, and non-military purposes. Specific written permission from the  '
  echo '# |    author must be obtained before any commercial use of this software is    '
  echo '# |    undertaken.                                                              '
  echo '# |                                                                             '
  echo '# | 3. The banners of this software, these conditions of use, and information   '
  echo '# |    about the author shall remain with this software and any descendent      '
  echo '# |    developed from and still based substantially upon this software.         '
  echo '# |                                                                             '
  echo '# | 4. The names of the institutions with which the author is or has been       '
  echo '# |    affiliated shall not be used to publicise any data and (or) results      '
  echo '# |    generated by this software. All findings and their interpretation are    ' 
  echo '# |    the opinions of the user and do not necessarily reflect those of the     '
  echo '# |    author nor the institutions with which the author is or has been         '
  echo '# |    affiliated.                                                              '
  echo '# |                                                                             '
  echo '# | The author makes no representations about the suitability of this software  '
  echo '# | for any purpose. Subject to the above conditions, this software are provided'
  echo '# | "as is" without any expressed or implied warranty.                          '
  echo '# |                                                                             '
  echo '# | Inquiries about the code, bug reports, constructive criticism, etc., can be '
  echo '# | directed to: alexander.menshchikov@cea.fr                                   '
  echo '# |                                                                             '
  echo '# |    CITATIONS AND ACKNOWLEDGEMENTS (*1)                                      '
  echo '# |                                                                             '
  echo '# | The algorithms used in this software are described in:                      '
  echo '# | Men'\''shchikov A. 2013, A&A, 560, A63                                         '
  echo '# | Palmeirim P., Andre Ph., Kirk J., et al. 2013, A&A, 550, A38                '
  echo '# | Men'\''shchikov A., Andre Ph., Didelon P., et al. 2012, A&A, 542, A81          '
  echo '# |                                                                             '
  echo '# | A very brief description of GETSOURCES has been published in:               '
  echo '# | Men'\''shchikov A., Andre Ph., Didelon P., et al. 2010, A&A, 518, L103         '
  echo '# |                                                                             '
  echo '# | It is requested that any publication reporting results obtained using this  '
  echo '# | software or any of its derivatives includes something like this:            '
  echo '# |                                                                             '
  echo '# | "Use of GETSOURCES (GETFILAMENTS), developed by A. Men'\''shchikov             '
  echo '# | at the SAp, IRFU, CEA Saclay, France, is hereby acknowledged." (*2)         '
  echo '# |_____________________________________________________________________________'
  echo '#                                                                               '
  echo '# (*1) Adapted from those written by David Clarke for the MHD code AZEuS.       '
  echo '# (*2) If only individual scripts or utilities (not the GETSOURCES script)      '
  echo '#      were used, substitute "GETSOURCES" with their names.                     '
}
#___________________________________________________________________________________________________________________________________
#
nparams=${#*}; commandline="$*"; scriptname=`basename "$0"`; script=`echo $scriptname | tr a-z A-Z`
scriptnameup=`echo $scriptname | tr a-z A-Z`; scriptnamelo=`echo $scriptname | tr A-Z a-z`; logfile='+log.'$scriptnameup
funo=$fun; fun=$scriptnameup
configname=$1; action=$2; number=$3; suffix=$4; nwavefin=$5; if [[ "$number" == "" ]]; then number=0; fi; verb='-verb1'
#___________________________________________________________________________________________________________________________________

chk_rc ()
{
# Checking return code and exiting if abnormal termination.
 
  local rc=$?
  if [[ "$rc" == "0" && "$1" != "" ]]; then rc=$1; fi
  if [[ "$rc" != "0" ]]
  then 
    where=$script; if [[ $fun != "" ]]; then where=$where': '$fun; fi; echo ' ERROR in '$where': RC='$rc': Aborted.'; echo
    exit "$rc"
  fi
}
#___________________________________________________________________________________________________________________________________

LEADZEROS ()
{
# Add $1 leading zeros to an integer given in $2 and return result in $3

  local zizi; local y; local numz=$1; zizi=$2

  if [[ $numz -eq 1 ]]; then if [[ $2 -lt 10 ]]; then zizi="0"$2  ; fi; fi
  if [[ $numz -eq 2 ]]; then if [[ $2 -lt 10 ]]; then zizi="00"$2 ; else if [[ $2 -lt 100  ]]; then zizi="0"$2 ; fi; fi; fi
  if [[ $numz -eq 3 ]]; then if [[ $2 -lt 10 ]]; then zizi='000'$2; else if [[ $2 -lt 100  ]]; then zizi='00'$2; else
                                                                         if [[ $2 -lt 1000 ]]; then zizi='0'$2 ; fi; fi; fi; fi
  y=\$"$3"; eval "$3=$zizi"
}
#_________________________________________________________________________________________________________________________________

POSITION ()
{
# Finding position of a substring in a string.

  local funo=$fun; fun='POSITION'
  if [[ $2 == "" ]]
  then echo; echo ' '$script': '$funo': '$fun': ERROR:'
    echo ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    echo ' INFO: Too few parameters: "'$1'" "'$2'"'
    echo ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    echo; fun=$funo; return 99
  fi
  local substr=$1; local string=$2 

  lensub=${#substr}; lenstr=${#string}; local dlen=$(( lenstr - lensub )); local i; returnposition=0
  if [[ $lensub -le $lenstr ]]
  then for (( i=0; i <= $dlen; i++ ))
  do if [[ "$substr" == "${string:$i:$lensub}" ]]; then returnposition=$(( i + 1 )); break; fi; done; fi

  fun=$funo; return $returnposition
}
#___________________________________________________________________________________________________________________________________

COMPARE_NUMBERS ()
{ 
# Comparing two either floating-point or integer numbers.

  local funo=$fun; fun='COMPARE_NUMBERS'
  if [[ $3 == "" ]]
  then echo; echo ' '$script': '$funo': '$fun': ERROR:'
    echo ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    echo ' INFO: Too few parameters: "'$1'" "'$2'"'
    echo ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    echo; fun=$funo; return 99
  fi
  if [[ $2 != ">" && $2 != "<" && $2 != "=" && $2 != "!=" && $2 != ">=" && $2 != "<=" ]]
  then echo; echo ' '$script': '$funo': '$fun': ERROR:'
    echo ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    echo ' INFO: Unsupported operator: "'$2'"'
    echo ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    echo; fun=$funo; return 99
  fi
  local num1=$1; local oper=$2; local num2=$3; local result=99; local i; local res1; local res2; local dot1; local dot2
  local sign1; local sign2
  if [[ ${num1:0:1} == "-" ]]; then sign1='-'; anum1=${num1:1}; else sign1=''; anum1=$num1; fi
  if [[ ${num2:0:1} == "-" ]]; then sign2='-'; anum2=${num2:1}; else sign2=''; anum2=$num2; fi
  'POSITION' '.' "$anum1"; dot1=$returnposition
  if [[ $dot1 -eq 0 ]]; then anum1=$anum1'.0'; 'POSITION' '.' "$anum1"; dot1=$returnposition; fi
  'POSITION' '.' "$anum2"; dot2=$returnposition
  if [[ $dot2 -eq 0 ]]; then anum2=$anum2'.0'; 'POSITION' '.' "$anum2"; dot2=$returnposition; fi
  local ddot=$(( dot1 - dot2 ))
  local len1=${#anum1}; local len2=${#anum2}; local dlen=$(( len1 - len2 ))
  local fra1=$(( len1 - dot1 )); local fra2=$(( len2 - dot2 )); local dfra=$(( fra1 - fra2 ))
  if [[ $ddot -lt 0 ]]; then for (( i=1; i <= -$ddot; i++ )); do anum1='0'$anum1; done; fi
  if [[ $ddot -gt 0 ]]; then for (( i=1; i <=  $ddot; i++ )); do anum2='0'$anum2; done; fi
  if [[ $dfra -lt 0 ]]; then for (( i=1; i <= -$dfra; i++ )); do anum1=$anum1'0'; done; fi
  if [[ $dfra -gt 0 ]]; then for (( i=1; i <=  $dfra; i++ )); do anum2=$anum2'0'; done; fi
  result='0'
  if [[ ${oper:0:1} == ">"  ]]
  then 
    if [[ $sign1 == ""  && $sign2 == "-" ]]; then result='1'; fi
    if [[ $sign1 == ""  && $sign2 == ""  && $anum1 > $anum2 ]]; then result='1'; fi
    if [[ $sign1 == "-" && $sign2 == "-" && $anum1 < $anum2 ]]; then result='1'; fi
    if [[ $oper == ">=" ]]; then if [[ $anum1 == $anum2 ]]; then result='1'; fi; fi
  fi
  if [[ ${oper:0:1} == "<"  ]]
  then 
    if [[ $sign1 == "-" && $sign2 == ""  ]]; then result='1'; fi
    if [[ $sign1 == ""  && $sign2 == ""  && $anum1 < $anum2 ]]; then result='1'; fi
    if [[ $sign1 == "-" && $sign2 == "-" && $anum1 > $anum2 ]]; then result='1'; fi
    if [[ $oper == "<=" ]]; then if [[ $anum1 == $anum2 ]]; then result='1'; fi; fi
  fi
  if [[ $oper == "="  && $anum1 == $anum2 ]]; then result='1'; fi
  if [[ $oper == "!=" && $anum1 != $anum2 ]]; then result='1'; fi

  fun=$funo; return $result
}
#___________________________________________________________________________________________________________________________________

# Read and store the header (version) information from the very beginning of this script.

l=0
{
while read linefromthisfile
do
  if [[ ${linefromthisfile:0:7} == "banner[" && $l -le 7 ]]
  then
    l=$(( l + 1 ))
    len=${#linefromthisfile}; last=$(( len - 12 ))
    version[l]="${linefromthisfile:11:$last}"
    'POSITION' ' Version' "${version[l]}"; ncv=$returnposition
    if [[ $ncv -gt 0 ]]; then versnum=${version[l]:ncv:16}; fi
  fi
  if [[ $l -eq 7 ]]; then break; fi
done
} < $0

headerinfo[1]='#  ___________________________________________________________________________________________________________________________________'
headerinfo[2]='# |                                                                                                                                   '
headerinfo[3]='# | HIRESCOLDENS = Compute High-Resolution Column Density Images = '$versnum' = Alexander Men'\''shchikov, SAp IRFU CEA Saclay        '
headerinfo[4]='# |___________________________________________________________________________________________________________________________________'

if [[ "$*" == ":" ]]
then
  echo "${headerinfo[3]:3:126}"
  exit 0
fi

'CONDITIONS' >> '+log.CONDITIONS'

{
if [[ $verb == "" ]]; then verb='-verb2'; fi

if [[ $verb != "-verb0" ]]
then
  for (( b=1; b <= 6; b++ ))
  do echo "${banner[b]}"; done
fi

if [[ ! -e $configname || ($action == "combine" && $suffix == '') || \
      ($action != "convolve" && $action != "prepare" && $action != "fitfluxes" && $action != "combine" && $action != "complete") ]]
then
  echo
  echo ' USAGE: '$scriptname' <configname> [<action>] [<number>] [<prefix>|<suffix>] [<nwave>]'
  echo
  echo '   where <action> is {convolve|prepare|fitfluxes|combine|complete}'
  echo '         <number> is a serial identifier of the input file names'
  echo '         <prefix> is {thin|gray} - the fitting model (thinbody, graybody)'
  echo '         <suffix> is {coldens|tempers|chisqnu} - the type of images to combine'
  echo '          <nwave> is the number of a waveband (from the config file), at whose'
  echo '                     resolution the output images will be created'
  echo
  echo ' This is a script that processes a multi-waveband set of images for deriving'
  echo ' column density and temperature images: prepares them by cross-convolving to'
  echo ' the resolutions of each waveband, creates catalogs of the images, shows command'
  echo ' lines to be used for fitting with FITFLUXES, and creates high-resolution column'
  echo ' density images from the column densities obtained with different resolutions'
  echo ' (using the algorithm described by Palmeirim et al. 2013).'
  echo
  if [[ $configname == '' || \
       ($action != '' && $action != "convolve" && $action != "prepare" && $action != "fitfluxes" && \
                         $action != "combine"  && $action != "complete") ]]
  then exit
  fi
  if [[ ! -e $configname ]]
  then
    echo "#"                                               > $configname
    echo "# $scriptnameup: CONFIGURATION PARAMETERS FOR FITFLUXES" >> $configname
    echo "#"                                              >> $configname
    echo "# PM  WAVE  BEAM  OFFSET  IMAGE  {+/-, microns, arcsecs, MJy/sr, files}" >> $configname
    echo "#"                                              >> $configname
    echo "   -   070   8.4  0.0000  'image.070.obs.fits'" >> $configname
    echo "   +   160  13.5  0.0000  'image.160.obs.fits'" >> $configname
    echo "   +   250  18.2  0.0000  'image.250.obs.fits'" >> $configname
    echo "   +   350  24.9  0.0000  'image.350.obs.fits'" >> $configname
    echo "   +   500  36.3  0.0000  'image.500.obs.fits'" >> $configname
    
    echo ' Default configuration file '\'$configname\'' has been saved.'; echo
  else
    if [[ $action == "combine" && $suffix == '' ]]
    then
      echo ' Specify suffix for this action.'; echo
    else
      echo ' Configuration file '\'$configname\'' already exists: choose an action.'; echo
    fi
  fi
  exit 0
fi

if [[ $action == "complete" && ($nwavefin == '' || $suffix == '') ]]
then
  echo; echo ' ERROR: Specify all parameters for this action (no optional parameters).'; echo
  exit 99
else
  prefix=$suffix
fi

if [[ $action != 'combine' && $action != '' ]]  # Read configuration file
then
  nwaves=0
  beam[nwaves]='0.0'
##  echo; echo ' Reading '$configname'...'
  {
  while read param1 param2 param3 param4 param5
  do
##    echo $param1 $param2 $param3 $param4 $param5

    if [[ ${param1:0:1} != "#" && ${param1:0:1} != "!" ]]
    then
      nwaves=$(( nwaves + 1 ))
      nwavesm1=$(( nwaves - 1 ))
      sign[nwaves]=$param1; wave[nwaves]=$param2; beam[nwaves]=$param3; offset[nwaves]=$param4; image[nwaves]=$param5
      image[nwaves]=${image[nwaves]/.fits/}; oimage[nwaves]=`basename "${image[nwaves]}"`
      image[nwaves]=${image[nwaves]}'.fits'; oimage[nwaves]=${oimage[nwaves]}'.fits'
      if [[ ${sign[nwaves]} != '+' && ${sign[nwaves]} != '-' && $action == "prepare" ]]
      then
        echo; echo ' '$scriptnameup': ERROR: 1st column in '\'$configname\'' is not '\'+\'' or '\'-\''.'; echo; exit 99
      fi
      if [[ ! -e "${image[nwaves]}" && $action != "complete" ]]
      then
        echo; echo ' '$scriptnameup': ERROR: File '\'${image[nwaves]}\'' not found.'; echo; exit 99
      fi
      'COMPARE_NUMBERS' ${beam[nwaves]} '<' ${beam[nwavesm1]}; rslt=$?
      if [[ $rslt -eq 1 ]]
      then 
        echo; echo ' '$scriptnameup': ERROR: Beams in '\'$configname\'' not increasing monotonically.'; echo; exit 99
      fi 
    fi
  done
  } < $configname
fi

# CONVOLVE: Resample *purely model* images to new pixel and convolve them to corresponding Herschel resolutions.

if [[ $action == 'convolve' ]]
then
  newpixel=2
  'POSITION' ".con." "$configname"; concfg=$returnposition
  if [[ $concfg -eq 0 ]]
  then echo; echo ' '$scriptnameup': ERROR: Special config file (*.con.cfg) for convolution not found.'; echo; exit 99
  fi
  if [[ `echo *'.p2as.'*` != '*.p2as.*' ]]; then \rm *'.p2as.'*; chk_rc; fi
    
  for (( i=1; i <= $nwaves; i++ ))
  do
# Make profiles of original model images after resampling to a 2" pixel and convolution.

    imagei=${image[i]/.fits/}
    oimagei=${oimage[i]/.fits/}
    'modfits' min 0 y $imagei -o $imagei'.p2as' $verb; chk_rc
    'resample' $imagei'.p2as' $newpixel 101 101; chk_rc
    'modfits' key CDELT1 -$newpixel $imagei'.p2as' -o $imagei'.p2as'; chk_rc
    'modfits' key CDELT2  $newpixel $imagei'.p2as' -o $imagei'.p2as'; chk_rc
    'modfits' key CD1_1  -$newpixel $imagei'.p2as' -o $imagei'.p2as'; chk_rc
    'modfits' key CD2_2   $newpixel $imagei'.p2as' -o $imagei'.p2as'; chk_rc
    'operate' $imagei'.p2as' h $imagei'.p2as.resamp' -o $imagei'.p2as' $verb; chk_rc
    \rm $imagei'.p2as.resamp.fits'; chk_rc
    'convolve' 0 ${beam[i]} $imagei'.p2as' $oimagei'.p2as.r'${beam[i]/./p} '-verb1'; chk_rc
    'modfits' min 0 y $oimagei'.p2as.r'${beam[i]/./p} -o $oimagei'.p2as.r'${beam[i]/./p} $verb; chk_rc
    'modfits' proavg '51,51,100,51' $oimagei'.p2as.r'${beam[i]/./p} $verb; chk_rc

# Make profiles of original model coldens images after resampling to a 2" pixel and convolution.

    imagei='runtem.cdn'
    oimagei='runtem.cdn'
    'modfits' min 0 y $imagei -o $imagei'.p2as' $verb; chk_rc
    'modfits' div 1e20 $imagei'.p2as' -o $imagei'.p2as' $verb; chk_rc
    'resample' $imagei'.p2as' $newpixel 101 101; chk_rc
    'modfits' key CDELT1 -$newpixel $imagei'.p2as' -o $imagei'.p2as'; chk_rc
    'modfits' key CDELT2  $newpixel $imagei'.p2as' -o $imagei'.p2as'; chk_rc
    'modfits' key CD1_1  -$newpixel $imagei'.p2as' -o $imagei'.p2as'; chk_rc
    'modfits' key CD2_2   $newpixel $imagei'.p2as' -o $imagei'.p2as'; chk_rc
    'operate' $imagei'.p2as' h $imagei'.p2as.resamp' -o $imagei'.p2as' $verb; chk_rc
    \rm $imagei'.p2as.resamp.fits'; chk_rc
    'modfits' mul 1e20 $imagei'.p2as' -o $imagei'.p2as' $verb; chk_rc
    'convolve' 0 ${beam[i]} $imagei'.p2as' $oimagei'.p2as.r'${beam[i]/./p} '-verb1'; chk_rc
    'modfits' min 0 y $oimagei'.p2as.r'${beam[i]/./p} -o $oimagei'.p2as.r'${beam[i]/./p} $verb; chk_rc
    'modfits' proavg '51,51,100,51' $oimagei'.p2as.r'${beam[i]/./p} $verb; chk_rc

# Make profiles of original model coldens images after convolution.

    'modfits' proavg '151,151,300,151' $imagei $verb; chk_rc
    'convolve' 0 ${beam[i]} $imagei $oimagei'.r'${beam[i]/./p} '-verb1'; chk_rc
    'modfits' min 0 y $oimagei'.r'${beam[i]/./p} -o $oimagei'.r'${beam[i]/./p} $verb; chk_rc
    'modfits' proavg '151,151,300,151' $oimagei'.r'${beam[i]/./p} $verb; chk_rc
  done
  imagei='runtem.cdn'
  'modfits' proavg '51,51,100,51' $imagei'.p2as' $verb; chk_rc
fi

# PREPARE: Add offsets, cross-convolve to various resolutions, and create configuration files for FITFLUXES.

if [[ $action == 'prepare' ]]
then
  nwavesm1=$(( nwaves - 1 ))
  
  for (( j=1; j <= $nwaves; j++ ))
  do
    'modfits' add ${offset[j]} ${image[j]} -o ${oimage[j]/.fits/}'.r'${wave[j]} $verb; chk_rc
    'modfits' min 0 y ${oimage[j]/.fits/}'.r'${wave[j]} -o ${oimage[j]/.fits/}'.r'${wave[j]} $verb; chk_rc
  done
  
  for (( i=$nwaves; i >= 1; i-- ))
  do
    im1=$(( i - 1 ))
    for (( j=1; j <= $im1; j++ ))
    do
      'modfits' add ${offset[j]} ${image[j]} -o ${oimage[j]/.fits/}'.r'${wave[i]} $verb; chk_rc
      'convolve' ${beam[j]} ${beam[i]} ${oimage[j]/.fits/}'.r'${wave[i]} ${oimage[j]/.fits/}'.r'${wave[i]} '-verb1'; chk_rc
      'modfits' min 0 y ${oimage[j]/.fits/}'.r'${wave[i]} -o ${oimage[j]/.fits/}'.r'${wave[i]} $verb; chk_rc
    done
  done
  
  for (( i=$nwaves; i >= 2; i-- ))
  do
    im1=$(( i - 1 ))
    if [[ ${sign[i]} == "+" &&  ${sign[im1]} == "+" ]]
    then
      config=${configname/.cfg/}'.r'${wave[i]}'.cfg'
      echo "#"                                                     > $config
      echo "# CONFIGURATION FOR FITFLUXES (CREATED BY HIRESCOLDENS)" >> $config
      echo "#"                                                    >> $config
      echo "#  WAVE  BEAM  OFFSET  IMAGE (WITH OFFSET ADDED)  {+/-, microns, arcsecs, MJy/sr, files}" >> $config
      echo "#"                                                    >> $config
      for (( j=1; j <= $i; j++ ))
      do
        lenw=`expr "${wave[j]}" : '.*'`
        if [[ $lenw -lt 4 ]]; then wavex='0'${wave[j]}; else wavex=${wave[j]}; fi
        lenb=`expr "${beam[i]}" : '.*'`
        if [[ $lenb -lt 4 ]]; then beamx=' '${beam[i]}; else beamx=${beam[i]}; fi
        echo "   "$wavex"  "$beamx"  "${offset[j]}"  '"${oimage[j]/.fits/}.r${wave[i]}'.fits'"'" >> $config
      done
    fi
  done
fi

# FITFLUXES: Proceed to the pixel-by-pixel fitting of images.
  
if [[ $action == 'fitfluxes' ]]
then
  parts=${number/-/}
  if [[ $parts -eq 0 || $parts -gt 1 ]]
  then
    echo
    echo ' Here are the commands to create high-resolution column density images:'
  fi
  for (( i=$nwaves; i >= 1; i-- ))
  do
    config=${configname/.cfg/}'.r'${wave[i]}'.cfg'
    if [[ $number -eq -1 ]]
    then
      if [[ `echo $config.0?'.log'` != $config'.0?.log' ]]; then \rm $config.0?'.log'; chk_rc; fi
      if [[ `echo $config.0?'.chisqnu'*` != $config'.0?.chisqnu*' ]]; then \rm $config.0?'.chisqnu'*; chk_rc; fi
      if [[ `echo $config.0?'.coldens'*` != $config'.0?.coldens*' ]]; then \rm $config.0?'.coldens'*; chk_rc; fi
      if [[ `echo $config.0?'.tempers'*` != $config'.0?.tempers*' ]]; then \rm $config.0?'.tempers'*; chk_rc; fi
    fi
    im1=$(( i - 1 ))
    if [[ ${sign[i]} == "+" && ${sign[im1]} == "+" ]]
    then
      echo; echo ' COLUMN DENSITY IMAGE AT '${wave[i]}' um RESOLUTION:'
      im=$i
      for (( k=1; k <= $im1; k++ ))
      do
        if [[ ${sign[k]} == "+" ]]
        then
          im=$(( im - 1 ))
          allwaves=''
          allerrs=''
          for (( j=1; j <= $i; j++ ))
          do
            if [[ $j -lt $im ]]
            then allwaves=$allwaves',-'$j':'${wave[j]}
            else allwaves=$allwaves','$j':'${wave[j]}
            fi
            allerrs=$allerrs'1:0.2,'
          done
          allwaves=${allwaves:1}
          len=`expr "$allerrs" : '.*'`
          lenm1=$(( len - 1 ))
          allerrs=${allerrs:0:lenm1}
          if [[ $parts -gt 1 ]]
          then
            for (( part=1; part <= $parts; part++ ))
            do
              echo ' fitfluxes' 0 $part $parts $allwaves $allerrs thinbody 0:1,0:1,2:0 1 10 $config '-verb0'
            done
          else
            if [[ $parts -eq 0 ]]
            then
              echo ' fitfluxes' 0 0 0 $allwaves $allerrs thinbody 0:1,0:1,2:0 1 10 $config '-verb0'
###              echo ' fitfluxes' 0 0 0 $allwaves $allerrs graybody 0:1,0:1,2:0,0:1 1 10 $config '-verb0'
            else
              'fitfluxes' 0 0 0 $allwaves $allerrs thinbody 0:1,0:1,2:0 1 10 $config '-verb0'; chk_rc
        
###              if [[ -e 'thin.'$config'.00.coldens.fits' ]]
###              then
###                'modfits' proavg '51,51,100,51' 'thin.'$config'.00.coldens' $verb; chk_rc
###                'modfits' proavg '51,51,100,51' 'thin.'$config'.00.tempers' $verb; chk_rc
###              fi
###              if [[ -e 'thin.'$config'.01.coldens.fits' ]]
###              then
###                'modfits' proavg '51,51,100,51' 'thin.'$config'.01.coldens' $verb; chk_rc
###                'modfits' proavg '51,51,100,51' 'thin.'$config'.01.tempers' $verb; chk_rc
###              fi
###              if [[ -e 'thin.'$config'.02.coldens.fits' ]]
###              then 
###                'modfits' proavg '51,51,100,51' 'thin.'$config'.02.coldens' $verb; chk_rc
###                'modfits' proavg '51,51,100,51' 'thin.'$config'.02.tempers' $verb; chk_rc
###              fi
        
###              'fitfluxes' 0 0 0 $allwaves $allerrs graybody 0:1,0:1,0:1,2:0 1 10 $config '-verb0'; chk_rc
        
###              if [[ -e 'gray.'$config'.00.coldens.fits' ]]
###              then
###                'modfits' proavg '51,51,100,51' 'gray.'$config'.00.coldens' $verb; chk_rc
###                'modfits' proavg '51,51,100,51' 'gray.'$config'.00.tempers' $verb; chk_rc
###              fi
###              if [[ -e 'gray.'$config'.01.coldens.fits' ]]
###              then
###                'modfits' proavg '51,51,100,51' 'gray.'$config'.01.coldens' $verb; chk_rc
###                'modfits' proavg '51,51,100,51' 'gray.'$config'.01.tempers' $verb; chk_rc
###              fi
###              if [[ -e 'gray.'$config'.02.coldens.fits' ]]
###              then 
###                'modfits' proavg '51,51,100,51' 'gray.'$config'.02.coldens' $verb; chk_rc
###                'modfits' proavg '51,51,100,51' 'gray.'$config'.02.tempers' $verb; chk_rc
###              fi
            fi
          fi
        fi
      done
    fi
  done
  if [[ $parts -eq 0 || $parts -gt 1 ]]
  then
    echo
    echo ' Start the above commands in parallel on different processors.'
    if [[ $parts -gt 1 ]]
    then
      echo ' When finished, sum the results up in combined final images at respective'
      echo ' resolutions and create the high-resolution column density by executing:'
      echo
      echo ' hirescoldens' ${configname/.cfg/.r500.cfg} combine 00 coldens
      echo ' hirescoldens' ${configname/.cfg/.r500.cfg} combine 01 coldens
      echo ' hirescoldens' ${configname/.cfg/.r500.cfg} combine 02 coldens
      echo ' hirescoldens' ${configname/.cfg/.r350.cfg} combine 00 coldens
      echo ' hirescoldens' ${configname/.cfg/.r350.cfg} combine 01 coldens
      echo ' hirescoldens' ${configname/.cfg/.r250.cfg} combine 00 coldens
      echo
    fi
    echo ' Finally, create three versions of high-resolution column densities'
    echo ' (at 250 um resolution) by executing:'
    echo
    echo ' hirescoldens '$configname' complete 00 thin 3'
    echo ' hirescoldens '$configname' complete 01 thin 3'
    echo ' hirescoldens '$configname' complete 02 thin 3'
###    echo ' hirescoldens '$configname' complete 00 gray'
###    echo ' hirescoldens '$configname' complete 01 gray'
###    echo ' hirescoldens '$configname' complete 02 gray'
    echo
    echo ' Compare the three versions of high-resolution column densities and'
    echo ' use the best one. It is recommended to use the version labeled "01".'
  else
    'hirescoldens' $configname complete 00 thin; chk_rc
    'hirescoldens' $configname complete 01 thin; chk_rc
    'hirescoldens' $configname complete 02 thin; chk_rc
###    'hirescoldens' $configname complete 00 gray; chk_rc
###    'hirescoldens' $configname complete 01 gray; chk_rc
###    'hirescoldens' $configname complete 02 gray; chk_rc
###    hires='18p2'
###    if [[ -e $configname'.00.coldens.r'$hires'as.fits' ]]
###    then 'modfits' proavg '51,51,100,51' $configname'.00.coldens.r'$hires'as' $verb; chk_rc; fi
###    if [[ -e $configname'.01.coldens.r'$hires'as.fits' ]]
###    then 'modfits' proavg '51,51,100,51' $configname'.01.coldens.r'$hires'as' $verb; chk_rc; fi
###    if [[ -e $configname'.02.coldens.r'$hires'as.fits' ]]
###    then 'modfits' proavg '51,51,100,51' $configname'.02.coldens.r'$hires'as' $verb; chk_rc; fi
  fi
fi

# COMBINE: Combine partial column density images.

if [[ $action == 'combine' ]]
then
  echo; echo ' COMBINING PARTIAL COLUMN DENSITY IMAGES...'

  nbeg='000'; nend='000'
  image=$configname'.'$number
  imagesum=$image'.'$suffix

  for (( j=1; j <= 999; j++ ))
  do
    'LEADZEROS' 2 $j zj
    if [[ $nbeg == '000' && -e $image'.'$zj'.'$suffix'.fits' ]]; then nbeg=$j; fi
  done
  for (( j=999; j >= 1; j-- ))
  do
    'LEADZEROS' 2 $j zj
    if [[ $nend == '000' && -e $image'.'$zj'.'$suffix'.fits' ]]; then nend=$j; fi
  done
  
  'LEADZEROS' 2 $nbeg znbeg; 'LEADZEROS' 2 $nend znend

  'modfits' multiply 0 $image'.'$znbeg'.'$suffix -o $imagesum '-verb0'; chk_rc
  
  for (( j=$nbeg; j <= $nend; j++ ))
  do
    'LEADZEROS' 2 $j zj
    if [[ -e $image'.'$zj'.'$suffix'.fits' ]]
    then
      'operate' $image'.'$zj'.'$suffix + $imagesum -o $imagesum '-verb1'; chk_rc
    else
      echo; echo ' '$scriptnameup': ERROR: File '\'$image'.'$zj'.fits'\'' not found.'; echo; exit 99
    fi
  done
  echo; echo ' COMBINED IMAGES WITH SERIAL IDENTIFIERS FROM '$znbeg' TO '$znend
fi

# COMPLETE: Create high-resolution column density image.

if [[ $action == 'complete' ]]
then
  nwavesm1=$(( nwaves - 1 ))
###  nwavefin=$(( nwaves - 1 ))
###  if [[ $beammin == '250' ]]; then nwavefin=$(( nwaves - 2 )); fi
###  if [[ $beammin == '160' ]]; then nwavefin=$(( nwaves - 3 )); fi
###  if [[ $beammin == '100' ]]; then nwavefin=$(( nwaves - 4 )); fi

  echo; echo ' CREATING AN IMAGE OF HIGHER-RESOLUTION ('${beam[nwavefin]}' ARCSEC) COLUMN DENSITIES...'

  len=`expr "$number" : '.*'`
  if [[ ${number:0:1} == '0' && $len -eq 2 ]]
  then number=${number:1}
  fi
  'LEADZEROS' 1 $number znumber
  coldenslo=$prefix'.'${configname/.cfg/}'.r'${wave[nwaves]}'.cfg.'$znumber'.coldens'
  coldensout=$prefix'.'$configname'.'$znumber'.coldens'
  hires=${beam[nwaves]}

  echo; echo ' COPYING: '$coldenslo' -> '$coldensout'...'
  
  'modfits' mul 1 $coldenslo -o $coldensout '-verb0'; chk_rc

  for (( i=$nwavesm1; i >= $nwavefin; i-- ))
  do
    ip1=$(( i + 1 ))
    largerbeam=$(echo "scale=20; sqrt((${beam[i]})^2+(${beam[i]}/3)^2)" | bc); chk_rc
    largerbeam=${largerbeam:0:6}

    for (( j=$number; j >= 0; j-- ))
    do
      'LEADZEROS' 1 $j zj
      coldens=$prefix'.'${configname/.cfg/}'.r'${wave[i]}'.cfg.'$zj'.coldens'
      if [[ -e $coldens'.fits' ]]
      then break
      fi
    done
    if [[ $j -lt 0 && ! -e $coldens'.fits' ]]
    then 
      echo; echo ' WARNING: File '\'$coldens'.fits'\'' not found... waveband '${wave[i]}' um ignored.'
      exit 99
    else
      echo; echo ' CONVOLVING: '$coldens' FROM '${beam[i]}' TO '${beam[ip1]}' ARCSEC RESOLUTION...'
  
      'convolve' ${beam[i]} ${beam[ip1]} $coldens $coldens'.r'${beam[ip1]/./p} '-verb0'; chk_rc

###      'convolve' ${beam[i]} $largerbeam $coldens $coldens'.r'${largerbeam/./p} '-verb0'; chk_rc
###      echo; echo ' CONVOLVING: '$coldens'.r'${largerbeam/./p}' FROM '$largerbeam' TO '${beam[ip1]}' ARCSEC RESOLUTION...'
###      'convolve' $largerbeam ${beam[ip1]} $coldens'.r'${largerbeam/./p} $coldens'.r'${beam[ip1]/./p} '-verb0'; chk_rc
###      echo; echo ' DIFFERENCING: '$coldens'.r'${largerbeam/./p}' - '$coldens'.r'${beam[ip1]/./p}' -> '$coldens'.delta...'
###      'operate' $coldens'.r'${largerbeam/./p} - $coldens'.r'${beam[ip1]/./p} $coldens'.delta' '-verb0'; chk_rc  ### $verb

      echo; echo ' DIFFERENCING: '$coldens' - '$coldens'.r'${beam[ip1]/./p}' -> '$coldens'.delta...'

      'operate' $coldens - $coldens'.r'${beam[ip1]/./p} -o $coldens'.delta' '-verb0'; chk_rc  ### $verb

      echo; echo ' ADDING DIFFERENCE: '$coldens'.delta + '$coldensout' -> '$coldensout'...'

      'operate' $coldensout + $coldens'.delta' -o $coldensout '-verb0'; chk_rc  ### $verb

      hires=${beam[i]}
    fi
  done

  echo; echo ' REMOVING NEGATIVES: '$coldensout' -> '$coldensout'.r'${hires/./p}'as...'

  'modfits' min 0 y $coldensout -o $coldensout'.r'${hires/./p}'as' '-verb0'; chk_rc   ### $verb

  echo; echo ' OUTPUT COLUMN DENSITY IMAGE HAS A RESOLUTION OF '$hires' ARCSEC'

  \rm $coldensout'.fits'; chk_rc
fi

cdate=`date +%d\ %b\ %Y\ %a\ %H:%M:%S\ %Z`
echo
echo ' '$scriptnameup' DONE. '$cdate
} | tee -a $logfile
exit
#_________________________________________________________________________________________________________________________________
