#!/bin/bash 
banner[1]='_______________________________________________________________________________'
banner[2]='                                                                               '
banner[3]='                                   GETIMAGES                                   '
banner[4]='              • Background Derivation & Image Flattening Method •              '
banner[5]='              Alexander Men’shchikov, DAp IRFU CEA Saclay, France              '
banner[6]='                                Version 2.180830                               '
banner[7]='_______________________________________________________________________________'
#
# Development initiated on 2008/07/20 by A.M. (alexander.menshchikov@cea.fr)
#___________________________________________________________________________________________________________________________________

HISTORY ()
{
  for (( b=1; b <= 7; b++ )); do print '#'"${banner[b]}"; done
  print '#'
  print '#                               HISTORY OF CHANGES'
  print '#                               (only last months)'
  print '#'
  print '# ------------------------------ VERSION 2.180830 ------------------------------'
  print '#'
  print '# [CHA] Modified FLATTEN_IMAGES by adding a second iteration only at the largest'
  print '#       scale, in both the background derivation and image flattening parts,'
  print '#       in conjunction with minimization of the median-filtered image. This'
  print '#       considerably improves both the background and scaling images, making it'
  print '#       both lower and smoother.'
  print '#'
  print '# ------------------------------ VERSION 2.180830 ------------------------------'
  print '#'
  print '# [CHA] Few minor usability improvements in FLATTEN_IMAGES. Completely removed'
  print '#       all variables and references to iterations, and modified config file.'
  print '#'
  print '# ------------------------------ VERSION 2.180825 ------------------------------'
  print '#'
  print '# [FIX] Modified FLATTEN_IMAGES by increasing initial median filtering radius'
  print '#       for stdev images equal to two beam sizes, instead of a single beam size'
  print '#       (introduced in v2.180727) that led to artificial peaks in the flattened'
  print '#       images, in some cases.'
  print '#'
  print '# ------------------------------ VERSION 2.180818 ------------------------------'
  print '#'
  print '# [CHA] GETIMAGES now defaults to a single pass of median filtering. The'
  print '#       original approach of iterating background-subtracted image works well'
  print '#       for detecting compact sources, but it does not allow sufficiently full'
  print '#       (deep) reconstruction of filaments. This leads to incomplete subtraction'
  print '#       of filamentary structures when detecting sources, as well as incorrect'
  print '#       filament profiles, due to oversubtraction of (overestimated) background.'
  print '#'
  print '# ------------------------------ VERSION 2.180817 ------------------------------'
  print '#'
  print '# [CHA] In FLATTEN_IMAGES, uniform number of 3 iterations on all spatial scales'
  print '#       is now done only on the 2nd and larger spatial scales.'
  print '#'
  print '# ------------------------------ VERSION 2.180816 ------------------------------'
  print '#'
  print '# [CHA] Modified FLATTEN_IMAGES making a uniform number of 3 iterations on all'
  print '#       spatial scales, when median-filtering images. Before processing, total'
  print '#       execution time of GETIMAGES is estimated and displayed.'
  print '#'
  print '# ------------------------------ VERSION 2.180814 ------------------------------'
  print '#'
  print '# [CHA] Modified FLATTEN_IMAGES reducing number of iterations at largest scale'
  print '#       from 3 to 2 (when median-filtering standard deviations images).'
  print '#'
  print '# ------------------------------ VERSION 2.180727 ------------------------------'
  print '#'
  print '# [CHA] Modified FLATTEN_IMAGES by removing the initial median filtering and'
  print '#       setting the first scale to beam size (previously, it was two beams). The'
  print '#       first scale is now not minimized, the minimization starts with the 2nd'
  print '#       scale and is done with all larger scales. Few other cosmetic changes.'
  print '#'
  print '# [CHA] Modified OPERATE: division of images is now done only in non-zero pixels'
  print '#       of the denominator image. The resulting divided image is not modified'
  print '#       (as before) by setting to NANs zero pixels of the denominator.'
  print '#'
  print '# ------------------------------ VERSION 2.180726 ------------------------------'
  print '#'
  print '# [FIX] Modified FLATTEN_IMAGES by setting the initial median filtering window'
  print '#       radius equal to 2 and 4 beam sizes, as they were before v2.180724.'
  print '#'
  print '# ------------------------------ VERSION 2.180724 ------------------------------'
  print '#'
  print '# [CHA] Modified FLATTEN_IMAGES by setting the initial median filtering window'
  print '#       radius equal to the beam size. Few other cosmetic changes.'
  print '#'
  print '# ------------------------------ VERSION 2.180622 ------------------------------'
  print '#'
  print '# [FIX] Corrected broken identification of the HIGHLIGHT utility in the system.'
  print '#       Removed unnecessary multiplications of images by observational mask. Now'
  print '#       the mask image should not necessarily be cut by 8 pixels on borders, it'
  print '#       may be unity in all pixels.'
  print '#'
  print '# ------------------------------ VERSION 2.180531 ------------------------------'
  print '#'
  print '# [CHA] Tuned the "stopitnow!" mechanism of prematurely finishing iterations.'
  print '#'
  print '# ------------------------------ VERSION 2.180409 ------------------------------'
  print '#'
  print '# [FIX] Minor fix in GETIMAGES: removal of z0.fits in concurrently running jobs'
  print '#       led to a crash; the zero image is now preserved.'
  print '#'
  print '# ------------------------------ VERSION 2.180210 ------------------------------'
  print '#'
  print '# [FIX] Minor fix and clarification in GETIMAGES.'
  print '#'
  print '# ------------------------------ VERSION 2.180209 ------------------------------'
  print '#'
  print '# [CHA] Minor usability improvements in function FLATTEN_IMAGES in GETIMAGES.'
  print '#'
  print '# ------------------------------ VERSION 2.180208 ------------------------------'
  print '#'
  print '# [CHA] Minor improvements in final median-filtered and flattened images.'
  print '#'
  print '# ------------------------------ VERSION 2.180114 ------------------------------'
  print '#'
  print '# [CHA] Modified FLATTEN_IMAGES removing change of minimum on last iteration,'
  print '#       in order to reduce negative areas in background-subtracted images.'
  print '#'
  print '# ------------------------------ VERSION 2.180112 ------------------------------'
  print '#'
  print '# [CHA] Minor improvements in GETIMAGES.'
  print '#'
  print '# ------------------------------ VERSION 2.180106 ------------------------------'
  print '#'
  print '# [FIX] Minor improvements in GETIMAGES for copying, moving, and deleting files.'
  print '#'
  print '# ------------------------------ VERSION 2.180105 ------------------------------'
  print '#'
  print '# [FIX] Minor fixes in GETIMAGES. Now proper error handling (SIGINT, SIGTERM).'
  print '#'
  print '# ------------------------------ VERSION 2.180103 ------------------------------'
  print '#'
  print '# [CHA] Many minor changes (mainly highlighting) in GETIMAGES.'
  print '#'
  print '# ------------------------------ VERSION 2.171226 ------------------------------'
  print '#'
  print '# [CHA] Minor modifications in GETSOURCES. Also cosmetic changes related to the'
  print '#       colorized output.'
  print '#'
  print '# ------------------------------ VERSION 2.171220 ------------------------------'
  print '#'
  print '# [CHA] Minor improvements in GETIMAGES to make the screen output colorized.'
  print '#       The PRINTF command is used with functions PRINT, PRINTL, and PRINTLH,'
  print '#       as well as the utility HIGHLIGHT, if the latter is available.'
  print '#'
  print '# ------------------------------ VERSION 2.171213 ------------------------------'
  print '#'
  print '# [CHA] Generalized median filtering in IMGSTAT to enable using both exact and'
  print '#       a faster approximate (binning) method for median calculation.'
  print '#'
  print '# ------------------------------ VERSION 2.171212 ------------------------------'
  print '#'
  print '# [CHA] Generalized FLATTEN_IMAGES in GETIMAGES to enable restarting from any'
  print '#       flattening iteration.'
  print '#'
  print '# ------------------------------ VERSION 2.171208 ------------------------------'
  print '#'
  print '# [CHA] Simplified GETSOURCES by removing the median-filtering and flattening,'
  print '#       which now is available only separately in GETIMAGES.'
  print '#'
  print '# ------------------------------ VERSION 2.171204 ------------------------------'
  print '#'
  print '# [CHA] Accelerated median filtering in IMGSTAT by using an approximate binning'
  print '#       method (function MEDAPPROX in TOOLS). The default number of bins (1000)'
  print '#       is optimal; the acceleration factor relative to SELECTK is ~1.5.'
  print '#'
  print '# ------------------------------ VERSION 2.171128 ------------------------------'
  print '#'
  print '# [CHA] Slightly improved screen output of FLATTEN_IMAGES.'
  print '#'
  print '# ------------------------------ VERSION 2.170428 ------------------------------'
  print '#'
  print '# [NEW] Added a possibility to restart iterations in FLATTEN_IMAGES.'
  print '#'
  print '# ------------------------------ VERSION 2.170420 ------------------------------'
  print '#'
  print '# [NEW] Introduced iterations in FLATTEN_IMAGES.'
  print '#'
  print '# ------------------------------ VERSION 2.170404 ------------------------------'
  print '#'
  print '# [CHA] Simplified GETIMAGES removing old stuff left over from GETSOURCES.'
  print '#'
  print '# ------------------------------ VERSION 2.170315 ------------------------------'
  print '#'
  print '# [NEW] FLATTEN_IMAGES from GETSOURCES script isolated into GETIMAGES script.'
  print '#    => IMPORTANT. When using GETIMAGES with GETSOURCES v1.140127, an updated'
  print '#       utility IMGSTAT from GETSOURCES v2 must replace the original IMGSTAT.'
  print '#'
  print '# ------------------------------ VERSION 2.170129 ------------------------------'
  print '#'
  print '# [CHA] Few minor improvements in FLATTEN_IMAGES.'
  print '#'
  print '# ------------------------------ VERSION 2.170128 ------------------------------'
  print '#'
  print '# [CHA] Recoded the function SELECTK (in IMAGESTAT) to possibly accelerate'
  print '#       median filtering.'
  print '#'
  print '# ------------------------------ VERSION 2.170127 ------------------------------'
  print '#'
  print '# [CHA] Median filtering is again in one iteration with two median filtering'
  print '#       iterations for the largest sliding window.'
  print '#'
  print '# ------------------------------ VERSION 2.170125 ------------------------------'
  print '#'
  print '# [CHA] Median filtering is again in one iteration with extra median filtering'
  print '#       iteration(s) for the last (largest) sliding window.'
  print '#'
  print '# ------------------------------ VERSION 2.170124 ------------------------------'
  print '#'
  print '# [CHA] Modified in FLATTEN_IMAGES the image to be median filtered. Now it is'
  print '#       not the original image, but already somewhat median filtered image.'
  print '#       Median filtering is now done in two iterations.'
  print '#'
  print '# ------------------------------ VERSION 2.170113 ------------------------------'
  print '#'
  print '# [CHA] Removed from FLATTEN_IMAGES an additional median filtering for stdev'
  print '#       images with a larger window compared to the usual images.'
  print '#'
  print '# ------------------------------ VERSION 2.170112 ------------------------------'
  print '#'
  print '# [CHA] Removed all additional mean filtering from FLATTEN_IMAGES and increased'
  print '#       maximum sliding window sizes from 3 x SRCMAXSIZE to 4 x SRCMAXSIZE.'
  print '#'
  print '# ------------------------------ VERSION 2.170111 ------------------------------'
  print '#'
  print '# [CHA] Minor improvements in FLATTEN_IMAGES. Changed *initial* filtering from'
  print '#       mean to median, as mean filtering may give some ring or edge artifacts.'
  print '#'
  print '# ------------------------------ VERSION 2.170110 ------------------------------'
  print '#'
  print '# [CHA] Added in FLATTEN_IMAGES mean filtering after median filtering to smooth'
  print '#       possible sharp features in derived background; few other improvements.'
  print '#'
  print '# ------------------------------ VERSION 2.170109 ------------------------------'
  print '#'
  print '# [CHA] Replaced in FLATTEN_IMAGES additional convolution with additional mean'
  print '#       filtering when deriving flattening image.'
  print '#'
  print '# ------------------------------ VERSION 2.170107 ------------------------------'
  print '#'
  print '# [CHA] Modified FLATTEN_IMAGES by computing flattening image from the stdev'
  print '#       image of background-subtracted image rather than of the original image.'
  print '#'
  print '# ------------------------------ VERSION 2.170106 ------------------------------'
  print '#'
  print '# [CHA] Modified FLATTEN_IMAGES in the definition of the first window radius for'
  print '#       median filtering of the original image.'
  print '#'
  print '# ------------------------------ VERSION 2.170104 ------------------------------'
  print '#'
  print '# [FIX] Generalized IMGSTAT by appending to the output name ".imgstat{1|2}" the'
  print '#       number of the wavelength being processed, to enable running the code in'
  print '#       the same current directory concurrently in different terminals.'
  print '#'
  print '# ------------------------------ VERSION 2.170103 ------------------------------'
  print '#'
  print '# [CHA] Generalized the naming of log files. Now the code can be run in the'
  print '#       same current directory concurrently in different terminals. If the waves'
  print '#       are processe in parallel, the log files are appended with the number of'
  print '#       the wavelength being processed in the run.'
  print '#'
  print '# ------------------------------ VERSION 2.170102 ------------------------------'
  print '#'
  print '# [FIX] Corrected creation of the observed filaments and filament-subtracted'
  print '#       measurement images for consistency with new image flattening algorithm.'
  print '#'
  print '# ------------------------------ VERSION 2.170101 ------------------------------'
  print '#'
  print '# [CHA] Accelerated FLATTEN_IMAGES by using a larger scale factor (~1.3) and'
  print '#       thus reducing the number of time-consuming median filterings. Images of'
  print '#       standard deviations are processed in windows larger than R = 3 x BEAM,'
  print '#       as small-scale minima in stdev images would lead to spurious peaks in'
  print '#       flattened images.'
  print '#'
  print '# ------------------------------ VERSION 2.161229 ------------------------------'
  print '#'
  print '# [CHA] Minor refinements of the image flattening procedure.'
  print '#'
  print '# ------------------------------ VERSION 2.161227 ------------------------------'
  print '#'
  print '# [CHA] Removed all items related to the old flattening algorithm and to the'
  print '#       second (final) extraction; now a single extraction is sufficient.'
  print '#'
  print '# ------------------------------ VERSION 2.161224 ------------------------------'
  print '#'
  print '# [CHA] Improved FLATTEN_IMAGES by determining an stdev image directly from the'
  print '#       original observed image. Both background and stdev images are processed'
  print '#       in exactly the same algorithm of median filtering.'
  print '#'
  print '# ------------------------------ VERSION 2.161222 ------------------------------'
  print '#'
  print '# [CHA] Modified FLATTEN_IMAGES by defining lastscale = 3 x SRCMAXSIZE.'
  print '#'
  print '# ------------------------------ VERSION 2.161220 ------------------------------'
  print '#'
  print '# [CHA] Simplified FLATTEN_IMAGES by removing unnecessary mask & image expansion'
  print '#'
  print '# ------------------------------ VERSION 2.161219 ------------------------------'
  print '#'
  print '# [CHA] Modified FLATTEN_IMAGES, restoring initial median filtering back to a'
  print '#       window radius of 1/3 x BEAM. This should give a better background.'
  print '#'
  print '# ------------------------------ VERSION 2.161217 ------------------------------'
  print '#'
  print '# [CHA] Further slight modification to improve flattening.'
  print '#'
  print '# ------------------------------ VERSION 2.161214 ------------------------------'
  print '#'
  print '# [CHA] Further slight modification of FLATTEN_IMAGES.'
  print '#'
  print '# ------------------------------ VERSION 2.161213 ------------------------------'
  print '#'
  print '# [CHA] Modified FLATTEN_IMAGES by shifting initial median filtering from window'
  print '#       radius of 1/3 x BEAM to BEAM.'
  print '#'
  print '# ------------------------------ VERSION 2.161212 ------------------------------'
  print '#'
  print '# [CHA] Modified FLATTEN_IMAGES: eliminated parameter MAXOBSBGSCALE and defined'
  print '#       lastscale = [(3 x BEAM)^4 + (SRCMAXSIZE)^4]^(1/4).'
  print '#'
  print '# ------------------------------ VERSION 2.161207 ------------------------------'
  print '#'
  print '# [CHA] Introduced small-scale noise flattening in FLATTEN_IMAGES. With this,'
  print '#       detection image filament subtraction from measurement images is removed.'
  print '#'
  print '# ------------------------------ VERSION 2.161205 ------------------------------'
  print '#'
  print '# [CHA] Added in FLATTEN_IMAGES another step: small-scale noise flattening.'
  print '#'
  print '# ------------------------------ VERSION 2.161128 ------------------------------'
  print '#'
  print '# [CHA] Default values of mfitermax and NSIGMACUTFLOOR were both set to 6.'
  print '#'
  print '# [CHA] Further changed FLATTEN_IMAGES. lastscale is increased to 3 x BEAM,'
  print '#       to truncate Gaussian sources approximately at its half-maximum. Also'
  print '#       MAXOBSBGSCALE is defined as max { 1/3 x SRCMAXSIZE, lastscale } and'
  print '#       EFFRES is decreased from 1.1 to 1.05 x BEAM.'
  print '#'
  print '# ------------------------------ VERSION 2.161123 ------------------------------'
  print '#'
  print '# [CHA] Modified behavior of FLATTEN_IMAGES depending on MAXOBSBGSCALE. If the'
  print '#       latter is set to 0, background image for that scale is not computed.'
  print '#       The original observed image (not background-subtracted one) should then'
  print '#       be used for measurements.'
  print '#'
  print '# ------------------------------ VERSION 2.161121 ------------------------------'
  print '#'
  print '# [CHA] Modified FlATTEN_IMAGES, changing the meaning of MAXOBSBGSCALE. It is'
  print '#       now a FWHM of an equivalent Gaussian, therefore MAXOBSBGSCALE is set to'
  print '#       6 x SRCMAXSIZE (which is a FWHM) at a wavelength. In other words, it is'
  print '#       3 x footprintl diameter.'
  print '#'
  print '# ------------------------------ VERSION 2.161118 ------------------------------'
  print '#'
  print '# [CHA] Further fine-tuning of background estimation using median filtering:'
  print '#       lastscale is increased from BEAM to twice the BEAM.'
  print '#'
  print '# ------------------------------ VERSION 2.161116 ------------------------------'
  print '#'
  print '# [CHA] Further improvements in preparation of images using median filtering.'
  print '#'
  print '# [FIX] Removed very minor fix of the accumulated smallest scale to be cleaned,'
  print '#       introduced in v.2.161014. The smallest scale is quite noisy and better'
  print '#       not to clean and not to use it for detection. The next scale is better.'
  print '#'
  print '# ------------------------------ VERSION 2.161111 ------------------------------'
  print '#'
  print '# [CHA] Further improvements in background estimation using median filtering.'
  print '#'
  print '# ------------------------------ VERSION 2.161101 ------------------------------'
  print '#'
  print '# [CHA] Default values of mfitermax and NSIGMACUTFLOOR were set to 7 and 6.'
  print '#'
  print '# ------------------------------ VERSION 2.161027 ------------------------------'
  print '#'
  print '# [NEW] Added in IMGSTAT an option of mean filtering in a sliding window.'
  print '#'
  print '# ------------------------------ VERSION 2.161026 ------------------------------'
  print '#'
  print '# [FIX] Minor fix in subtraction of the full reconstructed filaments from the'
  print '#       original measurement and detection images.'
  print '#'
  print '# ------------------------------ VERSION 2.161025 ------------------------------'
  print '#'
  print '# [CHA] Default values of mfitermax and NSIGMACUTFLOOR increased to 7 and 5.'
  print '#'
  print '# ------------------------------ VERSION 2.161024 ------------------------------'
  print '#'
  print '# [CHA] Minor improvements in background estimation using median filtering.'
  print '#'
  print '# ------------------------------ VERSION 2.161022 ------------------------------'
  print '#'
  print '# [NEW] Added background estimation using median filtering.'
  print '#_______________________________________________________________________________'
  print '#'
}
#___________________________________________________________________________________________________________________________________

CONDITIONS ()
{
  print '#  _____________________________________________________________________________'
  print '# |                                                                             '
  print '# |                                  GETIMAGES                                  '
  print '# |               Background Derivation & Image Flattening Method               '
  print '# |             Alexander Men’shchikov, DAp IRFU CEA Saclay, France             '
  print '# |                               '"$versnum"'                               '
  print '# |_____________________________________________________________________________'
  print '# |                                                                             '
  print '# |    CONDITIONS OF USE (*1)                                                   '
  print '# |                                                                             '
  print '# | In what follows, "this software" refers to the Bash scripts GETIMAGES,      '
  print '# | GETSOURCES, CONVOLVE, HIRESCOLDENS, IOSPEED, PREPAREOBS, RESAMPLE, and the  '
  print '# | associated FORTRAN utilities CLEANBG, ELLIPSES, EXPANDA, FFTCONV, FINALCAT, '
  print '# | FITFLUXES, FMEASURE, IMGSTAT, MODFITS, OPERATE, READHEAD, SEPARATE, SFINDER,'
  print '# | SMEASURE, SPLITCUBE, and the library TOOLS, and "the author" refers to      '
  print '# | Alexander Men’shchikov of the Département d’Astrophysique, IRFU, CEA Saclay,'
  print '# | France.                                                                     '
  print '# |                                                                             '
  print '# | It is assumed that anyone using this code ("the user") has read, understood,'
  print '# | and agreed to the following conditions of use:                              '
  print '# |                                                                             '
  print '# | 1. Distribution of this software shall remain the purview of the author. A  '
  print '# |    user is free to share this code with co-workers and students, but if it  '
  print '# |    is requested by a colleague not working directly with the user, the user '
  print '# |    is asked to redirect such requests to: alexander.menshchikov@cea.fr      '
  print '# |                                                                             '
  print '# | 2. This software shall be used exclusively for education, research, non-    '
  print '# |    profit, and non-military purposes. Specific written permission from the  '
  print '# |    author must be obtained before any commercial use of this software is    '
  print '# |    undertaken.                                                              '
  print '# |                                                                             '
  print '# | 3. The banners of this software, these conditions of use, and information   '
  print '# |    about the author shall remain with this software and any descendent      '
  print '# |    developed from and still based substantially upon this software.         '
  print '# |                                                                             '
  print '# | 4. The names of the institutions with which the author is or has been       '
  print '# |    affiliated shall not be used to publicise any data and (or) results      '
  print '# |    generated by this software. All findings and their interpretation are    ' 
  print '# |    the opinions of the user and do not necessarily reflect those of the     '
  print '# |    author nor the institutions with which the author is or has been         '
  print '# |    affiliated.                                                              '
  print '# |                                                                             '
  print '# | The author makes no representations about the suitability of this software  '
  print '# | for any purpose. Subject to the above conditions, this software is provided '
  print '# | "as is" without any expressed or implied warranty.                          '
  print '# |                                                                             '
  print '# | Inquiries about the code, bug reports, constructive criticism, etc., can be '
  print '# | directed to: alexander.menshchikov@cea.fr                                   '
  print '# |                                                                             '
  print '# |    CITATIONS AND ACKNOWLEDGEMENTS (*1)                                      '
  print '# |                                                                             '
  print '# | The algorithms used in this software are described in:                      '
  print '# | Men’shchikov A. 2017, A&A, 607, A64                                         '
  print '# | Men’shchikov A. 2013, A&A, 560, A63                                         '
  print '# | Men’shchikov A., Andre Ph., Didelon P., et al. 2012, A&A, 542, A81          '
  print '# | Men’shchikov A., Andre Ph., Didelon P., et al. 2010, A&A, 518, L103         '
  print '# |                                                                             '
  print '# | It is requested that any publication reporting results obtained using this  '
  print '# | software or any of its derivatives includes something like this:            '
  print '# |                                                                             '
  print '# | "Use of GETSOURCES (GETFILAMENTS), developed by A. Men’shchikov at the DAp, '
  print '# | IRFU, CEA Saclay, France, is hereby acknowledged." (*2)                     '
  print '# |_____________________________________________________________________________'
  print '#                                                                               '
  print '# (*1) Adapted from those written by David Clarke for AZEuS, an MHD code.       '
  print '# (*2) If only individual scripts or utilities (not the GETSOURCES script) were '
  print '#      used, substitute "GETSOURCES (GETFILAMENTS)" with their names.'
}
#___________________________________________________________________________________________________________________________________

USAGE ()
{
  printlh
  printlh '                         GETIMAGES: USAGE INFORMATION                         '
  printlh
  printlh ' '$scriptname' [<revision>|stopitnow|help]'
  printlh
  printlh ' Input data must be defined in a configuration file '\''+getimages.cfg'\''.' 
  printlh ' If no config file exists in the current directory, a new default version of '
  printlh ' the file is saved in that directory. One needs to edit the file as necessary,'
  printlh ' then run '$scriptname' again.'
  printlh
  printlh ' An optional command-line parameter "<revision>" can be given if the user wants'
  printlh ' to run a specific version (revision) of GETIMAGES. The parameter must be a'
  printlh ' six-digit number representing the revision in the format yymmdd (year-month-'
  printlh ' day). For example, if one wants to use the version 1.110320 of GETSOURCES and'
  printlh ' its utilities and scripts, then <revision> is 110320. Then the version from'
  printlh ' directory $GETSOURCES_BIN/v1.110320 will be used; here $GETSOURCES_BIN is the'
  printlh ' environment variable pointing to a directory where the binaries are installed.'
  printlh   
  printlh ' Another optional parameter "stopitnow" can be used to gently stop GETIMAGES.'
  printlh ' The command "getimages stopitnow" issued in a working directory will exit'
  printlh ' iterations in GETIMAGES as soon as possible (note that it may still take a'
  printlh ' long time). Another way of getting the same effect is to create a file named'
  printlh ' "+stopitnow!" in the extraction directory.'
  printlh
  printlh ' GETIMAGES calls a suite of FORTRAN utilities (written by A. Men’shchikov).'
  printlh ' They use the CFITSIO library (by William D Pence).'
  printlh
  printlh '    NOTE 1: PROCESSING EFFICIENCY GUIDELINES.'
  printlh
  printlh ' 1. Always run GETIMAGES on a *local* (internal) hard drive physically'
  printlh '    attached to the computer (CPU) used for extractions. The code performs'
  printlh '    *lots* of read/write operations with FITS files containing images.'
  printlh '    A user would benefit from the fastest possible I/O throughput.'
  printlh
  printlh ' 2. Never use hard drives attached over the local network, as the disk access'
  printlh '    over the network is *very* slow compared to the local disks. GETIMAGES'
  printlh '    may be *extremely* slow on the disks attached over networks.'
  printlh
  printlh ' 3. Processing speed depends on many circumstances in a given computer system.'
  printlh '    If one wants to optimize the processing times, it is generally a good idea'
  printlh '    to test available disk storage for speed.'
  printlh
  printlh ' 4. The script IOSPEED enables comparisons of available hard drives. One can'
  printlh '    perform tests of the local and network disks, as well as of the virtual'
  printlh '    RAM disks.'
  printlh
  printlh ' For the information on how to install GETSOURCES and to quickly start using it'
  printlh ' refer to the files "INSTALLATION.GUIDE", "QUICK.START.GUIDE", "USERS.CHECKLIST"'
  printlh ' included in the distribution archive.'
  printlh
  printlh ' The algorithms used in this software are described in:'
  printlh ' Men’shchikov A. 2017, A&A, 607, A64'
  printlh ' Men’shchikov A. 2013, A&A, 560, A63'
  printlh ' Men’shchikov A., Andre Ph., Didelon P., et al. 2012, A&A, 542, A81'
  printlh ' Men’shchikov A., Andre Ph., Didelon P., et al. 2010, A&A, 518, L103'
}
#___________________________________________________________________________________________________________________________________
#
system=`uname`; homepath=`printf "%s\n" ~`; originpath=`pwd`; executing="$0"; scriptname=`basename "$executing"`; hell='/dev/null'
script=`printf "%s\n" $scriptname | tr a-z A-Z`; configfile='+getimages.cfg'; nwmax=0; nwaves=0; nwave=0; cwavenum=''
export top_pid=$$; trap 'killhandler' SIGINT; trap 'errorhandler' SIGTERM; trapfirst='yes'; outsameline='no'
#
killhandler ()
{
  if [[ "$waitingforinput" == 'yes' ]]
  then 
    if [[ "$trapfirst" == "yes" ]]; then printlh; fi
    printlh ' HINT: To properly kill '$script': (1) hit ENTER, (2) press CTRL-C '
    trapfirst='no'
  else 
    if [[ "$outsameline" == 'yes' ]]; then printlh; fi
    printlh; printlh ' '$script': ABORTED BY USER REQUEST (CTRL-C)'; printlh; exit
  fi
}
errorhandler ()
{
  printlh; printlh ' '$script': ERRORS ENCOUNTERED: ABORTED'; printlh; exit
}
#___________________________________________________________________________________________________________________________________
#
# Default values of the Herschel beams adopted by GETIMAGES (averaged). 
# When modifying them here, don't forget to change them also in GETSOURCES and PREPAREOBS!
#
b070s10=' 5.4'; b070s20=' 5.6'; b070s20p=' 5.9'; b070s60=' 7.2'; b070s60p=' 8.4'
b100s10=' 6.7'; b100s20=' 6.8'; b100s20p=' 7.0'; b100s60=' 8.2'; b100s60p=' 9.4'
b160s10='11.2'; b160s20='11.3'; b160s20p='11.7'; b160s60='12.3'; b160s60p='13.5'; b250='18.2'; b350='24.9'; b500='36.3'
#___________________________________________________________________________________________________________________________________

DEFAULTPARAMS ()  
{ # Default parameters will be saved in a configuration file if no such file exists in the current directory.
  
  for (( v=1; v <= 4; v++ ))
  do 
  print '#'"${headerinfo[v]:3}"; done
  print '#                    |                                        !'
  print '# VALUES             | VARIABLES'\'' NAMES    {CHOICES/DEFAULTS} ! COMMENTS [NOTE: Default beams: Herschel fast (60"/s) parallel mode]'
  print '#____________________|________________________________________!_____________________________________________________________________'
  print '#'
  print ' name~of~your~field  | prefix ....( should not contain dots ) ! name that will be used as a prefix for naming output files'
  print ' 6                   | nwmax .................{ nwmax < 100 } ! maximum number of wavelengths in this configuration file'
  print ' 6 1 2 3 4 5 6       | nwaves nw[i] ......................... ! number of waves and indices of the wavelengths to process'
  print ' 070 '$b070s60p' <x070>     | ...................................... ! wavelength, beamsize ("FWHM), maxsize ("FWHM)'
  print ' /prepared/images .  | ...................................... ! paths to prepared images and to processing directory'
  print ' prepared.image.070  | ...................................... ! name of the prepared image for background subtracton and flattening'
  print ' 100 '$b100s60p' <x100>     | ...................................... ! wavelength, beamsize ("FWHM), maxsize ("FWHM)'
  print ' /prepared/images .  | ...................................... ! paths to prepared images and to processing directory'
  print ' prepared.image.100  | ...................................... ! name of the prepared image for background subtracton and flattening'
  print ' 160 '$b160s60p' <x160>     | ...................................... ! wavelength, beamsize ("FWHM), maxsize ("FWHM)'
  print ' /prepared/images .  | ...................................... ! paths to prepared images and to processing directory'
  print ' prepared.image.160  | ...................................... ! name of the prepared image for background subtracton and flattening'
  print ' 250 '$b250' <x250>     | ...................................... ! wavelength, beamsize ("FWHM), maxsize ("FWHM)'
  print ' /prepared/images .  | ...................................... ! paths to prepared images and to processing directory'
  print ' prepared.image.250  | ...................................... ! name of the prepared image for background subtracton and flattening'
  print ' 350 '$b350' <x350>     | ...................................... ! wavelength, beamsize ("FWHM), maxsize ("FWHM)'
  print ' /prepared/images .  | ...................................... ! paths to prepared images and to processing directory'
  print ' prepared.image.350  | ...................................... ! name of the prepared image for background subtracton and flattening'
  print ' 500 '$b500' <x500>     | ...................................... ! wavelength, beamsize ("FWHM), maxsize ("FWHM)'
  print ' /prepared/images .  | ...................................... ! paths to prepared images and to processing directory'
  print ' prepared.image.500  | ...................................... ! name of the prepared image for background subtracton and flattening'
  print ' y 1.41421  ..user.. | flattening sfac .........{ y 1.41421 } ! {y|y-|n}{#} background-subtract & flatten images and scale factor'
  print ' 2          ..user.. | savespace ..........{ 0~2; optimal 2 } ! {0|1|2} disk space saving level: clean up more and more files'
  print ' 0          ..user.. | verbosity ..........{ 0~2; optimal 0 } ! {0|1|2} verbosity level for the screen and log file output'
  print '#___________________________________________________________________________________________________________________________________'
} 
#___________________________________________________________________________________________________________________________________

CLEAR_STDIN ()  
# Discard all stdin buffer before prompting user input.
{
  old_tty_settings=`stty -g`
  stty -icanon min 0 time 0
  while read none; do :; done 
  stty "$old_tty_settings"
}
#___________________________________________________________________________________________________________________________________

print () 
{ 
  (printf "%s\n" "$*")
}
printl () 
{ 
  (printf "$*") | 'tee' -a $logfile
}
printlh () 
{ 
  (printf "%s\n" "$*") | "tee" -a $logfile | $highlight
}
#___________________________________________________________________________________________________________________________________

lenexe=${#executing}; lenscr=${#scriptname}
lastps=$(( lenexe - lenscr )); lastp=$(( lastps - 1 ))
if [[ "${executing:lastps:lenexe}" == "$scriptname" ]]
then path=${executing:0:lastp}
else path=''
fi
themefile=$path'/sasha.theme'
langfile=$path'/sasha.lang'
highlight='highlight'
$highlight --version 2>$hell >$hell; rchigh=$?
if [[ $rchigh -gt 0 && -e $langfile && -e $themefile ]]
then 
  highlight='/opt/local/bin/highlight'
  $highlight --version 2>$hell >$hell; rchigh=$?
fi
if [[ $rchigh -eq 0 ]]
then 
  highlight=$highlight" --out-format=xterm256 --config-file=$themefile --config-file=$langfile"
else
  highlight='tee'
fi
#___________________________________________________________________________________________________________________________________

tehi ()
{
  if [[ $verb != '-verb0' ]]
  then 'tee' -a $logfile | $highlight
  else 'tee'
  fi
  rc=`cat < .rc | { IFS= read -r first; printf "%s\n" "$first"; }`; chk_rc $rc
}
#___________________________________________________________________________________________________________________________________

chk_rc ()
{ 
  rc=$?
  if [[ $rc -eq 0 && "$1" != "" ]]; then rc=$1; fi
  if [[ $rc -gt 0 ]]
  then 
    where=$script
    if [[ $fun != "" ]]
    then where="$where: '$fun'"; fi
    printlh " $where: ERROR (RC=$rc)"
    kill -SIGTERM $top_pid
  fi
  return $rc
}
#___________________________________________________________________________________________________________________________________

ERROR ()
{
  local message=$1; local advice1=$2; local advice2=$3
  printlh
  printlh ' '$script': '\'$fun\'': ERROR:'
  printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
  printlh ' INFO: '$message
  printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
  if [[ "$advice1" != '' ]]; then printlh ' HINT: '$advice1
  if [[ "$advice2" != '' ]]; then printlh '       '$advice2; fi
  printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
  fi
  printlh ' ERROR in '$script': Aborted.'
  printlh
  exit 99
}
#___________________________________________________________________________________________________________________________________

LEADZEROS ()
{
  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 ()
{
  local funo=$fun; fun='POSITION'
  if [[ $2 == "" ]]
  then 
    'ERROR' 'Too few parameters: "'$1'" "'$2'"'
  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
}
#___________________________________________________________________________________________________________________________________

if [[ -e $configfile ]]
then
{
while read -r configline
do
  'POSITION' 'nwmax' "$configline"; ncv=$returnposition
  if [[ $ncv -gt 0 ]]
  then nwmax=${configline:0:2}; nwmax=${nwmax/\ /}
  fi
  'POSITION' 'nwaves' "$configline"; ncv=$returnposition
  if [[ $ncv -gt 0 ]]
  then nwaves=${configline:0:2}; nwave=${configline:2:2}; nwaves=${nwaves/\ /}; nwave=${nwave/\ /}
  fi
  if [[ $nwaves -gt 0 && $nwaves -lt $nwmax ]]
  then 'LEADZEROS' 1 $nwave znw; cwavenum='.'$znw; break
  fi
done
} < $configfile
fi
logfile='+log.'$script$cwavenum
logit='+log.'$script'.iter'$cwavenum
#___________________________________________________________________________________________________________________________________

CHECK_IMGSTAT ()
{
  local minimum=$1; local maximum=$2; local sigma=$3; local scriptname=$4; local funname=$5; local funoname=$6; local place=$7

  if [[ ${minimum:0:1} == "*" || ${maximum:0:1} == "*" || ${sigma:0:1} == "*" || $sigma == "0.0000000000000000000" ]]
  then 
    printlh
    printlh ' '$scriptname': '$funname': ERROR:'
    printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    printlh ' WARNING: Results from IMGSTAT ('$place'): min max std: '$minimum $maximum $sigma
    printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    printlh; fun=$funoname
    return 99
  else 
    \rm -f '.+imgstat2'$cwavenum 2>$hell
  fi
}
#___________________________________________________________________________________________________________________________________

COMPARE_NUMBERS ()
{
  local funo=$fun; fun='COMPARE_NUMBERS'
  if [[ $3 == "" ]]
  then 
    'ERROR' 'Too few parameters: "'$1'" "'$2'"'
  fi
  if [[ $2 != ">" && $2 != "<" && $2 != "=" && $2 != "!=" && $2 != ">=" && $2 != "<=" ]]
  then
    'ERROR' 'Unsupported operator: "'$2'"'
  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
}
#___________________________________________________________________________________________________________________________________

PROGRESS_BAR ()
{ 
  local funo=$fun; fun='PROGRESS_BAR'
  if [[ $1 != "complete" && $6 == "" ]]
  then 
    'ERROR' 'Too few parameters: "'$1'" "'$2'" "'$3'" "'$4'" "'$5'" "'$6'"'
  fi
  local i; local nbeg; local nend; local width; local ncur; local sym='•'; local nrange; local iter=0
  progresstitle=$1; ncur=$2; nbeg=$3; nend=$4; width=$5; char=$6; iter=$7

  if [[ "$1" != "complete" ]]
  then
    if [[ $ncur -eq $nbeg ]]
    then                
      lniter=${#iter}; lnend=${#nend}; ltitle=${#progresstitle}
      lbar=$(( width - 3 - ltitle - 1 - lnend - 1 ))
      if [[ $lbar -gt $nend ]]
      then lbar=$nend
      fi
      lbarp1=$(( lbar + 1 + lniter + 1 ))
      lback=$(( lbarp1 + lnend + 1 ))
      if [[ "$char" == "·" ]]
      then
        printl "%s\n"; printl "%s" "$progresstitle"' ['
      else 
        printf "%b" "] $nend $iter\033[${lback}D"
        printf "%b" "] $nend [" >> $logfile
      fi
      for (( i=1; i <= $lbar; i++ ))
      do printf "%b" "${char:0:1}"
      done
      printf "%b" "] $nend $iter\033[${lback}D"
      nprogress=1; npos=0; sleep 1
    else
      if [[ $nprogress -eq 0 ]]
      then
        'ERROR' 'Initial value /= lower bound: '$ncur $nbeg $nend $width
      fi
    fi
  fi
  nbegm1=$(( nbeg - 1 )); nrange=$(( nend - nbegm1 ))
  ncurm1=$(( ncur - 1 )); lncur=${#ncur}; lncurm1=${#ncurm1}; lncurm1p1=$(( lncurm1 + 1 ))
  lenr=$(print "scale=20; $nbegm1+($ncur-($nbeg-1))*($lbar-$nbeg)/$nrange" | bc); chk_rc
  if [[ ${lenr:0:1} != "-" ]]
  then leni=${lenr%.*}
  else leni=0
  fi

  if [[ $1 == "complete" ]]
  then 
    if [[ $ncur -ge $nend ]]
    then
      many=''
      for (( i=1; i <= $lncur; i++ ))
      do many=$many$sym
      done
      printf "%b" "\033[${lncur}D$many"
      printf "%b" "$sym" >> $logfile
      printf "%b" "] $nend $iter" >> $logfile
    else
      many=''
      ladd=$(( lbar - npos + 1))
      for (( i=1; i <= $ladd; i++ ))
      do many=$many"·"
      done
      printf "%b" "$many] $nend $iter" >> $logfile
    fi
    printlh
    fun=$funo
    return 0
  fi
  lenrold=$lenr
  leniold=$leni
  if [[ $leni -lt $nprogress ]]
  then
    if [[ $ncur -eq $nbeg ]]
    then 
      npos=$(( npos + lncur ))
      printf "%b" "$ncur"
      sleep 1
    else 
      npos=$(( npos - lncurm1 + lncur ))
      if [[ $npos -gt $lbar ]]
      then 
        lncorr=$(( npos - lbar ))
        npos=$(( npos - lncorr ))
        printf "%b" "\033[${lncorr}D"
      fi
      printf "%b" "\033[${lncurm1}D$ncur"
    fi
  else
    nprogress=$(( nprogress + 1 ))
    many=''
    for (( i=1; i <= $lncurm1p1; i++ ))
    do many=$many$sym
    done
    if [[ $npos -lt $lbar ]]
    then
      npos=$(( npos + 1 ))
      printf "%b" "$sym" | 'tee' -a $logfile
    fi
    printf "%b" "\033[${lncurm1p1}D$many\033[${lncur}D$ncur"
  fi
  
  fun=$funo
}
#___________________________________________________________________________________________________________________________________

GET_VERSION ()
{
  local funo='NO_FUN'; fun='GET_VERSION'
  cdate=`date +%d\ %b\ %Y\ %a\ %H:%M:%S\ %Z`
           
  if [[ "$1" == 'stopitnow' ]]
  then
    printlh ' ____________________________________________________________________________'
    printlh '                                                                             '
    printlh '  Telling '$script' to finish as soon as possible... you may have to wait... '
    printlh ' ____________________________________________________________________________'
    printlh
    print '#_______________________________________________________________________________'  > '+stopitnow!'
    print '#                                                                               ' >> '+stopitnow!'
    print '#              '$script' • Background Derivation & Image Flattening             ' >> '+stopitnow!'
    print '#              Alexander Men’shchikov, DAp IRFU CEA Saclay, France              ' >> '+stopitnow!'
    print '#_______________________________________________________________________________' >> '+stopitnow!'
    print '#                                                                               ' >> '+stopitnow!'
    print ' This file can be used to gently stop iterations in '$script                      >> '+stopitnow!'
    print '                                                                                ' >> '+stopitnow!'
    print ' When a user copies this file (or creates a file with the name "+stopitnow!") in' >> '+stopitnow!'
    print ' a directory where '$script' is running, the iterations in '$script' will be    ' >> '+stopitnow!'
    print ' exited as soon as possible.                                                    ' >> '+stopitnow!'
    print '                                                                                ' >> '+stopitnow!'
    print ' One gets the same effect by running '$script' with the command-line parameter  ' >> '+stopitnow!'
    print ' "stopitnow" in the directory where it is running, i.e. "'$script' stopitnow".  ' >> '+stopitnow!'
    exit 0
  fi
  
  if [[ "$1" == 'help' ]]
  then
    for (( b=1; b <= 7; b++ )); do printlh "${banner[b]}"; done
    printlh
    printlh '                         '$cdate
    
    'USAGE'
    
    printlh '_______________________________________________________________________________'
    printlh
    exit 0
  fi

  revision=$1

  if [[ $GETSOURCES_BIN == "" ]]
  then
    for (( b=1; b <= 7; b++ )); do printlh "${banner[b]}"; done
    printlh
    printlh '                         '$cdate
    
    'USAGE'
    
    printlh '_______________________________________________________________________________'
    printlh
    printlh ' Executing: '"$executing"
    printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    printlh ' ERROR: Environment variable GETSOURCES_BIN not set;'
    printlh '        it must point to the binary directory (bin/+GETSOURCES)'
    printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    printlh '  HINT: Make sure that your ~/.bashrc file contains a line like this:'
    printlh '        export GETSOURCES_BIN=/some/path/to/bin/+GETSOURCES'
    printlh '        and ask your computer sysadmin for help if necessary'
    printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    printlh
    return 99
  fi

  getsdefault="$GETSOURCES_BIN/getimages"                                              
  getsrevision="$GETSOURCES_BIN/v$revision/getimages"
                                                                                        
  if [[ $GETSOURCES_DEV != "" ]]
  then getsdefault="$GETSOURCES_DEV/getimages"
  fi
  if [[ $GETSOURCES_REL != "" ]]
  then 
    getsrevision="$GETSOURCES_REL/v$revision/getimages"
  fi
  
  if [[ "$revision" == "" ]]
  then 
    if [[ ! -e $getsdefault ]]
    then
      for (( b=1; b <= 7; b++ )); do printlh "${banner[b]}"; done
      printlh
      printlh '                         '$cdate

      'USAGE'

      printlh '_______________________________________________________________________________'
      printlh
      printlh ' Executing: '"$executing"

      'ERROR' 'Default installation of GETSOURCES problem: file not found: '$getsdefault
    fi
  else 
    if [[ -e $getsrevision ]]
    then
      $getsrevision; exit $?
    else
      for (( b=1; b <= 7; b++ )); do printlh "${banner[b]}"; done
      printlh
      printlh '                         '$cdate

      'USAGE'

      printlh '_______________________________________________________________________________'
      printlh
      printlh ' Executing: '"$executing"

      'ERROR' 'Revision '$revision' of GETSOURCES problem: file not found: '$getsrevision
    fi
  fi

  fun=$funo
}
#___________________________________________________________________________________________________________________________________

GET_PARAMETERS ()
{
  local funo=$fun; fun='GET_PARAMETERS'
  cdate=`date +%d\ %b\ %Y\ %a\ %H:%M:%S\ %Z`; local i; local j
  printlh
  printlh '                         '$cdate
  printlh
  printlh ' Executing: '"$executing"

  lenexe=${#executing}; lenscr=${#scriptname}
  lastps=$(( lenexe - lenscr )); lastp=$(( lastps - 1 ))
  
  if [[ "${executing:lastps:lenexe}" == "$scriptname" ]]
  then path=${executing:0:lastp}; else path=''; fi

  if [[ $path == "" ]]
  then 
    'ERROR' 'Trouble obtaining the path to '$script
  fi
  
  if [[ $path == '/Users/amenshch/Astronomy/+Computing/+GETSOURCES' ]]
  then
    path='/Users/amenshch/bin/+GETSOURCES'
  fi

# 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 8 ]]
    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
  } < $executing

  headerinfo[1]='#  ___________________________________________________________________________________________________________________________________'
  headerinfo[2]='# |                                                                                                                                   '
  headerinfo[3]='# | '$script' • Background Derivation & Image Flattening • '$versnum' • Alexander Men’shchikov, DAp IRFU CEA Saclay      '
  headerinfo[4]='# |___________________________________________________________________________________________________________________________________'

  print "${headerinfo[1]}"  > '+log.VERSION'
  print "${headerinfo[2]}" >> '+log.VERSION'
  print "${headerinfo[3]}" >> '+log.VERSION'
  print "${headerinfo[4]}" >> '+log.VERSION'

  'HISTORY'    > '+log.HISTORY'
  'CONDITIONS' > '+log.CONDITIONS'

  if [[ ! -e "$configfile" ]]
  then 
    'USAGE'; printlh '_______________________________________________________________________________'
    'DEFAULTPARAMS' > $configfile

    'ERROR' 'Configuration file '$configfile' does not exist' \
            'Default parameters saved in a new configuration file '$configfile
  fi
         
  if [[ "$1" != "" ]]; then printlh; printlh ' Parameters:' "$*"; fi

  \rm -f '+stopitnow!' 2>$hell

# Make sure that all my required utilities can be found in the system.
         
  printlh
  printlh ' '$script' WILL USE THE FOLLOWING FORTRAN UTILITIES:'
  {              
  $path'/expanda'  : 2>$hell; rc=$?                                                                                        
  if [[ $rc -eq 126 ]]; then 'ERROR' "Trouble executing '$path/expanda'..."; fi
  if [[ $rc -eq 127 ]]; then 'ERROR' "Utility '$path/expanda' not found or error loading shared libraries..."; fi
  if [[ $rc -eq 2 ]]; then expandaversion='2'; else expandaversion='1'; fi
  $path'/fftconv'  : 2>$hell; rc=$?
  if [[ $rc -eq 126 ]]; then 'ERROR' "Trouble executing '$path/fftconv'..."; fi
  if [[ $rc -eq 127 ]]; then 'ERROR' "Utility '$path/fftconv' not found or error loading shared libraries..."; fi
  if [[ $rc -eq 2 ]]; then fftconvversion='2'; else fftconvversion='1'; fi
  $path'/imgstat'  : 2>$hell; rc=$?
  if [[ $rc -eq 126 ]]; then 'ERROR' "Trouble executing '$path/imgstat'..."; fi
  if [[ $rc -eq 127 ]]; then 'ERROR' "Utility '$path/imgstat' not found or error loading shared libraries..."; fi
  if [[ $rc -eq 2 ]]; then imgstatversion='2'; else imgstatversion='1'; fi
  $path'/modfits'  : 2>$hell; rc=$?                                                                                        
  if [[ $rc -eq 126 ]]; then 'ERROR' "Trouble executing '$path/modfits'..."; fi
  if [[ $rc -eq 127 ]]; then 'ERROR' "Utility '$path/modfits' not found or error loading shared libraries..."; fi
  if [[ $rc -eq 2 ]]; then modfitsversion='2'; else modfitsversion='1'; fi
  $path'/operate'  : 2>$hell; rc=$?                                                                                        
  if [[ $rc -eq 126 ]]; then 'ERROR' "Trouble executing '$path/operate'..."; fi
  if [[ $rc -eq 127 ]]; then 'ERROR' "Utility '$path/operate' not found or error loading shared libraries..."; fi
  if [[ $rc -eq 2 ]]; then operateversion='2'; else operateversion='1'; fi
  $path'/readhead' : 2>$hell; rc=$?                                                                                        
  if [[ $rc -eq 126 ]]; then 'ERROR' "Trouble executing '$path/readhead'..."; fi
  if [[ $rc -eq 127 ]]; then 'ERROR' "Utility '$path/readhead' not found or error loading shared libraries..."; fi
  if [[ $rc -eq 2 ]]; then readheadversion='2'; else readheadversion='1'; fi
  if [[ $imgstatversion -eq 1 ]]
  then
    printlh; printlh
    printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    printlh " Incompatible version of '$path/imgstat'..."
    printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    printlh " FIX: Replace the old source file 'imgstat.for' file from GETSOURCES v1"
    printlh "      with new 'imgstat.for' and 'tools.for' from GETSOURCES v2, compile it,"
    printlh "      and run again. Send an email to the author for assistance."
    printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
    kill -SIGTERM $top_pid
  fi
  } | 'tee' -a $logfile | $highlight
  if [[ "| $highlight" == '| tee' ]]
  then printl "\n"
  else printf "\n" >> $logfile
  fi
  if [[ $rchigh -eq 0 && -e $langfile && -e $themefile ]]
  then
    $highlight --version 2>$hell >$hell; rc=$?; coloring='yes'
    if [[ $rc -eq 126 ]]; then printlh ' '$script': '$fun': ERROR:'; printlh; fi
    if [[ $rc -eq 126 ]]; then printlh " Trouble executing 'highlight'..."; coloring='no'; fi
    if [[ $rc -gt 0   ]]; then coloring='no'; fi
    if [[ $coloring == 'yes' ]]
    then
      lineh=`$highlight --version | awk 'BEGIN{i=1}{line[i++]=$0}END{j=1;while(j<4){printf"%s",line[j+1]line[j];j+=2};printf"\n"}'`
      lineh=${lineh/ highlight ver/HIGHLIGHT • Ver}
      lineh=${lineh/ Copyright / • }
      lineh=${lineh/Simon <andre.simon1 at gmx.de>/Simon (andre.simon1@gmx.de)}
      printlh
      printlh ' '$script' WILL ALSO USE AN ADDITIONAL UTILITY:'
      printlh
      printlh ' '${lineh/\,\ / • }
    fi
  fi
  printlh
  printlh ' '$script' UTILITIES AND SCRIPTS WILL BE TAKEN FROM:'
  printlh
  printlh ' '$path
  printlh
  waitingforinput='yes'
  printl "%s" '[10s] Press Enter to acknowledge...'; read -r -s -t 10 answer; printlh
  waitingforinput='no'
  printlh '_______________________________________________________________________________'

# Default parameters values.

  navecln=0   #<-- (= 0 seems optimal) tests showed that the method works best without averaging the end points of interpolation

# Default values for nsigma parameters (can be changed in the configuration file).
  
  verbosity=0; savespace=2; beam=10
  nscales=99; fscalemax=1.0; ffaintestscale=0
  k=1; nb=1; nwmax=0; l=0; la=0; nwa=0; nhead=0; lastline='no'; sfac=''

# Reading input parameters from '$configfile'.

  printlh
  printlh ' Reading input parameters from '$configfile'...'
  {
  while read -r -a v
  do
    if [[ "${v[0]:0:1}" != "#" ]]
    then 
      if [[ "${v[0]}" == "" ]]
      then 
        'ERROR' 'Line #1 of the file '$configfile' is empty.' \
                'Remove empty line(s) from the beginning of '$configfile' before running '$script
      fi
      l=$(( l + 1 ))
      if [[ $l -eq 1 ]]; then prefix=${v[0]}; fi
      if [[ $l -eq 2 ]]; then nwmax=${v[0]}; fi
      if [[ $l -eq 3 ]]
      then 
        nwaves=${v[0]}; nb=$(( l + 1 ))
        for (( i=1; i <= $nwaves; i++ ))
        do
          nw[i]=${v[i]}
          if [[ ${nw[i]:0:1} == "." || ${nw[i]:0:1} == "|" ]]
          then 
            'ERROR' 'Trouble reading line #'$(( nhead + 3 ))' from file '$configfile \
                    'Make sure that the line with numbers of wavelengths is correct before running '$script
          fi
        done
      fi
      if [[ $l -eq 4 ]]
      then
        sbem=0
        for (( i=1; i <= $nwmax; i++ ))
        do
          read -r vd1 vd2 vd3; l=$(( l + 1 )); read -r vo1 vo2; l=$(( l + 1 ))
          cfgwave[i]=${v[0]}; if [[ "${cfgwave[i]:0:1}" == "-" ]]; then cfgwave[i]=${cfgwave[i]:1}; fi
          cfgbeam[i]=${v[1]}; cfgsize[i]=${v[2]}; cfgdirp[i]=$vd1; cfgdirm[i]=$vd2
          if [[ ${v[1]:0:1} == "+" ]]; then cfgbeam[i]=${v[1]:1}; sbem=$i; fi
          if [[ ${cfgdirp[i]:0:2} == "~/" ]]; then cfgdirp[i]=~/${cfgdirp[i]:2}; fi
          if [[ ${cfgdirm[i]:0:2} == "~/" ]]; then cfgdirm[i]=~/${cfgdirm[i]:2}; fi
          if [[ ! -e ${cfgdirp[i]} ]]
          then
            'ERROR' 'Directory does not exist: '\'$vd1\' \
                    'Verify parameters in '$configfile' before running '$script
          fi
          cfgobsim[i]=${vo1/.fits/}
          if [[ $i -eq ${nw[k]} ]]
          then
            nwave[k]=$i; wave[k]=${v[0]}; beam[k]=${v[1]}
            if [[ ${v[1]:0:1} == "+" ]]; then beam[k]=${v[1]:1}; fi
            srcmaxsize[k]=${v[2]}; directory[k]=$vd1; dirmon[k]=$vd2; obsimage[k]=${vo1/.fits/}
            if [[ "${wave[k]:0:1}" == "-" ]]; then wave[k]=${wave[k]:1}; fi
            k=$(( k + 1 ))
          fi
          read -r -a v; l=$(( l + 1 ))
        done
      fi
      nx=$(( nb + nwmax * 3 - 1 ))
      if [[ ${v[0]} == "y"  ]]; then v[0]='yes'; fi
      if [[ ${v[0]} == "n"  ]]; then v[0]='no' ; fi
      if [[ $l -eq $(( nx + 1 )) ]]
      then 
        if [[ ${v[0]} == "y-" ]]; then v[0]='yes'; bgsub='-'; else bgsub=''; fi
        flattening=${v[0]}; sfac=${v[1]}
      fi
      if [[ $l -eq $(( nx + 2 )) ]]; then  savespace=${v[0]}; fi
      if [[ $l -eq $(( nx + 3 )) ]]; then  verbosity=${v[0]}; fi
    else
      nhead=$(( nhead + 1 )) 
    fi
    totalparams=3
    if [[ $l -eq $(( nx + totalparams )) && ${v[0]:0:4} == "#___" ]]; then lastline='yes'; fi
  done
  if [[ $lastline == "no" ]]
  then 
    'ERROR' 'Trouble finding end of the configuration file '$configfile \
            'Make sure that the total number of lines is '$(( nhead + nx + totalparams ))' before running '$script
  fi
  } < $configfile
#___________________________________________________________________________________________________________________________________

  verb=''
  if [[ $verbosity -eq 0 ]]; then verb='-verb0'; fi
  if [[ $verbosity -eq 1 ]]; then verb='-verb1'; fi
  if [[ $verbosity -eq 2 ]]; then verb='-verb2'; fi
  if [[ $verb == '' ]]
  then 
    'ERROR' 'Verbosity '\'$verbosity\'' out of range [0,1,2] ('$configfile')' \
            'Only these values of verbosity level are supported: 0, 1, or 2'
  fi

  if [[ $verb != '-verb0' ]]
  then uhigh=$highlight
  else uhigh='tee'
  fi
  
  if [[ ${sfac:0:3} != '<--' &&  ${sfac:0:2} != 'no' ]]
  then 
    'COMPARE_NUMBERS' $sfac '<' '1.05'; rslt=$?
  else
    rslt=1
  fi
  if [[ $rslt -eq 1 ]]
  then
    'ERROR' 'Invalid flattening parameter SFAC = '\'$sfac\'' in file '$configfile \
            'Make sure the quantity has a value of 1.05--1.5 before running '$script \
            'Save default configuration file to see valid options.'
  fi
    
  subdirbackups=$originpath'/+backups'
  subdircatalogs=$originpath'/+catalogs'
  subdirfilaments=$originpath'/+filaments'
  subdirinterface=$originpath'/+interface'
  subdiriterations=$originpath'/+iterations'
  subdirsources=$originpath'/+sources'
  subdirvisuals=$originpath'/+visuals'
  subdirannuli=$subdiriterations'/+annuli'
  subdirsprofiles=$subdirsources'/+profiles'

  finalcatalog=$prefix'.sw.final.relitent.cat'
  finellreliable=$subdirvisuals'/'$prefix'.sw.final.reliable.ell'
  finellrelitent=$subdirvisuals'/'$prefix'.sw.final.relitent.ell'
  findotreliable=$subdirvisuals'/'$prefix'.sw.final.reliable.dot'
  findotrelitent=$subdirvisuals'/'$prefix'.sw.final.relitent.dot'
  finfooreliable=$subdirvisuals'/'$prefix'.sw.final.reliable.foo'
  finfoorelitent=$subdirvisuals'/'$prefix'.sw.final.relitent.foo'

  for (( i=1; i <= $nwmax; i++ ))
  do
    for (( ii=1; ii <= $nwmax; ii++ ))
    do
      if [[ $ii != $i && "${cfgwave[ii]}" == "${cfgwave[i]}" ]]
      then 
        'ERROR' 'Configuration file contains images with the same wavelength: '${cfgwave[i]} \
                'Make sure the wavelengths are different in at least a single digit'
      fi
    done
    filmaxsize[i]=0
    filmaxsizeo[i]=0
  done

  for (( i=1; i <= $nwaves; i++ ))
  do
    if [[ ${directory[i]:0:2} == "~/" ]]; then directory[i]=~/${directory[i]:2}; fi
    if [[ ! -e ${directory[i]} ]]
    then 
      'ERROR' 'Directory not found: '${directory[i]} \
              'Correct the path in '$configfile
    fi
    if [[ ! -e ${directory[i]}'/'${obsimage[i]}'.fits' ]]
    then 
      'ERROR' 'File not found: '${directory[i]}'/'${obsimage[i]}'.fits' \
              'Correct the data in '$configfile
    fi
    imagepath[i]=`\cd ${directory[i]}; pwd`
  done

  imagesize=`\du -sk ${imagepath[1]}'/'${obsimage[1]}'.fits'`
  'POSITION' '/' "$imagesize"; ris=$returnposition
  if [[ $ris -eq 0 ]]
  then 
    'ERROR' 'Trouble finding size of ' ${imagepath[1]}'/'${obsimage[1]}'.fits'
  fi
  ris=$(( ris - 2 ))
  imagesize=${imagesize:0:ris}
  imagesize=$(print "scale=6; $imagesize/1024" | bc)
  if [[ ${imagesize:0:1} == "." ]]; then imagesize='0'$imagesize; fi
  printlh
  printlh ' Sizing input FITS images: '$imagesize' MB'

  'POSITION' "$homepath" "$originpath"; rchpath=$returnposition

  if [[ $rchpath -gt 0 ]]
  then
    lhpath=${#homepath}
    ihpath=$(( rchpath + lhpath - 1 ))
    tilde='~'
  else
    ihpath=0
    tilde=''
  fi
  for (( i=1; i <= $nwaves; i++ ))
  do
    if [[ ! -e ${dirmon[i]} ]]
    then 
      for (( i=1; i <= $nwaves; i++ ))
      do
        if [[ ! -e ${dirmon[i]} ]]
        then 
          'POSITION' '|' "${dirmon[i]}"; pos=$returnposition
          if [[ $pos -gt 0 ]]
          then 
            'ERROR' 'Invalid monochromatic directory in '$configfile': "'${dirmon[i]}'"' \
                    'Make sure the separator symbol "|" has a leading space: " |"'
          else
            'ERROR' 'Monochromatic directory in '$configfile' does not exist: "'${dirmon[i]}'"' \
                    'Make sure the entire directory structure has been created.'
          fi
        fi
      done
    fi
  done
  suffix=$prefix'.' 
  printlh
  printlh ' Current directory:'
  printlh ' '$tilde${originpath:ihpath}
  printlh
  printlh ' Flattening directory:'
  for (( i=1; i <= $nwaves; i++ ))
  do
    suffix=$suffix${nwave[i]}
    dirmono[i]=`\cd ${dirmon[i]}; pwd`
    printlh ' '$tilde${dirmono[i]:ihpath}
  done

  filmaxsizemax=0
  
  for (( ii=1; ii <= $nwmax; ii++ ))
  do
    if [[ "${cfgsize[ii]:0:1}" == "+" ]]
    then
      cfgsize[ii]=${cfgsize[ii]:1}
    fi
    cfgsize[ii]=${cfgsize[ii]%.*}   #<-- make the number integer
  done

  workingpath=$originpath

  'POSITION' '.' "$prefix"; dot=$returnposition
  
  if [[ $dot -gt 0 ]]
  then 
    'ERROR' 'Dots in the PREFIX '\'$prefix\'' not allowed ('$configfile')' \
            'PREFIX specifies a suitable name to be used for naming output files'
  fi
  if [[ $nwmax -ge 100 ]]
  then 
    'ERROR' 'Maximum number of wavelengths NWMAX exceeds '$nwmax' ('$configfile')' \
            'This version of '$script' is limited to maximum '$nwmax' wavebands'
  fi

# Initialize variables and arrays.

  nwavesm1=$(( nwaves - 1 )); nwavesm2=$(( nwaves - 2 ))

# Determine maximum beam size over all waves.

  cfgbeammin=999; cfgbeammax=0; cfgmaxsizemax=0

  for (( i=1; i <= $nwmax; i++ ))
  do
    'COMPARE_NUMBERS' ${cfgbeam[i]} '<' $cfgbeammin; rslt=$?
    if [[ $rslt -eq 1 ]]; then cfgbeammin=${cfgbeam[i]}; fi
    'COMPARE_NUMBERS' ${cfgbeam[i]} '>' $cfgbeammax; rslt=$?
    if [[ $rslt -eq 1 ]]; then cfgbeammax=${cfgbeam[i]}; fi
    'COMPARE_NUMBERS' ${cfgsize[i]} '>' $cfgmaxsizemax; rslt=$?
    if [[ $rslt -eq 1 ]]; then cfgmaxsizemax=${cfgsize[i]}; fi 
  done
    
  beamav=0; beammin=999; beammax=0; srcmaxsizebmax=0; filmaxsizebmax=0

  for (( i=1; i <= $nwaves; i++ ))
  do
    if [[ $sbem -gt 0 ]]; then beam[i]=${cfgbeam[sbem]}; fi

    if [[ "${beam[i]:0:1}" == "0" && "${beam[i]:0:2}" != "0." ]]
    then 
      'ERROR' 'Leading zeros in beam sizes not allowed ('$configfile')' \
              'Remove leading zeros from beam sizes; they should have a usual form'
    fi

    'COMPARE_NUMBERS' ${beam[i]} '<' $beammin; rslt=$?
    if [[ $rslt -eq 1 ]]; then beammin=${beam[i]}; fi
    'COMPARE_NUMBERS' ${beam[i]} '>' $beammax; rslt=$?
    if [[ $rslt -eq 1 ]]; then beammax=${beam[i]}; fi

    beamav=$(print "scale=10; $beamav+${beam[i]}/$nwaves" | bc); chk_rc
           
    if [[ "${srcmaxsize[i]:0:1}" == "+" ]]
    then
      srcmaxsize[i]=${srcmaxsize[i]:1}
    fi
    if [[ $sbem -gt 0 ]]
    then
      'COMPARE_NUMBERS' ${beam[i]} '>' ${srcmaxsize[i]}; rslt=$?
      if [[ $rslt -eq 1 ]]; then srcmaxsize[i]=${beam[i]}; fi
    fi
    srcmaxsizeb[i]=${srcmaxsize[i]}
    filmaxsizeb[i]=${filmaxsize[i]}
    beamb[i]=${beam[i]}

    'COMPARE_NUMBERS' ${srcmaxsizeb[i]} '>' $srcmaxsizebmax; rslt=$?
    if [[ $rslt -eq 1 ]]; then srcmaxsizebmax=${srcmaxsizeb[i]}; fi
    'COMPARE_NUMBERS' ${filmaxsizeb[i]} '>' $filmaxsizebmax; rslt=$?
    if [[ $rslt -eq 1 ]]; then filmaxsizebmax=${filmaxsizeb[i]}; fi
  done

# Diagnose problems with the lengths of wavelengths.
    
  lenwavo=0; lenmono=0
  for (( i=1; i <= $nwmax; i++ ))
  do
    lenwav=${#cfgwave[i]}
    if [[ $lenwavo -gt 0 && $lenwav -ne $lenwavo ]]
    then
      'ERROR' 'Lengths of the wavelengths differ in '$configfile
              'Use leading zeros to make sure the waves have the same length'
    fi
    lenwavo=$lenwav
    lenmon=${#cfgdirm[i]}
    if [[ $lenmono -gt 0 && $lenmon -ne $lenmono && $flattening != 'yes' ]]
    then 
      'ERROR' 'Lengths of the monochromatic directories differ in '$configfile
              'Use leading zeros to make sure the names have the same length'
    fi
    lenmono=$lenmon
  done
  lenbmx=${#beammax}
  lenszx=${#srcmaxsizebmax}
  lenfix=${#filmaxsizebmax}
  for (( i=1; i <= $nwaves; i++ ))
  do
    lenbm=${#beamb[i]}; ladbz=$(( lenbmx - lenbm ))
    lensz=${#srcmaxsizeb[i]}; ladsz=$(( lenszx - lensz ))
    lenfi=${#filmaxsizeb[i]}; ladfi=$(( lenfix - lenfi ))
    if [[ $ladbz -gt 0 ]]
    then 
      for (( l=1; l <= $ladbz; l++ )); do beamb[i]=' '"${beamb[i]}"; done
    fi
    if [[ $ladsz -gt 0 ]]
    then 
      for (( l=1; l <= $ladsz; l++ )); do srcmaxsizeb[i]=' '"${srcmaxsizeb[i]}"; done
    fi
    if [[ $ladfi -gt 0 ]]
    then 
      for (( l=1; l <= $ladfi; l++ )); do filmaxsizeb[i]=' '"${filmaxsizeb[i]}"; done
    fi
    measuremtclone[i]=$prefix'.'${wave[i]}'.obs'
    imageomask[i]=$prefix'.'${wave[i]}'.omask'
    backgroundclone[i]=$prefix'.'${wave[i]}'.background'
    flatfactorclone[i]=$prefix'.'${wave[i]}'.flattening'
  done
  
# Check the input images for correctness of the FITS keywords.
# Diagnose problems with the paths or files.

  for (( i=1; i <= $nwaves; i++ ))
  do 
    lendir=${#dirmono[i]}
    lendirm1=$(( lendir - 1 ))
    if [[ ${dirmono[i]:lendirm1} == "/" ]]
    then 
      'ERROR' 'Trailing slash "/" is found in the configuration path: '${dirmono[i]}
              'Trailing slashes should not be used in the names of directories'
    fi
    lendir=${#imagepath[i]}
    lendirm1=$(( lendir - 1 ))
    if [[ ${imagepath[i]:lendirm1} == "/" ]]
    then 
      'ERROR' 'Trailing slash "/" is found in the configuration path: '${imagepath[i]}
              'Trailing slashes should not be used in the names of directories'
    fi
    if [[ ${dirmono[i]:0:3} == "..." ]]
    then 
      'ERROR' 'Invalid name of the monochromatic extraction directory: '${dirmono[i]}
              'The name starts with at least three dots'
    fi
    if [[ ! -e ${dirmono[i]} ]]
    then 
      'ERROR' 'Monochromatic extraction directory not found: '${dirmono[i]}
              'The directory for source extraction must exist; check the path for typos'
    fi
    if [[ ! -e ${imagepath[i]} ]]
    then 
      'ERROR' 'Directory with input images not found: '${imagepath[i]}
              'The directory with prepared images must exist; check the path for typos'
    fi
    if [[ ! -e ${imagepath[i]}'/'${obsimage[i]}'.fits' ]]
    then 
      'ERROR' 'File with a prepared image not found: '${imagepath[i]}'/'${obsimage[i]}'.fits'
              'Prepared image must exist; check the path for typos'
    fi
    imagetype[i]='other'
  done

# Copy files to the current directory.
# Should only be done when processing images in the main steps.
  
  if [[ $flattening == 'yes' ]]
  then
    for (( i=1; i <= $nwaves; i++ ))
    do 
      im1=$(( i - 1 ))
      if [[ "$workingpath" == "${dirmono[i]}" ]]
      then
        printlh
        printlh ' Images directory:'
        printlh ' '${directory[i]}
        printlh ' Copying '${obsimage[i]}'.fits...'
        \cp -f ${imagepath[i]}'/'${obsimage[i]}'.fits' ${measuremtclone[i]}'.fits' 2>$hell
        printlh ' Copying '${obsimage[i]}'.omask.fits...'
        if [[ ! -e ${imagepath[i]}'/'${obsimage[i]}'.omask.fits' ]]
        then 
          'ERROR' 'Observational mask '${obsimage[i]}'.omask.fits not found.' \
                  ${imagepath[i]}
        fi
        \cp -f ${imagepath[i]}'/'${obsimage[i]}'.omask.fits' ${imageomask[i]}'.fits' 2>$hell
      else
        'ERROR' 'Current directory differs from monochromatic directory' \
                'Verify correctness of directory names in '$configfile
      fi
      
      ($path'/imgstat' savein -nomed ${imageomask[i]} ${imageomask[i]} $verb$cwavenum; echo $? > .rc) | 'tehi'
      
      read -r min max mea med tots sig skew kurt edgemean < '.+imgstat2'$cwavenum; chk_rc

      'CHECK_IMGSTAT' $min $max '' $script $fun $funo '1'
                                                 
      'COMPARE_NUMBERS' ${min:0:6} '=' '0.0000'; rslt0=$?
      'COMPARE_NUMBERS' ${max:0:6} '=' '1.0000'; rslt1=$?
      'COMPARE_NUMBERS' ${min:0:6} '=' ${max:0:6}; rslt2=$?
      
      if [[ ($rslt0 -eq 0 && $rslt2 -eq 0) || $rslt1 -eq 0 ]]
      then
        'ERROR' 'Invalid observational mask: min = '$min', max = '$max
                'Make a mask image with the minimum = 0.0000 and maximum = 1.0000'
      fi
    done
    
    'VERIFY_INPUT_IMAGES'; chk_rc
  fi
  
# Determine type of the images.

  answer=''
  
  for (( i=1; i <= $nwaves; i++ ))
  do
    if [[ ${wave[i]} == "070" || ${wave[i]} == "0070" || ${wave[i]} == "70" ]]
    then       
      if [[ "${beam[i]}" == "${b070s10:1}" || "${beam[i]}" == "${b070s20:1}" || "${beam[i]}" == "${b070s20p:1}" || \
            "${beam[i]}" == "${b070s60:1}" || "${beam[i]}" == "${b070s60p:1}" ]]
      then
        imagetype[i]='Herschel'
      else
        printlh
        printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        printlh ' INFO: The band at '${wave[i]}' microns is one of the Herschel PACS wavebands, but'
        printlh '       the beam ('${beam[i]}' arcsec) is not one of those adopted by GETIMAGES'
        printlh '       ('${b070s10:1}' '${b070s20:1}' '${b070s20p:1}' '${b070s60:1}' '${b070s60p:1}')'
        printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        waitingforinput='yes'
        printl "%s" 'QQQQ: Has the '${wave[i]}' micron image been observed by Herschel (y/N)?: '; read -r -s answer; rc=$?
        waitingforinput='no'
    
        if [[ ${answer:0:1} != "n" && ${answer:0:1} != "N" && ${answer:0:1} != "y" && ${answer:0:1} != "Y" && $answer != "" ]]
        then
          answer='n'; printl "%s" 'assuming '
        fi
        if [[ ${answer:0:1} == "n" || ${answer:0:1} == "N" || $answer == "" ]]
        then 
          answer='n'; printlh 'no!'
        fi
        if [[ ${answer:0:1} == "y" || ${answer:0:1} == "Y" ]]
        then 
          printlh 'yes!'; imagetype[i]='Herschel'
        fi
        sleep 1
      fi
    fi
    if [[ ${wave[i]} == "100" || ${wave[i]} == "0100" ]]
    then       
      if [[ "${beam[i]}" == "${b100s10:1}" || "${beam[i]}" == "${b100s20:1}" || "${beam[i]}" == "${b100s20p:1}" || \
            "${beam[i]}" == "${b100s60:1}" || "${beam[i]}" == "${b100s60p:1}" ]]
      then
        imagetype[i]='Herschel'
      else
        printlh
        printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        printlh ' INFO: The band at '${wave[i]}' microns is one of the Herschel PACS wavebands, but'
        printlh '       the beam ('${beam[i]}' arcsec) is not one of those adopted by GETIMAGES'
        printlh '       ('${b100s10:1}' '${b100s20:1}' '${b100s20p:1}' '${b100s60:1}' '${b100s60p:1}')'
        printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        waitingforinput='yes'
        printl "%s" 'QQQQ: Has the '${wave[i]}' micron image been observed by Herschel (y/N)?: '; read -r -s answer; rc=$?
        waitingforinput='no'
    
        if [[ ${answer:0:1} != "n" && ${answer:0:1} != "N" && ${answer:0:1} != "y" && ${answer:0:1} != "Y" && $answer != "" ]]
        then
          answer='n'; printl "%s" 'assuming '
        fi
        if [[ ${answer:0:1} == "n" || ${answer:0:1} == "N" || $answer == "" ]]
        then 
          answer='n'; printlh 'no!'
        fi
        if [[ ${answer:0:1} == "y" || ${answer:0:1} == "Y" ]]
        then 
          printlh 'yes!'; imagetype[i]='Herschel'
        fi
        sleep 1
      fi
    fi
    if [[ ${wave[i]} == "160" || ${wave[i]} == "0160" ]]
    then       
      if [[ "${beam[i]}" == "${b160s10}" || "${beam[i]}" == "${b160s20}" || "${beam[i]}" == "${b160s20p}" || \
            "${beam[i]}" == "${b160s60}" || "${beam[i]}" == "${b160s60p}" ]]
      then
        imagetype[i]='Herschel'
      else
        printlh
        printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        printlh ' INFO: The band at '${wave[i]}' microns is one of the Herschel PACS wavebands, but'
        printlh '       the beam ('${beam[i]}' arcsec) is not one of those adopted by GETIMAGES'
        printlh '       ('${b160s10}' '${b160s20}' '${b160s20p}' '${b160s60}' '${b160s60p}')'
        printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        waitingforinput='yes'
        printl "%s" 'QQQQ: Has the '${wave[i]}' micron image been observed by Herschel (y/N)?: '; read -r -s answer; rc=$?
        waitingforinput='no'
    
        if [[ ${answer:0:1} != "n" && ${answer:0:1} != "N" && ${answer:0:1} != "y" && ${answer:0:1} != "Y" && $answer != "" ]]
        then
          answer='n'; printl "%s" 'assuming '
        fi
        if [[ ${answer:0:1} == "n" || ${answer:0:1} == "N" || $answer == "" ]]
        then 
          answer='n'; printlh 'no!'
        fi
        if [[ ${answer:0:1} == "y" || ${answer:0:1} == "Y" ]]
        then 
          printlh 'yes!'; imagetype[i]='Herschel'
        fi
        sleep 1
      fi
    fi
    if [[ ${wave[i]} == "250" || ${wave[i]} == "0250" ]]
    then       
      if [[ "${beam[i]}" == "$b250" ]]
      then 
        imagetype[i]='Herschel'
      else
        printlh
        printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        printlh ' INFO: The band at '${wave[i]}' microns is one of the Herschel SPIRE wavebands, but'
        printlh '       the beam ('${beam[i]}' arcsec) is not the one adopted by GETIMAGES ('${b250}')'
        printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        waitingforinput='yes'
        printl "%s" 'QQQQ: Has the '${wave[i]}' micron image been observed by Herschel (y/N)?: '; read -r -s answer; rc=$?
        waitingforinput='no'
    
        if [[ ${answer:0:1} != "n" && ${answer:0:1} != "N" && ${answer:0:1} != "y" && ${answer:0:1} != "Y" && $answer != "" ]]
        then
          answer='n'; printl "%s" 'assuming '
        fi
        if [[ ${answer:0:1} == "n" || ${answer:0:1} == "N" || $answer == "" ]]
        then 
          answer='n'; printlh 'no!'
        fi
        if [[ ${answer:0:1} == "y" || ${answer:0:1} == "Y" ]]
        then 
          printlh 'yes!'; imagetype[i]='Herschel'
        fi
        sleep 1
      fi
    fi
    if [[ ${wave[i]} == "350" || ${wave[i]} == "0350" ]]
    then       
      if [[ "${beam[i]}" == "$b350" ]]
      then 
        imagetype[i]='Herschel'
      else
        printlh
        printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        printlh ' INFO: The band at '${wave[i]}' microns is one of the Herschel SPIRE wavebands, but'
        printlh '       the beam ('${beam[i]}' arcsec) is not the one adopted by GETIMAGES ('${b350}')'
        printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        waitingforinput='yes'
        printl "%s" 'QQQQ: Has the '${wave[i]}' micron image been observed by Herschel (y/N)?: '; read -r -s answer; rc=$?
        waitingforinput='no'
    
        if [[ ${answer:0:1} != "n" && ${answer:0:1} != "N" && ${answer:0:1} != "y" && ${answer:0:1} != "Y" && $answer != "" ]]
        then
          answer='n'; printl "%s" 'assuming '
        fi
        if [[ ${answer:0:1} == "n" || ${answer:0:1} == "N" || $answer == "" ]]
        then 
          answer='n'; printlh 'no!'
        fi
        if [[ ${answer:0:1} == "y" || ${answer:0:1} == "Y" ]]
        then 
          printlh 'yes!'; imagetype[i]='Herschel'
        fi
        sleep 1
      fi
    fi
    if [[ ${wave[i]} == "500" || ${wave[i]} == "0500" ]]
    then       
      if [[ "${beam[i]}" == "$b500" ]]
      then 
        imagetype[i]='Herschel'
      else
        printlh
        printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        printlh ' INFO: The band at '${wave[i]}' microns is one of the Herschel SPIRE wavebands, but'
        printlh '       the beam ('${beam[i]}' arcsec) is not the one adopted by GETIMAGES ('${b500}')'
        printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
        waitingforinput='yes'
        printl "%s" 'QQQQ: Has the '${wave[i]}' micron image been observed by Herschel (y/N)?: '; read -r -s answer; rc=$?
        waitingforinput='no'
    
        if [[ ${answer:0:1} != "n" && ${answer:0:1} != "N" && ${answer:0:1} != "y" && ${answer:0:1} != "Y" && $answer != "" ]]
        then
          answer='n'; printl "%s" 'assuming '
        fi
        if [[ ${answer:0:1} == "n" || ${answer:0:1} == "N" || $answer == "" ]]
        then 
          answer='n'; printlh 'no!'
        fi
        if [[ ${answer:0:1} == "y" || ${answer:0:1} == "Y" ]]
        then 
          printlh 'yes!'; imagetype[i]='Herschel'
        fi
        sleep 1
      fi
    fi
  done
  
# Useful zero file for later initializations.

  if [[ $flattening == "yes" ]]
  then 
    $path'/modfits' multiply 0 ${measuremtclone[1]} -o 'z0' $verb > $hell; chk_rc

    for (( i=1; i <= $nwaves; i++ ))
    do
      fileflattening[i]=${dirmono[i]}'/'${measuremtclone[i]}'.flattening'
    done
  fi
#___________________________________________________________________________________________________________________________________

# Cut off the '.fits' name extension, if present.

  image=${image/.fits/}
  prefix=${prefix/.fits/}

  'CLEAR_STDIN'
  printlh '_______________________________________________________________________________'
  printlh
  printlh '        GETIMAGES: INPUT PARAMETERS'
  printlh  
  printlh '       flattening:' $flattening$bgsub $sfac
  printlh '        savespace:' $savespace
  printlh '        verbosity:' ${verb:5:1}
  for (( i=1; i <= $nwaves; i++ ))
  do
    'LEADZEROS' 1 $i zi
    printlh '       obsimage'$zi': '${obsimage[i]}' ~> type: '"${imagetype[i]}"
  done
  for (( i=1; i <= $nwaves; i++ ))
  do
    'LEADZEROS' 1 $i zi
    printlh '     imageomask'$zi': ./'${imageomask[i]}
  done
  for (( i=1; i <= $nwaves; i++ ))
  do
    'LEADZEROS' 1 $i zi
    printlh ' measuremtclone'$zi': ./'${measuremtclone[i]} '~> beamsize: '"${beamb[i]}" '~> maxsize: '"${srcmaxsizeb[i]}"
  done
  printlh
  waitingforinput='yes'
  printl "%s" 'QQQQ: Are the above parameters correct (Y/n)?: '; answer=''; read -r -s -t 30 answer; rc=$?
  waitingforinput='no'

  if [[ ${answer:0:1} != "n" && ${answer:0:1} != "N" && ${answer:0:1} != "y" && ${answer:0:1} != "Y" && $answer != "" ]]
  then answer='n'; printl 'assuming '
  fi
  if [[ ${answer:0:1} == "n" || ${answer:0:1} == "N" ]]
  then answer='n'; printlh 'no!'; printlh; printlh ' Answer '\'$answer\'' received - '$script' aborted.'; printlh; exit 99
  fi
  if [[ ${answer:0:1} == "y" || ${answer:0:1} == "Y" || $answer == "" ]]
  then printlh 'yes!'
  fi
  printlh '_______________________________________________________________________________'

  fun=$funo
}
#___________________________________________________________________________________________________________________________________

VERIFY_FITS_HEADERS ()
{
  local funo=$fun; fun='VERIFY_FITS_HEADERS'
  if [[ $4 == "" ]]
  then 
    printlh
    'ERROR' 'Too few parameters: "'$1'" "'$2'" "'$3'" "'$4'"'
  fi
  local bunit; local bunitc; local bzero; local bzeroc; local bscale; local bscalec; local abort
  local fpath=$1; local origimage=$2; local cloneimage=$3; local wavelength=$4

  abort='no'

  if [[ ! -e $fpath'/'$origimage'.fits' ]]
  then 
    printlh
    'ERROR' 'File not found: '$fpath'/'$origimage'.fits'
  fi
    
# Check the image for correctness of the keyword BUNIT.
  
  $path'/readhead' save BUNIT $fpath'/'$origimage $verb &> $hell; rch=$?
  
  if [[ $rch -eq 0 ]]; then read -r bunitc bunit < '.+readhead'; chk_rc; \rm -f '.+readhead' 2>$hell; fi
  
  if [[ ${bunit:1:8} == '       ' || ${bunit:1:8} == 'UNKNOWN' || ${bunit:1:8} == 'unknown' || $rch -eq 2 ]]
  then 
    printlh
    'ERROR' "Missing/blank/bad value ${bunit:1:8} of BUNIT in '$origimage.fits'" \
            'Use MODFITS to correct BUNIT keyword and then run '$script
  fi
  if [[ ${bunit:1:6} != "MJy/sr" && ${bunit:1:8} != "no-units" ]]
  then 
    printlh
    'ERROR' "Invalid value ${bunit:1:8} of BUNIT in '$origimage.fits'" \
            'Use MODFITS to set BUNIT to MJy/sr and then run '$script
  fi

# Check the image for correctness of the keyword BZERO.
  
  $path'/readhead' save BZERO $fpath'/'$origimage $verb &> $hell; rch=$?
  
  if [[ $rch -eq 0 ]]; then read -r bzeroc bzero < '.+readhead'; chk_rc; \rm -f '.+readhead' 2>$hell; fi
       
  if [[ ${bzero:0:19} != '0.0000000000000E+00' && ${bzero:0:18} != '0.000000000000E+00' && \
        ${bzero:0:17} != '0.00000000000E+00' && ${bzero:0:16} != '0.000000000E+00' && ${bzero:0:15} != '0.0000000E+00' && \
        ${bzero:0:14} != '0.00000000E+00' && ${bzero:0:13} != '0.0000000E+00' && ${bzero:0:12} != '0.000000E+00' && \
        ${bzero:0:11} != '0.00000E+00' && ${bzero:0:10} != '0.0000E+00' && ${bzero:0:9} != '0.000E+00' ]]
  then 
    printlh
    'ERROR' "Invalid value $bzero of BZERO in '$origimage.fits'" \
            'Use MODFITS to set BZERO to 0.0 and then run '$script
  fi

# Check the image for correctness of the keyword BSCALE.

  $path'/readhead' save BSCALE $fpath'/'$origimage $verb &> $hell; rch=$?

  if [[ $rch -eq 0 ]]; then read -r bscalec bscale < '.+readhead'; chk_rc; \rm -f '.+readhead' 2>$hell; fi

  if [[ ${bscale:0:19} != '1.0000000000000E+00' && ${bscale:0:18} != '1.000000000000E+00' && \
        ${bscale:0:17} != '1.00000000000E+00' && ${bscale:0:16} != '1.000000000E+00' && ${bscale:0:15} != '1.0000000E+00' && \
        ${bscale:0:14} != '1.00000000E+00' && ${bscale:0:13} != '1.0000000E+00' && ${bscale:0:12} != '1.000000E+00' && \
        ${bscale:0:11} != '1.00000E+00' && ${bscale:0:10} != '1.0000E+00' && ${bscale:0:9} != '1.000E+00' ]]
  then 
    printlh
    'ERROR' "Invalid value $bscale of keyword BSCALE in '$origimage.fits'" $bscale \
            'Use MODFITS to set BSCALE to 1.0 and then run '$script
  fi

  $path'/readhead' save CDELT1 CDELT2 $fpath'/'$origimage $verb &> $hell; rch=$?
  
  if [[ $rch -eq 0 ]]
  then
    read -r cdeltc cdelt1 cdeltc cdelt2 < '.+readhead'; chk_rc; \rm -f '.+readhead' 2>$hell
  else
    printlh
    'ERROR' "Trouble reading the keyword CDELT1 or CDELT2 in '$origimage.fits'" \
            'Use READHEAD, MODFITS to correct the header and then run '$script
  fi

  $path'/readhead' save CRPIX1 CRPIX2 $fpath'/'$origimage $verb &> $hell; rch=$?
  
  if [[ $rch -eq 0 ]]
  then
    read -r crpixc crpix1 crpixc crpix2 < '.+readhead'; chk_rc; \rm -f '.+readhead' 2>$hell
  else
    printlh
    'ERROR' "Trouble reading the keyword CRPIX1 or CRPIX2 in '$origimage'" \
            'Use READHEAD, MODFITS to correct the header and then run '$script
  fi

  $path'/readhead' save CRVAL1 CRVAL2 $fpath'/'$origimage $verb &> $hell; rch=$?
  
  if [[ $rch -eq 0 ]]
  then
    read -r crvalc crval1 crvalc crval2 < '.+readhead'; chk_rc; \rm -f '.+readhead' 2>$hell
  else
    printlh
    'ERROR' "Trouble reading the keyword CRVAL1 or CRVAL2 in '$origimage.fits'" \
            'Use READHEAD, MODFITS to correct the header and then run '$script
  fi

  $path'/readhead' save CTYPE1 CTYPE2 $fpath'/'$origimage $verb &> $hell; rch=$?
  
  if [[ $rch -eq 0 ]]
  then
    read -r ctypec ctype1 ctypec ctype2 < '.+readhead'; chk_rc; \rm -f '.+readhead' 2>$hell
  else
    printlh
    'ERROR' "Trouble reading the keyword CTYPE1 or CTYPE2 in '$origimage.fits'" \
            'Use READHEAD, MODFITS to correct the header and then run '$script
  fi

  $path'/readhead' save RA DEC $fpath'/'$origimage $verb &> $hell; rch=$?
  
  if [[ $rch -eq 0 ]]
  then
    read -r rascenc rascen declinc declin < '.+readhead'; chk_rc; \rm -f '.+readhead' 2>$hell
  else
    printlh
    'ERROR' "Trouble reading the keyword RA or DEC in '$fpath/$origimage.fits'" \
            'Use READHEAD, MODFITS to correct the header and then run '$script
  fi

  if [[ $crval1 != $rascen ]]
  then
    'ERROR' "The keywords CRVAL1: $crval1 and RA: $rascen do not match in '$fpath/$origimage.fits'"  \
            "Use READHEAD, MODFITS to correct the header and then run $script. The easiest fix is" \
            "to multiply the image by 1 using MODFITS: modfits multiply 1 $origimage -o $origimage"
  fi

  if [[ $crval2 != $declin ]]
  then
    printlh 
    'ERROR' "The keywords CRVAL2: $crval2 and DEC: $declin do not match in '$fpath/$origimage.fits'" \
            "Use READHEAD, MODFITS to correct the header and then run $script. The easiest fix is" \
            "to multiply the image by 1 using MODFITS: modfits multiply 1 $origimage -o $origimage"
  fi

  if [[ ${ctype1:1:8} != "RA---TAN" && ${ctype1:1:8} != "RA---GLS" ]]
  then 
    printlh 
    'ERROR' "Invalid value of CTYPE1: $ctype1 in '$fpath/$origimage.fits'" \
            "Use MODFITS to make CTYPE1 = RA---TAN and then run $script."
  fi

  if [[ ${ctype2:1:8} != "DEC--TAN" && ${ctype2:1:8} != "DEC--GLS" ]]
  then 
    printlh 
    'ERROR' "Invalid value of CTYPE2: $ctype2 in '$fpath/$origimage.fits'" \
            "Use MODFITS to make CTYPE2 = DEC--TAN and then run $script."
  fi

  if [[ $abort == "yes" ]]; then printlh; fun=$funo; return 99; fi

  $path'/readhead' save WAVE $fpath'/'$origimage $verb &> $hell; rch=$?; \rm -f '.+readhead' 2>$hell
 
  if [[ $rch -eq 2 ]]
  then 
    printlh
    printlh ' '$script': '$fun': WARNING: Inserting missing WAVE keyword...'
 
    ($path'/modfits' key WAVE $wavelength $cloneimage -o $cloneimage $verb; echo $? > .rc) | 'tehi'
  fi

  fun=$funo
}
#___________________________________________________________________________________________________________________________________

VERIFY_INPUT_IMAGES ()
{
  local funo=$fun; fun='VERIFY_INPUT_IMAGES'

  printlh; printl "%s" 'Verifying that all images have been correctly prepared: '

  obscd1o=''; obscd2o=''; obscrp1o=''; obscrp2o=''; obscrv1o=''; obscrv2o=''; obsraso=''; obsdeto=''
  detcd1o=''; detcd2o=''; detcrp1o=''; detcrp2o=''; detcrv1o=''; detcrv2o=''; detraso=''; detdeto=''
  mskcd1o=''; mskcd2o=''; mskcrp1o=''; mskcrp2o=''; mskcrv1o=''; mskcrv2o=''; mskraso=''; mskdeto=''
  abort='no'

  for (( i=1; i <= $nwmax; i++ ))
  do 
    printl "%s"$i
    if [[ ${cfgdirp[i]:0:2} == "~/" ]]
    then cfgimpath[i]=`\cd ~/${cfgdirp[i]:2}; pwd`
    else cfgimpath[i]=`\cd ${cfgdirp[i]}; pwd`
    fi
    im1=$(( i - 1 ))

    'VERIFY_FITS_HEADERS' ${cfgimpath[i]} ${cfgobsim[i]} $prefix'.'${cfgwave[i]}'.obs' ${cfgwave[i]}; chk_rc

    obscd1=$cdelt1; obscd2=$cdelt2; obscrp1=$crpix1; obscrp2=$crpix2; obscrv1=$crval1; obscrv2=$crval2; obsras=$rascen; obsdec=$declin

    'VERIFY_FITS_HEADERS' ${cfgimpath[i]} ${cfgobsim[i]}'.omask' ${cfgobsim[i]}'.omask' ${cfgwave[i]}; chk_rc
    
    mskcd1=$cdelt1; mskcd2=$cdelt2; mskcrp1=$crpix1; mskcrp2=$crpix2; mskcrv1=$crval1; mskcrv2=$crval2; mskras=$rascen; mskdec=$declin

    if [[ $obscd1 != $mskcd1 || $obscd2 != $mskcd2 ]]
    then 
      printlh
      printlh ' '$script': '$fun': ERROR:'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' INFO: Keywords CDELT1 or CDELT2 do not match in input images in the directory'
      printlh '       '${cfgimpath[i]}
      printlh '       1st image: '${cfgobsim[i]}
      printlh '       2nd image: '${cfgobsim[i]}'.omask'
      printlh '       1st image: CDELT1:' $obscd1 'CDELT2:' $obscd2
      printlh '       2nd image: CDELT1:' $mskcd1 'CDELT2:' $mskcd2
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' HINT: Use READHEAD, MODFITS, or OPERATE to make images have the same CDELT'
      printlh '       The easiest fix is to copy the header using OPERATE, e.g.:'
      printlh '       operate '${cfgobsim[i]}' h '${cfgobsim[i]}'.omask -o '${cfgobsim[i]}'.omask'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      abort='yes'
    fi
    if [[ $obscrp1 != $mskcrp1 || $obscrp2 != $mskcrp2 ]]
    then 
      printlh
      printlh ' '$script': '$fun': ERROR:'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' INFO: Keywords CRPIX1 or CRPIX2 do not match in input images in the directory'
      printlh '       '${cfgimpath[i]}
      printlh '       1st image: '${cfgobsim[i]}
      printlh '       2nd image: '${cfgobsim[i]}'.omask'
      printlh '       1st image: CRPIX1:' $obscrp1 'CRPIX2:' $obscrp2
      printlh '       2nd image: CRPIX1:' $mskcrp1 'CRPIX2:' $mskcrp2
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' HINT: Use READHEAD, MODFITS, or OPERATE to make images have the same CRPIX'
      printlh '       The easiest fix is to copy the header using OPERATE, e.g.:'
      printlh '       operate '${cfgobsim[i]}' h '${cfgobsim[i]}'.omask -o '${cfgobsim[i]}'.omask'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      abort='yes'
    fi
    if [[ $obscrv1 != $mskcrv1 || $obscrv2 != $mskcrv2 ]]
    then 
      printlh
      printlh ' '$script': '$fun': ERROR:'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' INFO: Keywords CRVAL1 or CRVAL2 do not match in input images in the directory'
      printlh '       '${cfgimpath[i]}
      printlh '       1st image: '${cfgobsim[i]}
      printlh '       2nd image: '${cfgobsim[i]}'.omask'
      printlh '       1st image: CRVAL1:' $obscrv1 'CRVAL2:' $obscrv2
      printlh '       2nd image: CRVAL1:' $mskcrv1 'CRVAL2:' $mskcrv2
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' HINT: Use READHEAD, MODFITS, or OPERATE to make images have the same CRVAL'
      printlh '       The easiest fix is to copy the header using OPERATE, e.g.:'
      printlh '       operate '${cfgobsim[i]}' h '${cfgobsim[i]}'.omask -o '${cfgobsim[i]}'.omask'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      abort='yes'
    fi
    if [[ $obsras != $mskras || $obsdec != $mskdec ]]
    then 
      printlh
      printlh ' '$script': '$fun': ERROR:'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' INFO: Keywords RA or DEC do not match in input images in the directory'
      printlh '       '${cfgimpath[i]}
      printlh '       1st image: '${cfgobsim[i]}
      printlh '       2nd image: '${cfgobsim[i]}'.omask'
      printlh '       1st image: RA:' $obsras 'DEC:' $obsdec
      printlh '       2nd image: RA:' $mskras 'DEC:' $mskdec
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' HINT: Use READHEAD, MODFITS, or OPERATE to make images have the same RA DEC'
      printlh '       The easiest fix is to copy the header using OPERATE, e.g.:'
      printlh '       operate '${cfgobsim[i]}' h '${cfgobsim[i]}'.omask -o '${cfgobsim[i]}'.omask'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      abort='yes'
    fi
    if [[ $im1 -gt 0 && ($obscd1 != $obscd1o || $obscd2 != $obscd2o) ]]
    then 
      printlh
      printlh ' '$script': '$fun': ERROR:'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' INFO: Keywords CDELT1 or CDELT2 do not match in input images in the directory'
      printlh '       '${cfgimpath[i]}
      printlh '       1st image: '${cfgobsim[i]}
      printlh '       2nd image: '${cfgobsim[im1]}
      printlh '       1st image: CDELT1:' $obscd1  'CDELT2:' $obscd2
      printlh '       2nd image: CDELT1:' $obscd1o 'CDELT2:' $obscd2o
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' HINT: Use READHEAD, MODFITS, or OPERATE to make images have the same CDELT'
      printlh '       The easiest fix is to copy the header using OPERATE, e.g.:'
      printlh '       operate '${cfgobsim[i]}' h '${cfgobsim[im1]}' -o '${cfgobsim[im1]}
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      abort='yes'
    fi
    if [[ $im1 -gt 0 && ($mskcd1 != $mskcd1o || $mskcd2 != $mskcd2o) ]]
    then 
      printlh
      printlh ' '$script': '$fun': ERROR:'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' INFO: Keywords CDELT1 or CDELT2 do not match in input images in the directory'
      printlh '       '${cfgimpath[i]}
      printlh '       1st image: '${cfgobsim[i]}'.omask'
      printlh '       2nd image: '${cfgobsim[im1]}'.omask'
      printlh '       1st image: CDELT1:' $mskcd1  'CDELT2:' $mskcd2
      printlh '       2nd image: CDELT1:' $mskcd1o 'CDELT2:' $mskcd2o
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' HINT: Use READHEAD, MODFITS, or OPERATE to make images have the same CDELT#'
      printlh '       The easiest fix is to copy the header using OPERATE, e.g.:'
      printlh '       operate '${cfgobsim[i]}'.omask h '${cfgobsim[im1]}'.omask -o '${cfgobsim[im1]}'.omask'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      abort='yes'
    fi
    if [[ $im1 -gt 0 && ($obscrp1 != $obscrp1o || $obscrp2 != $obscrp2o) ]]
    then 
      printlh
      printlh ' '$script': '$fun': ERROR:'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' INFO: Keywords CRPIX1 or CRPIX2 do not match in input images in the directory'
      printlh '       '${cfgimpath[i]}
      printlh '       1st image: '${cfgobsim[i]}
      printlh '       2nd image: '${cfgobsim[im1]}
      printlh '       1st image: CRPIX1:' $obscrp1  'CRPIX2:' $obscrp2
      printlh '       2nd image: CRPIX1:' $obscrp1o 'CRPIX2:' $obscrp2o
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' HINT: Use READHEAD, MODFITS, or OPERATE to make images have the same CRPIX#'
      printlh '       The easiest fix is to copy the header using OPERATE, e.g.:'
      printlh '       operate '${cfgobsim[i]}' h '${cfgobsim[im1]}' -o '${cfgobsim[im1]}
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      abort='yes'
    fi
    if [[ $im1 -gt 0 && ($mskcrp1 != $mskcrp1o || $mskcrp2 != $mskcrp2o) ]]
    then 
      printlh
      printlh ' '$script': '$fun': ERROR:'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' INFO: Keywords CRPIX1 or CRPIX2 do not match in input images in the directory'
      printlh '       '${cfgimpath[i]}
      printlh '       1st image: '${cfgobsim[i]}'.omask'
      printlh '       2nd image: '${cfgobsim[im1]}'.omask'
      printlh '       1st image: CRPIX1:' $mskcrp1  'CRPIX2:' $mskcrp2
      printlh '       2nd image: CRPIX1:' $mskcrp1o 'CRPIX2:' $mskcrp2o
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' HINT: Use READHEAD, MODFITS, or OPERATE to make images have the same CRPIX#'
      printlh '       The easiest fix is to copy the header using OPERATE, e.g.:'
      printlh '       operate '${cfgobsim[i]}'.omask h '${cfgobsim[im1]}'.omask -o '${cfgobsim[im1]}'.omask'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      abort='yes'
    fi
    if [[ $im1 -gt 0 && ($obscrv1 != $obscrv1o || $obscrv2 != $obscrv2o) ]]
    then 
      printlh
      printlh ' '$script': '$fun': ERROR:'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' INFO: Keywords CRVAL1 or CRVAL2 do not match in input images in the directory'
      printlh '       '${cfgimpath[i]}
      printlh '       1st image: '${cfgobsim[i]}
      printlh '       2nd image: '${cfgobsim[im1]}
      printlh '       1st image: CRVAL1:' $obscrv1  'CRVAL2:' $obscrv2
      printlh '       2nd image: CRVAL1:' $obscrv1o 'CRVAL2:' $obscrv2o
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' HINT: Use READHEAD, MODFITS, or OPERATE to make images have the same CRVAL#'
      printlh '       The easiest fix is to copy the header using OPERATE, e.g.:'
      printlh '       operate '${cfgobsim[i]}' h '${cfgobsim[im1]}' -o '${cfgobsim[im1]}
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      abort='yes'
    fi
    if [[ $im1 -gt 0 && ($mskcrv1 != $mskcrv1o || $mskcrv2 != $mskcrv2o) ]]
    then 
      printlh
      printlh ' '$script': '$fun': ERROR:'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' INFO: Keywords CRVAL1 or CRVAL2 do not match in input images in the directory'
      printlh '       '${cfgimpath[i]}
      printlh '       1st image: '${cfgobsim[i]}'.omask'
      printlh '       2nd image: '${cfgobsim[im1]}'.omask'
      printlh '       1st image: CRVAL1:' $mskcrv1  'CRVAL2:' $mskcrv2
      printlh '       2nd image: CRVAL1:' $mskcrv1o 'CRVAL2:' $mskcrv2o
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' HINT: Use READHEAD, MODFITS, or OPERATE to make images have the same CRVAL#'
      printlh '       The easiest fix is to copy the header using OPERATE, e.g.:'
      printlh '       operate '${cfgobsim[i]}'.omask h '${cfgobsim[im1]}'.omask -o '${cfgobsim[im1]}'.omask'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      abort='yes'
    fi
    if [[ $im1 -gt 0 && ($obsras != $obsraso || $obsdec != $obsdeco) ]]
    then 
      printlh
      printlh ' '$script': '$fun': ERROR:'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' INFO: Keywords RA or DEC do not match in input images in the directory'
      printlh '       '${cfgimpath[i]}
      printlh '       1st image: '${cfgobsim[i]}
      printlh '       2nd image: '${cfgobsim[im1]}
      printlh '       1st image: RA:' $obsras  'DEC:' $obsdec
      printlh '       2nd image: RA:' $obsraso 'DEC:' $obsdeco
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' HINT: Use READHEAD, MODFITS, or OPERATE to make images have the same RA# DEC#'
      printlh '       The easiest fix is to copy the header using OPERATE, e.g.:'
      printlh '       operate '${cfgobsim[i]}' h '${cfgobsim[im1]}' -o '${cfgobsim[im1]}
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      abort='yes'
    fi
    if [[ $im1 -gt 0 && ($mskras != $mskraso || $mskdec != $mskdeco) ]]
    then 
      printlh
      printlh ' '$script': '$fun': ERROR:'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' INFO: Keywords RA or DEC do not match in input images in the directory'
      printlh '       '${cfgimpath[i]}
      printlh '       1st image: '${cfgobsim[i]}'.omask'
      printlh '       2nd image: '${cfgobsim[im1]}'.omask'
      printlh '       1st image: RA:' $mskras  'DEC:' $mskdec
      printlh '       2nd image: RA:' $mskraso 'DEC:' $mskdeco
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      printlh ' HINT: Use READHEAD, MODFITS, or OPERATE to make images have the same RA# DEC#'
      printlh '       The easiest fix is to copy the header using OPERATE, e.g.:'
      printlh '       operate '${cfgobsim[i]}'.omask h '${cfgobsim[im1]}'.omask -o '${cfgobsim[im1]}'.omask'
      printlh ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      abort='yes'
    fi

    obscd1o=$obscd1; obscd2o=$obscd2; obscrp1o=$obscrp1; obscrp2o=$obscrp2; obscrv1o=$obscrv1; obscrv2o=$obscrv2; obsraso=$obsras; obsdeco=$obsdec
    mskcd1o=$mskcd1; mskcd2o=$mskcd2; mskcrp1o=$mskcrp1; mskcrp2o=$mskcrp2; mskcrv1o=$mskcrv1; mskcrv2o=$mskcrv2; mskraso=$mskras; mskdeco=$mskdec
  done
  
  printlh ' OK'
  if [[ $abort == "yes" ]]; then printlh; fun=$funo; return 99; fi

  fun=$funo
}
#___________________________________________________________________________________________________________________________________

GET_IMAGE_SIZES ()
{
  local funo=$fun; fun='GET_IMAGE_SIZES'; title=$script': => '\'$fun\'; printlh; printlh ' '"$title"
  if [[ $flattening == "no" ]]; then fun=$funo; return 0; fi
  local image=$1; local i

# Get the numbers of pixels in both dimensions (from sky at the 1st wavelength).
  
  ($path'/readhead' save NAXIS1 NAXIS2 $image $verb; echo $? > .rc) | 'tehi'

  read -r naxis1 nx naxis2 ny < '.+readhead'; chk_rc; \rm -f '.+readhead' 2>$hell

  if [[ $nx -le $ny ]]; then pixmin=$nx; else pixmin=$ny; fi             
  if [[ $nx -ge $ny ]]; then pixmax=$nx; else pixmax=$ny; fi             

# Get the pixel sizes in both dimensions (in degrees).
  
  ($path'/readhead' save CDELT1 CDELT2 $image $verb; echo $? > .rc) | 'tehi'

  read -r cdelt1 dx cdelt2 dy < '.+readhead'; chk_rc; \rm -f '.+readhead' 2>$hell

# Avoid negative pixel sizes.

  if [[ ${dx:0:1} == "-" ]]; then dx=${dx:1}; fi                        

# Convert pixel sizes from scientific notation to the fixed floating-point format.

  lendx=${#dx}; lendy=${#dy}; dexpx=${dx:(-4)}; dexpy=${dy:(-4)}

  if [[ ${dexpx:0:1} == "E" || ${dexpx:0:1} == "e" ]]
  then dexpx=${dexpx:1}; dx=$(print "scale=20; ${dx:0:(lendx-4)}*10^$dexpx" | bc); chk_rc
  fi
  if [[ ${dexpy:0:1} == "E" || ${dexpy:0:1} == "e" ]]
  then dexpy=${dexpy:1}; dy=$(print "scale=20; ${dy:0:(lendy-4)}*10^$dexpy" | bc); chk_rc
  fi

# Convert pixel and image sizes in arcseconds and degrees.

     dxas=$(print "scale=20; $dx*3600"  | bc); chk_rc;    dyas=$(print "scale=20; $dy*3600"  | bc); chk_rc
  sidexas=$(print "scale=20; $nx*$dxas" | bc); chk_rc; sideyas=$(print "scale=20; $ny*$dyas" | bc); chk_rc
  sidexdg=$(print "scale=20; $nx*$dx"   | bc); chk_rc; sideydg=$(print "scale=20; $ny*$dy"   | bc); chk_rc
  sqpixas=$(print "scale=20; $dxas*$dyas" | bc); chk_rc
  if [[ ${dxas:0:1} == "." ]]; then dxas='0'$dxas; fi
  if [[ ${dyas:0:1} == "." ]]; then dyas='0'$dyas; fi
  if [[ ${sidexas:0:1} == "." ]]; then sidexas='0'$sidexas; fi
  if [[ ${sideyas:0:1} == "." ]]; then sideyas='0'$sideyas; fi
  if [[ ${sidexdg:0:1} == "." ]]; then sidexdg='0'$sidexdg; fi
  if [[ ${sideydg:0:1} == "." ]]; then sideydg='0'$sideydg; fi
  printlh
  printlh ' '"$script"': Images: '$nx' x '$ny' pixels, dx = '${dxas:0:8}', dy = '${dyas:0:8}' arcsec'
  printlh ' '"$script"': Images: '${sidexas:0:6}' x '${sideyas:0:6}' arcsec = '${sidexdg:0:9}' x '${sideydg:0:9}' degrees'

# Observational beam should not be smaller than 2 pixels in well-sampled images.

  for (( i=1; i <= $nwaves; i++ ))
  do
    dxas2=$(print "scale=20; 2.0*$dxas" | bc); chk_rc
    dyas2=$(print "scale=20; 2.0*$dyas" | bc); chk_rc
    'COMPARE_NUMBERS' ${beam[i]} '<=' $dxas2; rsltx=$?
    'COMPARE_NUMBERS' ${beam[i]} '<=' $dyas2; rslty=$?

    if [[ $rsltx -eq 1 || $rslty -eq 1 ]]
    then 
      printlh
      printlh ' '$script': WARNING: Image pixel is too large at '${wave[i]}' um:'
      printlh ' '$script': WARNING: beam: '${beam[i]}' arcsec, pixel: '${dxas:0:11}' arcsec'
      printlh ' '$script':  ADVICE: Think about it now! (timeout 1 minute)...'
      waitingforinput='yes'
      read -r -s -t 60 answer; rc=$?
      waitingforinput='no'
    fi
  done
  
  fun=$funo
}
#___________________________________________________________________________________________________________________________________

CLEAR_IMAGE_BORDERS ()
{
  local funo=$fun; fun='CLEAR_IMAGE_BORDERS'; title=$script': => '\'$fun\'
  local clrtype=$1; local width=$2; local imagenamein=$3; local imagenameout=$4
  if [[ $4 == "" ]]
  then 
    'ERROR' 'Too few parameters: "'$1'" "'$2'" "'$3'" "'$4'"'
  fi
  if [[ $verb != "-verb0" ]]
  then 
    printlh
    printlh ' '$fun': '$clrtype' '$width' '$imagenamein
  fi
  if [[ ! -e $imagenamein && ! -e $imagenamein'.fits' ]]
  then 
    'ERROR' 'File not found: '$imagenamein
  fi

# Clean several-pixel border around the image to clean it of border effects.

  ($path'/modfits' border  -$width $imagenamein  -o $imagenameout $verb; echo $? > .rc) | 'tehi'
  ($path'/modfits' $clrtype $width $imagenameout -o $imagenameout $verb; echo $? > .rc) | 'tehi'

  fun=$funo
}
#___________________________________________________________________________________________________________________________________

CONVOLVE_IMAGE ()
{                          
  local funo=$fun; fun='CONVOLVE_IMAGE'; title=$script': => '\'$fun\'
  local what1=$1; local what2=$2; local conkernexp=$3; local addpixels=$4; local beamsize=$5; local imagein=$6; local imageout=$7
  local maskimage=$8
  if [[ $7 == "" ]]
  then 
    'ERROR' 'Too few parameters: "'$1'" "'$2'" "'$3'" "'$4'" "'$5'" "'$6
  fi
  addpixels=${addpixels:0:10}
  beamsize=${beamsize:0:15}
  if [[ $verb != "-verb0" ]]
  then 
    printlh
    printlh ' '$fun': '$addpixels' '$beamsize' '$imagein
  fi
  imagein=${imagein/.fits/}; imageout=${imageout/.fits/}
  if [[ ! -e $imagein'.fits' ]]
  then 
    'ERROR' 'File not found: '$imagein
  fi

# Convolve the stdev image with a large smoothing beam; if requested, expand masked area of the image before convolution.

  if [[ $maskimage == "" ]]
  then
    if [[ $imagein != $imageout ]]
    then
      \cp -f $imagein'.fits' $imageout'.fits' 2>$hell
    fi
  else
    ($path'/expanda' 1 $imagein $maskimage -o $imageout $verb; echo $? > .rc) | 'tehi'
  fi
  if [[ $conkernexp == "0" || $conkernexp == '' ]]
  then addpix=$addpixels
  else addpix=$(print "scale=20; $addpixels+3*e((1/sqrt(($conkernexp)^2))*l(100))" | bc -l); chk_rc
  fi
  ($path'/modfits' $what1 $addpixels $imageout -o $imageout $verb; echo $? > .rc) | 'tehi'
  ($path'/fftconv' $beamsize $conkernexp $imageout -o $imageout $verb; echo $? > .rc) | 'tehi'
  ($path'/modfits' $what2 -$addpixels $imageout -o $imageout $verb; echo $? > .rc) | 'tehi'

  fun=$funo
}
#___________________________________________________________________________________________________________________________________

FLATTEN_IMAGES ()
{                          
  local funo=$fun; fun='FLATTEN_IMAGES'; title=$script': => (1) '\'$fun\'
  if [[ $flattening == "no" ]]
  then printlh; printlh ' '"$title"' ~~~~~~~> skipping...'; fun=$funo; return 0
  else printlh; printlh ' '"$title"
  fi
  median_beg=`date +%s`; cdate=`date +%d\ %b\ %Y\ %a\ %H:%M:%S\ %Z`
  local nscales=0; local sfactor; local usm=''; local msu=''; local con=''; local pos=''
  local origpath; local image; local maskimage; local i; local j; local zj; local nparams=${#*}
  local nmust=2
  if [[ $nparams -lt $nmust ]]
  then 
    'ERROR' 'Too few parameters: "'$1'" "'$2'"'
  fi
  nprogress=0
  if [[ "$1" != "0" ]]; then nscales=$1; fi
  if [[ "$2" != "0" ]]; then sfactor=$2; fi

# Estimate total computation time necessary for processing.

  printlh
  printl ' '$script': Estimating approximate execution time: '
  itmax=2; total_time=0

  for (( i=1; i <= $nwaves; i++ ))
  do    
    mfscale[0]=$(print "scale=20; ${beam[i]}/$sfactor" | bc); chk_rc
    'COMPARE_NUMBERS' ${srcmaxsizeb[i]} '<' ${beam[i]}; rsltdet=$?
    if [[ $rsltdet -eq 1 ]]
    then lastscale=$(print "scale=20; 4.0*${beam[i]}" | bc); chk_rc
    else lastscale=$(print "scale=20; 4.0*${srcmaxsizeb[i]}" | bc); chk_rc
    fi
    lastscalem=$(print "scale=20; 0.9*$lastscale" | bc); chk_rc
    j1std=1
    beamstdx=$(print "scale=20; 2.0*${beam[i]}" | bc); chk_rc
    'COMPARE_NUMBERS' $beamstdx '>' $lastscale; rslt=$?
    if [[ $rslt -eq 1 ]]
    then beamstdx=$lastscale
    fi
    for (( j=1; j <= 100; j++ ))
    do 
      jm1=$(( j - 1 )); nscales=$j
      mfscale[j]=$(print "scale=20; $sfactor*${mfscale[jm1]}" | bc); chk_rc
      mfscale[j]=${mfscale[j]}'0000'; mfscale[j]=${mfscale[j]:0:6}
      'COMPARE_NUMBERS' ${mfscale[j]} '>' $lastscalem; rsltd=$?
      if [[ $rsltd -eq 1 ]]
      then 
        mfscale[j]=$lastscale'0000'; mfscale[j]=${mfscale[j]:0:6}
        break
      fi
      'COMPARE_NUMBERS' ${mfscale[j]} '<=' $beamstdx; rsltbeamstdx=$?
      if [[ $rsltbeamstdx -eq 1 ]]
      then j1std=$(( j + 1 ))
      fi
    done
    for (( j=1; j <= $nscales; j++ ))
    do 
      printl '·'
      mf_beg=`date +%s`
      radius=$(print "scale=20; ${mfscale[j]}/$dxas" | bc); chk_rc; radius=${radius:0:15}
      ($path'/imgstat' median $radius ${measuremtclone[i]} -o ${measuremtclone[i]}'.test' $verb; echo $? > .rc) | 'tehi'
      timeratio=$(print "scale=20; (${mfscale[j]}/${mfscale[1]})^2" | bc); chk_rc
      mf_end=`date +%s`; mf_time=$(( $mf_end - $mf_beg ))
      if [[ $mf_time -ge 9 || $j -eq $nscales ]]
      then break
      fi
    done
    printf "%b" "\033["$j"D\033[K"
    \rm -f ${measuremtclone[i]}'.test.fits' 2>$hell
    mf_time1=$(print "scale=20; $mf_time/$timeratio" | bc); chk_rc
    backgrd_time=0
    for (( j=1; j <= $nscales; j++ ))
    do 
      scale_time[j]=$(print "scale=20; $mf_time1*(${mfscale[j]}/${mfscale[1]})^2" | bc); chk_rc
      backgrd_time=$(print "scale=20; $backgrd_time+${scale_time[j]}" | bc); chk_rc
    done
    backgrd_time=$(print "scale=20; $backgrd_time+($itmax-1)*${scale_time[nscales]}" | bc); chk_rc
    if [[ $bgsub == '' ]]
    then total_time=$(print "scale=20; $total_time+$backgrd_time" | bc); chk_rc
    fi
    unused_time=0
    for (( j=1; j < $j1std; j++ ))
    do 
      unused_time=$(print "scale=20; $unused_time+${scale_time[j]}" | bc); chk_rc
    done
    total_time=$(print "scale=20; $total_time+($backgrd_time-$unused_time)" | bc); chk_rc
  done
  total_minutes=$(print "scale=2; $total_time/60" | bc); chk_rc
  total_hours=$(print "scale=2; $total_time/3600" | bc); chk_rc
  total_days=$(print "scale=2; $total_time/3600/24" | bc); chk_rc
  if [[ ${total_minutes:0:1} == '.' ]]; then total_minutes='0'$total_minutes; fi
  if [[ ${total_hours:0:1} == '.' ]]; then total_hours='0'$total_hours; fi
  if [[ ${total_days:0:1} == '.' ]]; then total_days='0'$total_days; fi
  printf "%b" "\033[900D\033[K"
  printlh ' '$script': Expected runtime: '${total_time%.*}' seconds = '$total_minutes' minutes = '$total_hours' hours = '$total_days' days'

# Median filter all images across spatial scales.

  for (( i=1; i <= $nwaves; i++ ))
  do    
    if [[ $originpath != ${dirmono[i]} ]]
    then origpath=${dirmono[i]}
    else origpath='.'
    fi
    mfscale[0]=$(print "scale=20; ${beam[i]}/$sfactor" | bc); chk_rc
    'COMPARE_NUMBERS' ${srcmaxsizeb[i]} '<' ${beam[i]}; rsltdet=$?
    if [[ $rsltdet -eq 1 ]]
    then lastscale=$(print "scale=20; 4.0*${beam[i]}" | bc); chk_rc
    else lastscale=$(print "scale=20; 4.0*${srcmaxsizeb[i]}" | bc); chk_rc
    fi
    lastscalem=$(print "scale=20; 0.9*$lastscale" | bc); chk_rc
    for (( j=1; j <= 100; j++ ))
    do 
      jm1=$(( j - 1 )); nscales=$j
      mfscale[j]=$(print "scale=20; $sfactor*${mfscale[jm1]}" | bc); chk_rc
      mfscale[j]=${mfscale[j]}'0000'; mfscale[j]=${mfscale[j]:0:6}
      sbeamx=${mfscale[j]}; sbeamx=${sbeamx%.*}
      'COMPARE_NUMBERS' ${mfscale[j]} '>' $lastscalem; rsltd=$?
      if [[ $rsltd -eq 1 ]]
      then 
        mfscale[j]=$lastscale'0000'; mfscale[j]=${mfscale[j]:0:6}
        sbeamx=${mfscale[j]}; sbeamx=${sbeamx%.*}
        break
      fi
    done
    image=${measuremtclone[i]}
    maskimage=${imageomask[i]}
    bgimage=$image'.bg.s'$sbeamx'as'
    bgsimage=$image'.bgs.s'$sbeamx'as'
    stdimage=$bgsimage'.std'
    bgsiflat=$bgsimage'.flat'
    
    beampixs=$(print "scale=20; ${beam[i]}/$dxas" | bc); chk_rc
    if [[ $bgsub == '' ]]
    then print '#'"${headerinfo[1]:3}" > ${fileflattening[i]}
    else print '#'"${headerinfo[1]:3}" >> ${fileflattening[i]}
    fi
    for (( v=2; v <= 4; v++ )); do print '#'"${headerinfo[v]:3}" >> ${fileflattening[i]}; done
    print '#' >> ${fileflattening[i]}
    print '#' $cdate >> ${fileflattening[i]}
    print '#' >> ${fileflattening[i]}
    print '# Wavelength: '${wave[i]}' um, beam size: '${beam[i]}' arcsec = '$beampixs' pixels' >> ${fileflattening[i]}
    print '# Median filtering from the smallest scale up to the maximum scale of '$lastscale' arcsec' >> ${fileflattening[i]}
    print '# Observed images are median-filtered with a scale factor of '$sfactor >> ${fileflattening[i]}
    print '#' >> ${fileflattening[i]}
    print '# Median filtering ~> '$bgimage'.fits' >> ${fileflattening[i]}

    if [[ $verb == "-verb0" ]]
    then 
      if [[ $nwaves -eq 1 ]]; then iw=''; else iw=$i; fi
      'PROGRESS_BAR' "$script: FLATTENING$iw" 0 0 $nscales 79 "·" ''
    fi
    ($path'/modfits' multiply 1 $maskimage -o $maskimage $verb; echo $? > .rc) | 'tehi'
    ($path'/expanda' 1 $image $maskimage -o $image $verb; echo $? > .rc) | 'tehi'

# Perform median filtering to estimate background.

    for (( j=1; j <= $nscales; j++ ))
    do 
      'LEADZEROS' 1 $j zj
      imagessc=$image'.'$zj'.c'${mfscale[j]}'as'
      if [[ $verb == "-verb0" ]]
      then 
        if [[ $nwaves -eq 1 ]]; then iw=''; else iw=$i; fi
        'PROGRESS_BAR' "$script: FLATTENING$iw" $j 0 $nscales 79 "·" ''
      fi
      if [[ $bgsub == '' ]]
      then
        radius=$(print "scale=20; ${mfscale[j]}/$dxas" | bc); chk_rc; radius=${radius:0:15}
        ($path'/expanda' 1 $image $maskimage -o $imagessc'.mf' $verb; echo $? > .rc) | 'tehi'
        for (( it=1; it <= $itmax; it++ ))
        do 
          ($path'/imgstat' median $radius $imagessc'.mf' -o $imagessc'.mf' $verb; echo $? > .rc) | 'tehi'
          print \'${dirmono[i]}'/'$imagessc'.mf.fits'\' \
                 'command: '\''imgstat'\'' median' $radius $imagessc'.mf' -o $imagessc'.mf' >> ${fileflattening[i]}
          if [[ $j -gt 1 ]]
          then ($path'/operate' $imagessc'.mf' l $image'.bg' -o $imagessc'.mf' $verb; echo $? > .rc) | 'tehi'
          fi
          if [[ $j -lt $nscales ]]; then break; fi
        done
        if [[ $j -eq 1 ]]
        then \cp -f $imagessc'.mf.fits' $image'.bg.fits' 2>$hell
        else ($path'/operate' $imagessc'.mf' l $image'.bg' -o $image'.bg' $verb; echo $? > .rc) | 'tehi'
        fi
        \cp -f $image'.bg.fits' $imagessc'.bg.fits' 2>$hell
        ($path'/operate' $image - $imagessc'.bg' -o $imagessc'.bgs' $verb; echo $? > .rc) | 'tehi'
      fi
    done                               

# Final median-filtered background and background-subtracted images.

    if [[ $bgsub == '' ]]
    then
      \cp -f $image'.bg.fits' $bgimage'.fits' 2>$hell
      ($path'/operate' $image - $bgimage -o $bgsimage $verb; echo $? > .rc) | 'tehi'
      if [[ $savespace -gt 0 ]]
      then \rm -f $image.??.c??????as*.[mb][fg]*.fits 2>$hell; \rm -f $image'.bg.fits' 2>$hell
      fi
    fi
    if [[ -e '+stopitnow!' ]]
    then print '# User requested to exit prematurely... done.' >> ${fileflattening[i]}
    fi

# Flatten background-subtracted images obtained above by median filtering.

    print '# Flattening ~> '$bgsiflat'.fits' >> ${fileflattening[i]}

    if [[ $verb == "-verb0" ]]
    then 
      if [[ $nwaves -eq 1 ]]; then iw=''; else iw=$i; fi
      'PROGRESS_BAR' "$script: FLATTENING$iw" 0 0 $nscales 79 "-" ''
    fi
    \cp -f $bgsimage'.fits' $stdimage'.fits' 2>$hell
    ($path'/imgstat' stdev 1 1 $stdimage $maskimage -o $stdimage $verb; echo $? > .rc) | 'tehi'

    for (( j=1; j <= $nscales; j++ ))
    do 
      'LEADZEROS' 1 $j zj
      imagessc=$stdimage'.'$zj'.c'${mfscale[j]}'as'
      if [[ $verb == "-verb0" ]]
      then 
        if [[ $nwaves -eq 1 ]]; then iw=''; else iw=$i; fi
        if [[ $j -le $nscales ]]
        then 'PROGRESS_BAR' "$script: FLATTENING$iw" $j 0 $nscales 79 "·"
        fi
      fi
      if [[ $j -ge $j1std ]]
      then
        radius=$(print "scale=20; ${mfscale[j]}/$dxas" | bc); chk_rc; radius=${radius:0:15}
        ($path'/expanda' 1 $stdimage $maskimage -o $imagessc'.std.mf' $verb; echo $? > .rc) | 'tehi'
        for (( it=1; it <= $itmax; it++ ))
        do 
          ($path'/imgstat' median $radius $imagessc'.std.mf' -o $imagessc'.std.mf' $verb; echo $? > .rc) | 'tehi'
          print \'${dirmono[i]}'/'$imagessc'.std.mf.fits'\' \
                 'command: '\''imgstat'\'' median' $radius $imagessc'.std.mf' -o $imagessc'.std.mf' >> ${fileflattening[i]}
          if [[ $j -gt $j1std ]]
          then ($path'/operate' $imagessc'.std.mf' l $stdimage'.bg' -o $imagessc'.std.mf' $verb; echo $? > .rc) | 'tehi'
          fi
          if [[ $j -lt $nscales ]]; then break; fi
        done
        if [[ $j -eq $j1std ]]
        then \cp -f $imagessc'.std.mf.fits' $stdimage'.bg.fits' 2>$hell
        else ($path'/operate' $imagessc'.std.mf' l $stdimage'.bg' -o $stdimage'.bg' $verb; echo $? > .rc) | 'tehi'
        fi
        \cp -f $stdimage'.bg.fits' $imagessc'.std.bg.fits' 2>$hell
        ($path'/operate' $stdimage - $imagessc'.std.bg' -o $imagessc'.std.bgs' $verb; echo $? > .rc) | 'tehi'
      fi
    done                               

# Rename and create background-subtracted detection and stdev images.
    
    \mv -f $stdimage'.bg.fits' $stdimage'.bg.s'$sbeamx'as.fits' 2>$hell
    ($path'/operate' $stdimage - $stdimage'.bg.s'$sbeamx'as' -o $stdimage'.bgs.s'$sbeamx'as' $verb; echo $? > .rc) | 'tehi'
    ($path'/operate' $bgsimage % $stdimage'.bg.s'$sbeamx'as' -o $bgsiflat $verb; echo $? > .rc) | 'tehi'
    ($path'/imgstat' stdev 1 1 $bgsiflat $maskimage -o $bgsiflat'.std' $verb; echo $? > .rc) | 'tehi'
  
    if [[ $savespace -gt 0 ]]
    then \rm -f $stdimage.??.c??????as*.[mb][fg]*.fits 2>$hell
    fi

# Smooth the background-subtracted and flattened detection image and name the observational mask correspondingly.

    effres=$(print "scale=20; 1.05*${beam[i]}" | bc); chk_rc
    smbeam=$(print "scale=20; sqrt($effres^2-${beam[i]}^2)" | bc); chk_rc
    addpix=$(print "scale=20; 3.0*$smbeam/$dxas" | bc); chk_rc
  
    if [[ ${effres:0:1} == '.' ]]; then effres='0'${effres:0:4}; fi
    if [[ ${effres:0:4} == '0.00' ]]
    then effres=${effres:0:7}; elif [[ ${effres:0:3} == '0.0' ]]
    then effres=${effres:0:6}; elif [[ ${effres:0:2} == '0.' ]]
    then effres=${effres:0:5}; else effres=${effres:0:5}
    fi
    effres=${effres/./p}

    'CONVOLVE_IMAGE' expand border 0 $addpix $smbeam $bgsiflat $bgsiflat'.r'$effres $maskimage; chk_rc

    \cp -f $maskimage'.fits' $bgsiflat'.r'$effres'.omask.fits' 2>$hell

    if [[ $verb == "-verb0" ]]
    then 'PROGRESS_BAR' 'complete' $nscales 0 $nscales 79 "·"
    fi
  done
  
  median_end=`date +%s`; median_time=$(( median_end - median_beg ))
  printlh; printlh ' GETIMAGES: '$median_time' SECONDS'
  fun=$funo
} 
#___________________________________________________________________________________________________________________________________

FINISH ()
{                                         
  local funo=$fun; fun='FINISH'; title=$script': => (2) '\'$fun\'
  printlh; printlh ' '$title; if [[ $verb == '-verb2' ]]; then printlh; printlh '   Removing temporary files...' ; fi
  finish_beg=`date +%s`
  local i; local j; local noisebgr=0; local position=0; local notdetec=0; local spurious=0; local substruc=0
  local totalbad=0; local totalgod=0
                  
# Make sure all the files are readable by all (sometimes on some systems they are not).

  chmod a+r *
  \rm -f '.rc' 2>$hell

  finish_end=`date +%s`; finish_time=$(( finish_end - finish_beg ))
  printlh; printlh ' GETIMAGES: '$finish_time' SECONDS'
  fun=$funo
}
#___________________________________________________________________________________________________________________________________
#___________________________________________________________________________________________________________________________________
#___________________________________________________________________________________________________________________________________

'GET_VERSION' $*; chk_rc

for (( b=1; b <= 7; b++ )); do printlh "${banner[b]}"; done

getimages_beg=`date +%s`

'GET_PARAMETERS' $*; chk_rc

'GET_IMAGE_SIZES' ${measuremtclone[1]}; chk_rc

'FLATTEN_IMAGES' $nscales $sfac; chk_rc

'FINISH'; chk_rc

getimages_end=`date +%s`
getimages_time=$(( getimages_end - getimages_beg ))
cdate=`date +%d\ %b\ %Y\ %a\ %H:%M:%S\ %Z`
gets_minutes=$(print "scale=3; $getimages_time/60" | bc); chk_rc
gets_hours=$(print "scale=3; $getimages_time/3600" | bc); chk_rc
gets_days=$(print "scale=3; $getimages_time/3600/24" | bc); chk_rc
if [[ ${gets_minutes:0:1} == '.' ]]; then gets_minutes='0'$gets_minutes; fi
if [[ ${gets_hours:0:1} == '.' ]]; then gets_hours='0'$gets_hours; fi
if [[ ${gets_days:0:1} == '.' ]]; then gets_days='0'$gets_days; fi
printlh
printlh ' '$script': DONE IN '$getimages_time' SECONDS = '$gets_minutes' MINUTES = '$gets_hours' HOURS = '$gets_days' DAYS'
printlh
printlh ' '$cdate
printlh
exit
#___________________________________________________________________________________________________________________________________
