/* -*- Mode: C; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
/* vim:set sw=2 sts=2 et cin: */
/*
 * This file is part of the MUSE Instrument Pipeline
 * Copyright (C) 2005-2018 European Southern Observatory
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

/*----------------------------------------------------------------------------*
 *                              Includes                                      *
 *----------------------------------------------------------------------------*/
#include <stdio.h>
#include <float.h>
#include <math.h>
#include <string.h>
#include <cpl.h>

#include "muse_basicproc.h"

#include "muse_combine.h"
#include "muse_pfits.h"
#include "muse_quadrants.h"
#include "muse_quality.h"
#include "muse_utils.h"
#include "muse_data_format_z.h"

/*----------------------------------------------------------------------------*
 *                             Debugging Macros                               *
 *         Set these to 1 or higher for (lots of) debugging output            *
 *----------------------------------------------------------------------------*/
#define GENERATE_TEST_IMAGES 0 /* generate FITS file(s) to be used for testing */

/*---------------------------------------------------------------------------*/
/**
 * @defgroup muse_basicproc    Basic-processing function
 *
 * This group contains functions that are common to some basic-processing
 * recipes.
 */
/*---------------------------------------------------------------------------*/

/**@{*/

/*---------------------------------------------------------------------------*/
/**
  @brief  Create a new structure of basic processing parameters.
  @param  aParameters   the list of parameters
  @param  aPrefix       the prefix of the recipe
  @return The basic processing parameters as muse_basicproc_params or NULL on
          error.

  Create a new set of basic processing parameters from a recipe parameter list.
  It takes the parameters "overscan", "ovscreject", "ovscsigma", and
  "ovscignore" and converts them into the structure.
  If aPrefix contains "muse_scibasic", it also takes the parameters "cr",
  "xbox", "ybox", "passes", and "thres".

  The error handling is done by the cpl_parameterlist functions that this
  function calls.

  The component "keepflat" is set to false, so that by default the flat-field
  image is not stored in the "flatimage" component.

  Use muse_basicproc_params_delete() to free memory allocated by this function.
 */
/*---------------------------------------------------------------------------*/
muse_basicproc_params *
muse_basicproc_params_new(cpl_parameterlist *aParameters, const char *aPrefix)
{
  muse_basicproc_params *bpars = cpl_calloc(1, sizeof(muse_basicproc_params));
  cpl_parameter *param;
  param = muse_cplparamerterlist_find_prefix(aParameters, aPrefix, "overscan");
  bpars->overscan = cpl_strdup(cpl_parameter_get_string(param));
  param = muse_cplparamerterlist_find_prefix(aParameters, aPrefix, "ovscreject");
  bpars->rejection = cpl_strdup(cpl_parameter_get_string(param));
  param = muse_cplparamerterlist_find_prefix(aParameters, aPrefix, "ovscsigma");
  cpl_errorstate state = cpl_errorstate_get();
  bpars->ovscsigma = cpl_parameter_get_double(param);
  if (!cpl_errorstate_is_equal(state)) { /* try again as int, may be misidentified */
    cpl_errorstate_set(state);
    bpars->ovscsigma = cpl_parameter_get_int(param);
  }
  param = muse_cplparamerterlist_find_prefix(aParameters, aPrefix, "ovscignore");
  bpars->ovscignore = cpl_parameter_get_int(param);

  if (strstr(aPrefix, "muse_scibasic")) {
    param = muse_cplparamerterlist_find_prefix(aParameters, aPrefix, "cr");
    bpars->crmethod = cpl_strdup(cpl_parameter_get_string(param));
    param = muse_cplparamerterlist_find_prefix(aParameters, aPrefix, "xbox");
    bpars->dcrxbox = cpl_parameter_get_int(param);
    param = muse_cplparamerterlist_find_prefix(aParameters, aPrefix, "ybox");
    bpars->dcrybox = cpl_parameter_get_int(param);
    param = muse_cplparamerterlist_find_prefix(aParameters, aPrefix, "passes");
    bpars->dcrpasses = cpl_parameter_get_int(param);
    param = muse_cplparamerterlist_find_prefix(aParameters, aPrefix, "thres");
    state = cpl_errorstate_get();
    bpars->dcrthres = cpl_parameter_get_double(param);
    if (!cpl_errorstate_is_equal(state)) { /* try again as int */
      cpl_errorstate_set(state);
      bpars->dcrthres = cpl_parameter_get_int(param);
    }
  } /* if scibasic prefix */

  /* components keepflat and flatimage are automatically CPL_FALSE and NULL */
  return bpars;
} /* muse_basicproc_params_new() */

/*---------------------------------------------------------------------------*/
/**
  @brief  Create a structure of basic processing parameters from a FITS header.
  @param  aHeader   the FITS header
  @return The basic processing parameters as muse_basicproc_params or NULL on
          error.

  This creates a new parameter list from the ESO.PRO.REC1 keywords present in
  the input FITS header, and then gives the result to @ref
  muse_basicproc_params_new() to actually create the output structure.

  @note This only works, if the input header is that of a product produced by
        the MUSE pipeline.

  Use muse_basicproc_params_delete() to free memory allocated by this function.
 */
/*---------------------------------------------------------------------------*/
muse_basicproc_params *
muse_basicproc_params_new_from_propertylist(const cpl_propertylist *aHeader)
{
  cpl_ensure(aHeader, CPL_ERROR_NULL_INPUT, NULL);

  /* create a parameter list from the input header, so  *
   * that we can give it to muse_basicproc_params_new() */
  cpl_parameterlist *parlist = muse_cplparameterlist_from_propertylist(aHeader, 1);
  cpl_ensure(parlist, CPL_ERROR_ILLEGAL_INPUT, NULL);
  const char *recipe = cpl_propertylist_get_string(aHeader, "ESO PRO REC1 ID");
  char *prefix = cpl_sprintf("muse.%s", recipe);
  muse_basicproc_params *bpars = muse_basicproc_params_new(parlist, prefix);
  cpl_parameterlist_delete(parlist);
  cpl_free(prefix);
  return bpars;
} /* muse_basicproc_params_new_from_propertylist() */

/*---------------------------------------------------------------------------*/
/**
  @brief  Free a structure of basic processing parameters.
  @param  aBPars   the structure of basic processing parameters

  This deallocates the string and muse_image components and the structure
  itself.
 */
/*---------------------------------------------------------------------------*/
void
muse_basicproc_params_delete(muse_basicproc_params *aBPars)
{
  if (!aBPars) {
    return;
  }
  cpl_free(aBPars->overscan);
  cpl_free(aBPars->rejection);
  cpl_free(aBPars->crmethod);
  muse_image_delete(aBPars->flatimage);
  aBPars->flatimage = NULL;
  cpl_free(aBPars);
}

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Verify the chip/instrument setup of image and reference.
  @param  aImage   the image to compare
  @param  aRef     the reference to use
  @return CPL_ERROR_NONE on success or a CPL error code on failure or if
          there is a grave problem.

  @error{return CPL_ERROR_NULL_INPUT, aImage and/or aRef are NULL}
  @error{return CPL_ERROR_NULL_INPUT,
         header components of aImage and/or aRef are NULL}
  @error{output error message\, return CPL_ERROR_ILLEGAL_INPUT,
         aRef does not contain a PRO.CATG}
  @error{output error message\, return CPL_ERROR_TYPE_MISMATCH,
         any grave problem occurs when comparing aImage to aRef}
  @error{output warning message only,
         a not critical problem occurs when comparing aImage to aRef}
 */
/*---------------------------------------------------------------------------*/
static cpl_error_code
muse_basicproc_verify_setup(const muse_image *aImage, const muse_image *aRef)
{
  cpl_ensure_code(aImage && aRef, CPL_ERROR_NULL_INPUT);
  cpl_ensure_code(aImage->header && aRef->header, CPL_ERROR_NULL_INPUT);
  /* shortcuts to the headers */
  cpl_propertylist *him = aImage->header,
                   *href = aRef->header;
  /* reference image need to have a processed category (PRO.CATG) */
  const char *fn1 = cpl_propertylist_get_string(him, MUSE_HDR_TMP_FN),
             *fn2 = cpl_propertylist_get_string(href, MUSE_HDR_TMP_FN),
             *catg = muse_pfits_get_pro_catg(href);
  if (!catg) {
    cpl_msg_error(__func__, "\"%s\" does not contain a category (ESO.PRO.CATG)!",
                  fn2);
    return CPL_ERROR_ILLEGAL_INPUT;
  }

  cpl_boolean ignoreread = getenv("MUSE_DEBUG_IGNORE_READOUT")
                         && atoi(getenv("MUSE_DEBUG_IGNORE_READOUT")) > 0,
              ignoremode = getenv("MUSE_DEBUG_IGNORE_INSMODE")
                         && atoi(getenv("MUSE_DEBUG_IGNORE_INSMODE")) > 0;

  muse_ins_mode mode1 = muse_pfits_get_mode(him),
                mode2 = muse_pfits_get_mode(href);
  const char *modestr1 = muse_pfits_get_insmode(him),
             *modestr2 = muse_pfits_get_insmode(href),
             *rawtag = cpl_propertylist_get_string(him, MUSE_HDR_TMP_INTAG);
  int binx1 = muse_pfits_get_binx(him),
      biny1 = muse_pfits_get_biny(him),
      readid1 = muse_pfits_get_read_id(him),
      binx2 = muse_pfits_get_binx(href),
      biny2 = muse_pfits_get_biny(href),
      readid2 = muse_pfits_get_read_id(href);
  const char *readname1 = muse_pfits_get_read_name(him),
             *readname2 = muse_pfits_get_read_name(href),
             *chipname1 = muse_pfits_get_chip_name(him),
             *chipid1 = muse_pfits_get_chip_id(him),
             *chipname2 = muse_pfits_get_chip_name(href),
             *chipid2 = muse_pfits_get_chip_id(href);
  cpl_boolean chipinfo = chipname1 && chipid1
                       && chipname2 && chipid2;
  if (!chipinfo) {
    cpl_msg_warning(__func__, "CHIP information is missing (ESO.DET.CHIP."
                    "{NAME,ID}) from \"%s\" (%s, %s) or \"%s\" (%s, %s)",
                    fn1, chipname1, chipid1, fn2, chipname2, chipid2);
  }
  /* Everything should fail for non-matching binning. */
  if (binx1 != binx2 || biny1 != biny2) {
    cpl_msg_error(__func__, "Binning of \"%s\" (%dx%d) and \"%s\" (%dx%d) does "
                  "not match", fn1, binx1, biny1, fn2, binx2, biny2);
    return CPL_ERROR_TYPE_MISMATCH;
  }

  /* The pipeline should refuse to work if the bias is not taken in the   *
   * same read-out as the image that is being bias-subtracted. Hence it   *
   * should give an ERROR message when searching for bias files and stop. */
  if (!strncmp(catg, "MASTER_BIAS", 12)) {
    if (!ignoreread && (readid1 != readid2)) {
      cpl_msg_error(__func__, "Read-out mode of \"%s\" (%d: %s) and \"%s\" (%d:"
                    " %s) does not match", fn1, readid1, readname1, fn2,
                    readid2, readname2);
      return CPL_ERROR_TYPE_MISMATCH;
    }
    if (chipinfo && (strcmp(chipname1, chipname2) || strcmp(chipid1, chipid2))) {
      cpl_msg_error(__func__, "CHIP information (ESO.DET.CHIP.{NAME,ID}) "
                    "does not match for \"%s\" (%s, %s) and \"%s\" (%s, %s)",
                    fn1, chipname1, chipid1, fn2, chipname2, chipid2);
      return CPL_ERROR_TYPE_MISMATCH;
    }
  } /* if ref is bias */

  /* We probably need similar guards for the instrument mode, so that one  *
   * cannot use a flat-field in WFM-NOAO-N for data taken with WFM-NOAO-E. */
  if (!strncmp(catg, "MASTER_FLAT", 12)) {
    if (!ignoremode && (mode1 != mode2)) {
      if (rawtag && !strncmp(rawtag, MUSE_TAG_ILLUM, strlen(MUSE_TAG_ILLUM) + 1)) {
        /* ignore mode differences between ILLUM and other calibrations */
        cpl_msg_debug(__func__, "Instrument modes for \"%s\" (%s, is %s) and \"%s\""
                      " (%s) do not match", fn1, modestr1, rawtag, fn2, modestr2);
      } else {
        cpl_msg_error(__func__, "Instrument modes for \"%s\" (%s) and \"%s\" (%s)"
                      " do not match", fn1, modestr1, fn2, modestr2);
        return CPL_ERROR_TYPE_MISMATCH;
      } /* else */
    } /* if modes different */
  } /* if ref is flat */

  /* XXX add check to not mix WFM and NFM for illuminated exposures */

  /* It should output WARNINGs (but continue processing), when other *
   * calibrations are not in the same read-out mode or if the images *
   * originate from different chips or chip installation dates.      */
  if (!ignoreread && (readid1 != readid2)) {
    cpl_msg_warning(__func__, "Read-out mode of \"%s\" (%d: %s) and \"%s\" (%d:"
                    " %s) does not match", fn1, readid1, readname1, fn2,
                    readid2, readname2);
  }
  if (chipinfo && (strcmp(chipname1, chipname2) || strcmp(chipid1, chipid2))) {
    cpl_msg_warning(__func__, "CHIP information (ESO.DET.CHIP.{NAME,ID,DATE}) "
                    "does not match for \"%s\" (%s, %s) and \"%s\" (%s, %s)",
                    fn1, chipname1, chipid1, fn2, chipname2, chipid2);
  }

  return CPL_ERROR_NONE;
} /* muse_basicproc_verify_setup() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Compute overscan statistics and reject cosmic rays in the overscans.
  @param  aList    the image list
  @param  aBPars   basic processing parameters
  @return CPL_ERROR_NONE on success or a CPL error code on failure

  Use muse_quadrants_overscan_stats() on each image in aList to compute overscan
  statistics and save them into the header of that image.

  @error{set and return CPL_ERROR_NULL_INPUT, aList is NULL}
  @error{call muse_quadrants_overscan_stats() without rejection and with zero ignore,
         aBPars is NULL}
 */
/*---------------------------------------------------------------------------*/
static cpl_error_code
muse_basicproc_overscans_compute_stats(muse_imagelist *aList,
                                       muse_basicproc_params *aBPars)
{
  cpl_ensure_code(aList, CPL_ERROR_NULL_INPUT);
  unsigned int k;
  for (k = 0; k < aList->size; k++) {
    muse_image *image = muse_imagelist_get(aList, k);
    if (muse_quadrants_overscan_stats(image, aBPars ? aBPars->rejection : NULL,
                                      aBPars ? aBPars->ovscignore : 0)
        != CPL_ERROR_NONE) {
      cpl_msg_warning(__func__, "Could not compute overscan statistics in IFU "
                      "%hhu of exposure %u: %s", muse_utils_get_ifu(image->header),
                      k+1, cpl_error_get_message());
    } /* if */
  } /* for k (all images) */
  return CPL_ERROR_NONE;
} /* muse_basicproc_overscans_compute_stats() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Fit and subtract overscan polynomial model.
  @param  aList    the image list
  @param  aBPars   basic processing parameters
  @return CPL_ERROR_NONE on success or a CPL error code on failure

  @error{set and return CPL_ERROR_NULL_INPUT, aList is NULL}
  @error{return CPL_ERROR_NONE, aBPars->overscan does not start with "vpoly"}
  @error{propagate error from muse_quadrants_overscan_polyfit_vertical(),
         polynomial correction failed}
 */
/*---------------------------------------------------------------------------*/
static cpl_error_code
muse_basicproc_correct_overscans_vpoly(muse_imagelist *aList,
                                       muse_basicproc_params *aBPars)
{
  cpl_ensure_code(aList, CPL_ERROR_NULL_INPUT);
  cpl_boolean ovscvpoly = aBPars && aBPars->overscan
                        && !strncmp(aBPars->overscan, "vpoly", 5);
  if (!ovscvpoly) {
    cpl_msg_debug(__func__, "not vpoly: %s!", aBPars ? aBPars->overscan : "");
    return CPL_ERROR_NONE;
  }
  /* vertical polyfit requested, see if there are more parameters */
  unsigned char ovscvorder = 5;
  double frms = 1.01,
         fchisq = 1.04;
  char *rest = strchr(aBPars->overscan, ':');
  if (strlen(aBPars->overscan) > 6 && rest++) { /* try to access info after "vpoly:" */
    ovscvorder = strtol(rest, &rest, 10);
    if (strlen(rest++) > 0) { /* ++ to skip over the comma */
      frms = strtod(rest, &rest);
      if (strlen(rest++) > 0) {
        fchisq = strtod(rest, &rest);
      }
    }
  } /* if */
  cpl_msg_debug(__func__, "vpoly: %s (vorder=%hhu, frms=%f, fchisq=%f, ignore=%u,"
                " sigma=%.3f)", aBPars->overscan, ovscvorder, frms, fchisq,
                aBPars->ovscignore, aBPars->ovscsigma);

  cpl_error_code rc = CPL_ERROR_NONE;
  unsigned int k;
  for (k = 0; k < aList->size; k++) {
    muse_image *image = muse_imagelist_get(aList, k);
    rc = muse_quadrants_overscan_polyfit_vertical(image, aBPars->ovscignore,
                                                  ovscvorder, aBPars->ovscsigma,
                                                  frms, fchisq);
    if (rc != CPL_ERROR_NONE) {
      unsigned char ifu = muse_utils_get_ifu(image->header);
      cpl_msg_error(__func__, "Could not correct quadrants levels using vertical"
                    " overscan fit in IFU %hhu: %s", ifu, cpl_error_get_message());
    } /* if */
  } /* for k (all images) */
  return rc;
} /* muse_basicproc_correct_overscans_vpoly() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Trim prescan and overscan regions off images in a muse_image list.
  @param  aList    the image list
  @return CPL_ERROR_NONE on success or a CPL error code on failure

  This just loops through all images in the list and trims off the pre- and
  overscan regions using muse_quadrants_trim_image().

  @error{set and return CPL_ERROR_NULL_INPUT, aList is NULL}
  @error{propagate error from muse_quadrants_trim_image(),
         trimming of an image failed}
 */
/*---------------------------------------------------------------------------*/
static cpl_error_code
muse_basicproc_trim_images(muse_imagelist *aList)
{
  cpl_ensure_code(aList, CPL_ERROR_NULL_INPUT);

  unsigned int k;
  for (k = 0; k < aList->size; k++) {
    muse_image *image = muse_imagelist_get(aList, k);
    muse_image *trimmed = muse_quadrants_trim_image(image);
    cpl_ensure_code(trimmed, cpl_error_get_code());

    /* setting a new one deletes the old, so no need to free |image| */
    muse_imagelist_set(aList, trimmed, k);
  } /* for k (all images) */
  return muse_imagelist_is_uniform(aList) == 0 ? CPL_ERROR_NONE
                                               : CPL_ERROR_ILLEGAL_OUTPUT;
} /* muse_basicproc_trim_images() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Use overscan levels to offset bias data levels to reference bias.
  @param  aList         the image list
  @param  aProcessing   the processing structure
  @param  aBPars        basic processing parameters
  @return CPL_ERROR_NONE on success or a CPL error code on failure

  This is for lists of bias images where all images should be offset in level
  to the first bias in the list. For non-bias raw images, this function does
  nothing.

  @error{set and return CPL_ERROR_NULL_INPUT, aList or aProcessing are NULL}
  @error{output debug message\, return CPL_ERROR_NONE,
         intags of aProcessing is not MUSE_TAG_BIAS}
  @error{return CPL_ERROR_NONE, aBPars is NULL}
 */
/*---------------------------------------------------------------------------*/
static cpl_error_code
muse_basicproc_correct_overscans_offset(muse_imagelist *aList,
                                        muse_processing *aProcessing,
                                        muse_basicproc_params *aBPars)
{
  cpl_ensure_code(aList && aProcessing, CPL_ERROR_NULL_INPUT);
  /* this only makes sense to do for inputs that are bias */
  if (!muse_processing_check_intags(aProcessing, MUSE_TAG_BIAS, 5)) {
    return CPL_ERROR_NONE;
  }
  cpl_boolean ovscoffset = aBPars && aBPars->overscan
                         && !strncmp(aBPars->overscan, "offset", 7);
  if (!ovscoffset) {
    return CPL_ERROR_NONE; /* no correction necessary */
  }
  unsigned char ifu = muse_utils_get_ifu(muse_imagelist_get(aList, 0)->header);
  cpl_msg_info(__func__, "Running overscan correction using %u %s images in IFU"
               " %hhu", aList->size, MUSE_TAG_BIAS, ifu);
  muse_image *ref = muse_imagelist_get(aList, 0);
  unsigned int k;
  for (k = 1; k < aList->size; k++) {
    muse_quadrants_overscan_correct(ref, muse_imagelist_get(aList, k));
  } /* for k (all images except first) */
  return CPL_ERROR_NONE;
} /* muse_basicproc_correct_overscans_offset() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Check that overscans of all biases in the list have similar levels.
  @param  aList         the image list
  @param  aProcessing   the processing structure
  @param  aBPars        basic processing parameters
  @return CPL_ERROR_NONE on success or a CPL error code on failure;

  This function only does something, if overscan differences between exposures
  are not already treated in some way, i.e. only if aBPars->overscan == "none".

  Use the overscan mean+/-stdev of the first frame as reference to compare to in
  a loop over all other frames.  Store the final averaged value (and its
  external mean error computed from individual standard deviations and the
  deviation from the mean) in updated FITS header of the first input image, so
  that it can be propagated into the final combined master frame.

  @error{set and return CPL_ERROR_NULL_INPUT, aList or aProcessing are NULL}
  @error{return CPL_ERROR_NONE, intags of aProcessing is not MUSE_TAG_BIAS}
  @error{return CPL_ERROR_NONE, there are less than 2 images in aList}
  @error{use defaults (sigma = 1 and no overscan correction), aBPars is NULL}
  @error{set and return CPL_ERROR_INCOMPATIBLE_INPUT,
         overscan of images in the list do not match the overscans of the first image}
 */
/*---------------------------------------------------------------------------*/
static cpl_error_code
muse_basicproc_check_overscans(muse_imagelist *aList,
                               muse_processing *aProcessing,
                               muse_basicproc_params *aBPars)
{
  cpl_ensure_code(aList && aProcessing, CPL_ERROR_NULL_INPUT);
  /* files other than bias are already checked when  *
   * subtracting the bias, so we can skip this check */
  if (!muse_processing_check_intags(aProcessing, MUSE_TAG_BIAS, 5)) {
    return CPL_ERROR_NONE;
  }
  if (aList->size < 2) { /* doesn't make sense for single image lists */
    return CPL_ERROR_NONE;
  }
  cpl_boolean ovscnone = aBPars && aBPars->overscan
                       && !strncmp(aBPars->overscan, "none", 5);
  if (!ovscnone) { /* check not needed */
    return CPL_ERROR_NONE;
  }
  double sigma = aBPars ? aBPars->ovscsigma : 1.;

  /* header of the first image */
  muse_image *refimage = muse_imagelist_get(aList, 0);
  cpl_propertylist *refheader = refimage->header;
  cpl_error_code rc = CPL_ERROR_NONE;
  unsigned char n, ifu = muse_utils_get_ifu(refheader);
  for (n = 1; n <= 4; n++) {
    /* create correct header keyword names */
    char *keywordmean = cpl_sprintf(MUSE_HDR_OVSC_MEAN, n),
         *keywordstdev = cpl_sprintf(MUSE_HDR_OVSC_STDEV, n);

    /* overscan stats for first image in list */
    const char *reffn = cpl_propertylist_get_string(refheader, MUSE_HDR_TMP_FN);
    float refmean = cpl_propertylist_get_float(refheader, keywordmean),
          refstdev = cpl_propertylist_get_float(refheader, keywordstdev),
          hilimit = refmean + sigma * refstdev,
          lolimit = refmean - sigma * refstdev;
    /* variables to create average output values */
    double summean = refmean,
           sumstdev = pow(refstdev, 2.);

    /* compare with the other images */
    unsigned int k;
    for (k = 1; k < aList->size; k++) {
      cpl_propertylist *h = muse_imagelist_get(aList, k)->header;
      float mean = cpl_propertylist_get_float(h, keywordmean),
            stdev = cpl_propertylist_get_float(h, keywordstdev);
      summean += mean;
      sumstdev += pow(stdev, 2.) + pow(mean - refmean, 2.);
      const char *fn = cpl_propertylist_get_string(h, MUSE_HDR_TMP_FN);
      if (mean > hilimit || mean < lolimit) {
        cpl_msg_error(__func__, "Overscan of IFU %hhu, quadrant %1hhu of image %u [%s] "
                      "(%.3f+/-%.3f) differs from first image [%s] (%.3f+/-%.3f)!",
                      ifu, n, k+1, fn, mean, stdev, reffn, refmean, refstdev);
        rc = cpl_error_set(__func__, CPL_ERROR_INCOMPATIBLE_INPUT);
        continue; /* debug output only if there was no error */
      }
      cpl_msg_debug(__func__, "Overscan of IFU %hhu, quadrant %1hhu of image %u [%s] "
                    "%.3f+/-%.3f (first image [%s] %.3f+/-%.3f)",
                    ifu, n, k+1, fn, mean, stdev, reffn, refmean, refstdev);
    } /* for k (all images except first) */

    /* update values in the 1st header to be the averaged values *
     * which should be propagated to the final combined frame    */
    summean /= aList->size;
    sumstdev = sqrt(sumstdev) / aList->size;
    cpl_msg_debug(__func__, "Averaged overscan values in IFU %hhu, quadrant %1hhu: "
                  "%.3f +/- %.3f (%d images)", ifu, n, summean, sumstdev, aList->size);
    cpl_propertylist_update_float(refheader, keywordmean, summean);
    cpl_propertylist_update_float(refheader, keywordstdev, sumstdev);

    cpl_free(keywordmean);
    cpl_free(keywordstdev);
  } /* for n (quadrants) */

  return rc;
} /* muse_basicproc_check_overscans() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Apply external bad pixel table(s) to the muse_image list.
  @param  aList         the image list
  @param  aProcessing   the processing structure
  @param  aIFU          the IFU/channel number
  @return CPL_ERROR_NONE on success or a CPL error code on failure

  This just loops through all bad pixel tables in the input list and all images
  in the list, and ORs the values in the DQ extension of each image with what it
  finds in the badpix table.

  @error{set and return CPL_ERROR_NULL_INPUT, aList or aProcessing are NULL}
  @error{return CPL_ERROR_NONE,
         a table tagged with MUSE_TAG_BADPIX_TABLE is not passed in aProcessing}
  @error{propagate error from muse_cpltable_check(), bad pixel table has wrong format}
  @error{propagate error code, cpl_image_set() fails to set bad pixel}
 */
/*---------------------------------------------------------------------------*/
static cpl_error_code
muse_basicproc_apply_badpix(muse_imagelist *aList, muse_processing *aProcessing,
                            unsigned char aIFU)
{
  cpl_ensure_code(aList && aProcessing, CPL_ERROR_NULL_INPUT);
  cpl_frameset *frames = muse_frameset_find(aProcessing->inframes,
                                            MUSE_TAG_BADPIX_TABLE, aIFU,
                                            CPL_FALSE);
  if (!frames) { /* should not happen... */
    return CPL_ERROR_NONE;
  }

  /* loop through all frames that were found to have bad pixels for this IFU */
  cpl_errorstate prestate = cpl_errorstate_get();
  cpl_size iframe, nframes = cpl_frameset_get_size(frames);
  for (iframe = 0; iframe < nframes; iframe++) {
    cpl_frame *frame = cpl_frameset_get_position(frames, iframe);
    const char *fn = cpl_frame_get_filename(frame);
    char *extname = cpl_sprintf("CHAN%02hhu", aIFU);
    cpl_table *table = muse_cpltable_load(fn, extname, muse_badpix_table_def);
    cpl_free(extname);
    if (!table) {
      continue; /* file couldn't be loaded or is faulty */
    }
    muse_processing_append_used(aProcessing, frame, CPL_FRAME_GROUP_CALIB, 1);

    /* valid table found, loop through all entries and mark them in all images */
    int nrow = cpl_table_get_nrow(table);
    unsigned int k, nbadpix = 0, /* could transferred bad pixels */
                 nbadpos = 0; /* and entries with bad positions */
    for (k = 0; k < aList->size; k++) {
      muse_image *image = muse_imagelist_get(aList, k);
      /* XXX need to verify the detector/chip properties here, too! */
      int i;
      for (i = 0; i < nrow; i++) {
        int x = cpl_table_get(table, MUSE_BADPIX_X, i, NULL),
            y = cpl_table_get(table, MUSE_BADPIX_Y, i, NULL);

        /* get value from table and check that getting *
         * the value from the DQ image succeeds        */
        uint32_t dq = cpl_table_get(table, MUSE_BADPIX_DQ, i, NULL);
        cpl_errorstate state = cpl_errorstate_get();
        int err;
        uint32_t value = cpl_image_get(image->dq, x, y, &err);
        if (err != 0 || !cpl_errorstate_is_equal(state)) {
          cpl_errorstate_set(state); /* swallow the error, nothing was changed */
          if (k == 0) {
            nbadpos++;
          } /* if first image */
          continue;
        } /* if error occurred */

        /* OR the data in the DQ extension with what's in the badpix table */
        cpl_image_set(image->dq, x, y, dq | value);
        if (k == 0) {
          nbadpix++; /* count bad pixels transferred to 1st image */
        } /* if first image */
      } /* for i (all table rows) */
    } /* for k (all images) */
    cpl_table_delete(table);

    cpl_msg_debug(__func__, "Applied %u bad pixels from %s \"%s\" in IFU %hhu.",
                  nbadpix, MUSE_TAG_BADPIX_TABLE, fn, aIFU);
    /* warn, if there were bad entries */
    if (nbadpos > 0) {
      cpl_msg_warning(__func__, "\"%s\" contained %u entries outside the CCD in"
                      " IFU %hhu!", fn, nbadpos, aIFU);
    } /* if bad entries */
  } /* for iframe (all bad pixel tables found) */
  cpl_frameset_delete(frames);
  return cpl_errorstate_is_equal(prestate) ? CPL_ERROR_NONE
                                           : cpl_error_get_code();
} /* muse_basicproc_apply_badpix() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Mark all pixels above the saturation limit in the bad pixel mask.
  @param  aList   the image list
  @return CPL_ERROR_NONE on success or a CPL error code on failure

  @error{set and return CPL_ERROR_NULL_INPUT, aList is NULL}
 */
/*---------------------------------------------------------------------------*/
static cpl_error_code
muse_basicproc_check_saturation(muse_imagelist *aList)
{
  cpl_ensure_code(aList, CPL_ERROR_NULL_INPUT);
  unsigned int k;
  for (k = 0; k < aList->size; k++) {
    muse_image *image = muse_imagelist_get(aList, k);
    unsigned char ifu = muse_utils_get_ifu(image->header);
    int nsaturated = muse_quality_set_saturated(image);
    /* if we have more than 10% of saturated pixels then something is wrong */
    int npix = cpl_image_get_size_x(image->data)
             * cpl_image_get_size_y(image->data);
    if (nsaturated > (0.01 * npix)) {
      const char *fn = cpl_propertylist_get_string(image->header,
                                                   MUSE_HDR_TMP_FN);
      cpl_msg_error(__func__, "Raw exposure %u [%s] is strongly saturated in "
                    "IFU %hhu (%d of %d pixels)!", k+1, fn, ifu, nsaturated,
                    npix);
    } else if (nsaturated > (0.001 * npix)) {
      const char *fn = cpl_propertylist_get_string(image->header,
                                                   MUSE_HDR_TMP_FN);
      cpl_msg_warning(__func__, "Raw exposure %u [%s] is probably saturated in "
                      "IFU %hhu (%d of %d pixels)!", k+1, fn, ifu, nsaturated,
                      npix);
    } else {
      cpl_msg_debug(__func__, "Raw exposure %u in IFU %hhu (%d of %d pixels "
                    "saturated)!", k+1, ifu, nsaturated, npix);
    }
    cpl_propertylist_update_int(image->header, MUSE_HDR_TMP_NSAT, nsaturated);
  } /* for k (all images) */
  return CPL_ERROR_NONE;
} /* muse_basicproc_check_saturation() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Compute quadrant statistics in bias images.
  @param  aList         the list of bias images
  @param  aProcessing   the processing structure
  @return CPL_ERROR_NONE on success or a CPL error code on failure

  The "statistics" are simply the median value, they are recorded as temporary
  keywords (macro MUSE_HDR_TMP_QUADnMED) in the header of each image of aList.

  @error{set and return CPL_ERROR_NULL_INPUT, aList is NULL}
 */
/*---------------------------------------------------------------------------*/
static cpl_error_code
muse_basicproc_quadrant_statistics(muse_imagelist *aList,
                                   muse_processing *aProcessing)
{
  cpl_ensure_code(aList, CPL_ERROR_NULL_INPUT);
  /* this only makes sense to do for inputs that are bias */
  if (!muse_processing_check_intags(aProcessing, MUSE_TAG_BIAS, 5)) {
    return CPL_ERROR_NONE;
  }
  unsigned char ifu = muse_utils_get_ifu(muse_imagelist_get(aList, 0)->header);
  cpl_msg_info(__func__, "Computing per-quadrant medians for %u exposures in "
               "IFU %hhu", aList->size, ifu);
  unsigned int k;
  for (k = 0; k < aList->size; k++) {
    muse_image *image = muse_imagelist_get(aList, k);
    unsigned char n;
    for (n = 1; n <= 4; n++) {
      cpl_size *w = muse_quadrants_get_window(muse_imagelist_get(aList, k), n);
      float median = cpl_image_get_median_window(image->data, w[0], w[2],
                                                 w[1], w[3]);
      cpl_free(w);
      char *kw = cpl_sprintf(MUSE_HDR_TMP_QUADnMED, n);
      cpl_propertylist_append_float(image->header, kw, median);
      cpl_free(kw);
    } /* for n (quadrants) */
  } /* for k (all images) */
  return CPL_ERROR_NONE;
} /* muse_basicproc_quadrant_statistics() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Apply the bias to the muse_image list.
  @param  aList         the image list
  @param  aProcessing   the processing structure
  @param  aIFU          the IFU/channel number
  @param  aBPars        basic processing parameters
  @return CPL_ERROR_NONE on success or a CPL error code on failure

  In addition to subtracting the bias, this function also creates the variance
  images in the stat component of each input image.

  @error{set and return CPL_ERROR_NULL_INPUT, aList or aProcessing are NULL}
  @error{return CPL_ERROR_NONE, no bias found in inframes of aProcessing}
  @error{propagate error from muse_image_load()/muse_image_load_from_extensions(),
         bias image could not be loaded}
  @error{do not use overscan correction, aBPars is NULL}
  @error{set and return CPL_ERROR_INCOMPATIBLE_INPUT,
         overscan levels are 100 sigma off comparing bias against raw image}
  @error{propagate error from muse_image_subtract(), bias subtraction fails}
 */
/*---------------------------------------------------------------------------*/
static cpl_error_code
muse_basicproc_apply_bias(muse_imagelist *aList, muse_processing *aProcessing,
                          unsigned char aIFU, muse_basicproc_params *aBPars)
{
  cpl_ensure_code(aList && aProcessing, CPL_ERROR_NULL_INPUT);
  cpl_frame *biasframe = muse_frameset_find_master(aProcessing->inframes,
                                                   MUSE_TAG_MASTER_BIAS, aIFU);
  cpl_errorstate prestate = cpl_errorstate_get();
  if (!biasframe) return CPL_ERROR_NONE;
  cpl_errorstate_set(prestate);
  const char *biasname = cpl_frame_get_filename(biasframe);
  muse_image *biasimage = muse_image_load(biasname);
  if (!biasimage) {
    /* remember error message for a while, but reset the state, *
     * before trying to load from channel extensions            */
    char *errmsg = cpl_strdup(cpl_error_get_message());
    cpl_errorstate_set(prestate);
    biasimage = muse_image_load_from_extensions(biasname, aIFU);
    if (!biasimage) {
      /* now display both the older and the new error messages, *
       * so that they cannot be swallowed by parallelization    */
      cpl_msg_error(__func__, "%s", errmsg);
      cpl_msg_error(__func__, "%s", cpl_error_get_message());
      cpl_free(errmsg);
      cpl_frame_delete(biasframe);
      return cpl_error_get_code();
    } /* if 2nd load failure */
    cpl_free(errmsg);
    cpl_msg_info(__func__, "Bias correction in IFU %hhu using \"%s[CHAN%02hhu."
                 "DATA]\"", aIFU, biasname, aIFU);
  } else {
    cpl_msg_info(__func__, "Bias correction in IFU %hhu using \"%s[DATA]\"",
                 aIFU, biasname);
  }
  /* add temporary input filename for diagnostics */
  cpl_propertylist_append_string(biasimage->header, MUSE_HDR_TMP_FN, biasname);

  /* cross-check with the parameters of the master bias */
  muse_basicproc_params *mbpars
    = muse_basicproc_params_new_from_propertylist(biasimage->header);
  cpl_boolean parmatch = mbpars && aBPars; /* both sets exist */
  if (parmatch) {
    parmatch = parmatch && mbpars->overscan && aBPars->overscan
             && !strncmp(aBPars->overscan, mbpars->overscan,
                         CPL_MIN(strlen(aBPars->overscan), strlen(mbpars->overscan) + 1));
    parmatch = parmatch && mbpars->rejection && aBPars->rejection
             && !strncmp(aBPars->rejection, mbpars->rejection,
                         CPL_MIN(strlen(aBPars->rejection), strlen(mbpars->rejection) + 1));
    parmatch = parmatch && fabs(aBPars->ovscsigma - mbpars->ovscsigma)
                           < 10 * DBL_EPSILON;
    parmatch = parmatch && aBPars->ovscignore == mbpars->ovscignore;
  }
  if (!parmatch) {
    cpl_msg_warning(__func__, "overscan parameters differ between %s and recipe"
                    " parameters!", MUSE_TAG_MASTER_BIAS);
  } else {
    cpl_msg_debug(__func__, "overscan parameters between %s and given "
                  "parameters nicely match", MUSE_TAG_MASTER_BIAS);
  }
  muse_basicproc_params_delete(mbpars);

  cpl_boolean ovscoffset = aBPars && aBPars->overscan
                         && !strncmp(aBPars->overscan, "offset", 7),
              ovscvpoly = aBPars && aBPars->overscan
                        && !strncmp(aBPars->overscan, "vpoly", 5); /* up to the : */
  muse_processing_append_used(aProcessing, biasframe, CPL_FRAME_GROUP_CALIB, 0);
  cpl_error_code rc = CPL_ERROR_NONE;
  unsigned int k;
  for (k = 0; k < aList->size && rc == CPL_ERROR_NONE; k++) {
    muse_image *image = muse_imagelist_get(aList, k);
    rc = muse_basicproc_verify_setup(image, biasimage);
    if (rc != CPL_ERROR_NONE) {
      break;
    }
    rc = muse_image_variance_create(image, biasimage);
    if (ovscoffset) {
      muse_quadrants_overscan_correct(image, biasimage);
    } else if (ovscvpoly) {
      /* create error message and code, if the bias was not    *
       * vpoly-handled, i.e. still has the original bias level *
       * of ~1000 adu or > 100xsigma above ~zero               */
      cpl_boolean good = muse_quadrants_overscan_check(image, biasimage, 100.);
      if (!good) {
        cpl_msg_error(__func__, "Very different overscan levels found between "
                      "%s and raw exposure %u in IFU %hhu: incompatible overscan"
                      " parameters?", MUSE_TAG_MASTER_BIAS, k + 1, aIFU);
        rc = cpl_error_set(__func__, CPL_ERROR_INCOMPATIBLE_INPUT);
        break;
      } /* if not good */
    } else {
      /* just warn above indicated sigma level, ignore failure */
      muse_quadrants_overscan_check(image, biasimage,
                                    aBPars ? aBPars->ovscsigma : 1.);
    }
    rc = muse_image_subtract(image, biasimage);
#if GENERATE_TEST_IMAGES
    /* if we need to generate images for automated tests, here is a good *
     * place to write them out, as they are trimmed and bias corrected   */
    if (k == 0) {
      muse_image_save(image, "trimmed_bias_sub.fits");
    }
#endif
  } /* for k (all images) */
  muse_image_delete(biasimage);
  return rc;
} /* muse_basicproc_apply_bias() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Compute an estimate of the gain and cross-check values in the header.
  @param  aList         the image list
  @param  aProcessing   the processing structure
  @param  aIFU          the IFU/channel number
  @return CPL_ERROR_NONE on success or a CPL error code on failure

  This function should be called on data that was trimmed but not yet bias
  corrected nor converted to electrons. Due to the non-uniform illumination it
  probably gives an inaccurate estimate of the true gain.

  @error{set and return CPL_ERROR_NULL_INPUT, aList or aProcessing are NULL}
  @error{return CPL_ERROR_NONE, intags of aProcessing is not MUSE_TAG_FLAT}
  @error{return CPL_ERROR_INCOMPATIBLE_INPUT,
         there are less than 3 images in aList}
  @error{return CPL_ERROR_NONE, no bias found in inframes of aProcessing}
 */
/*---------------------------------------------------------------------------*/
static cpl_error_code
muse_basicproc_check_gain(muse_imagelist *aList, muse_processing *aProcessing,
                          unsigned char aIFU)
{
  cpl_ensure_code(aList && aProcessing, CPL_ERROR_NULL_INPUT);
  /* this only makes sense to do for input flat-fields */
  if (!muse_processing_check_intags(aProcessing, MUSE_TAG_FLAT, 5)) {
    return CPL_ERROR_NONE;
  }
  cpl_ensure_code(muse_imagelist_get_size(aList) >= 2,
                  CPL_ERROR_INCOMPATIBLE_INPUT);

  cpl_frame *fbias = muse_frameset_find_master(aProcessing->inframes,
                                               MUSE_TAG_MASTER_BIAS, aIFU);
  if (!fbias) {
    /* it's OK if there is no bias frame, probably we are not in muse_bias */
    return CPL_ERROR_NONE;
  }
  const char *biasname = cpl_frame_get_filename(fbias);
  cpl_propertylist *hbias = cpl_propertylist_load(biasname, 0);
  if (!cpl_propertylist_has(hbias, QC_BIAS_MASTER_RON)) {
    cpl_propertylist_delete(hbias);
    /* try to load from channel extension */
    char *extname = cpl_sprintf("CHAN%02hhu.%s", aIFU, EXTNAME_DATA);
    int extension = cpl_fits_find_extension(biasname, extname);
    hbias = cpl_propertylist_load(biasname, extension);
    cpl_free(extname);
  }
  cpl_frame_delete(fbias);
  cpl_image *f1 = muse_imagelist_get(aList, 0)->data,
            *f2 = muse_imagelist_get(aList, 1)->data;
  cpl_propertylist *header = muse_imagelist_get(aList, 0)->header;
  cpl_image *diff = cpl_image_subtract_create(f1, f2);

  unsigned char n;
  for (n = 1; n <= 4; n++) {
    cpl_size *w = muse_quadrants_get_window(muse_imagelist_get(aList, 0), n);
    double m1 = cpl_image_get_mean_window(f1, w[0], w[2], w[1], w[3]),
           m2 = cpl_image_get_mean_window(f2, w[0], w[2], w[1], w[3]),
           sf = cpl_image_get_stdev_window(diff, w[0], w[2], w[1], w[3]);
    char *keyword = cpl_sprintf(QC_BIAS_MASTERn_PREFIX" MEAN", n);
    float mb = cpl_propertylist_get_float(hbias, keyword);
    cpl_free(keyword);
    keyword = cpl_sprintf(QC_BIAS_MASTER_RON, n);
    /* the RON formula taken from Howell inverted to give sigma(b1-b2) */
    double gainheader = muse_pfits_get_gain(header, n),
           sb = cpl_propertylist_get_float(hbias, keyword) * sqrt(2.)
              / gainheader, /* gain in count/adu */
           /* the gain formula taken from Howell: */
           gain = (m1 + m2 - 2*mb) / (sf*sf - sb*sb),
           dgain = 1. - gain / gainheader;
    /* output warning for difference larger than 20%, info message otherwise */
    if (dgain > 0.2) {
      cpl_msg_warning(__func__, "IFU %hhu, quadrant %1hhu: estimated gain %.3f "
                      "count/adu but header gives %.3f!", aIFU, n, gain,
                      gainheader);
    } else {
      cpl_msg_info(__func__, "IFU %hhu, quadrant %1hhu: estimated gain %.3f "
                   "count/adu (header %.3f, delta %.3f)", aIFU, n, gain,
                   gainheader, dgain);
    }
    cpl_free(keyword);
    cpl_free(w);
  } /* for n (quadrants) */

  cpl_image_delete(diff);
  cpl_propertylist_delete(hbias);

  return CPL_ERROR_NONE;
} /* muse_basicproc_check_gain() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Override the gain in the raw header from that of a separate file.
  @param  aList         the image list
  @param  aProcessing   the processing structure
  @param  aIFU          the IFU/channel number
  @return CPL_ERROR_NONE on success or a CPL error code on failure

  This function transfers the gain values from the headers of the file tagged
  MUSE_TAG_NONLINGAIN.

  If the environment variable MUSE_BASICPROC_SKIP_GAIN_OVERRIDE is set to a
  positive integer, this function returns without overriding the gain.

  @error{set and return CPL_ERROR_NULL_INPUT, aList or aProcessing are NULL}
  @error{return CPL_ERROR_NONE,
         no nonlinearity file found in inframes of aProcessing}
 */
/*---------------------------------------------------------------------------*/
static cpl_error_code
muse_basicproc_gain_override(muse_imagelist *aList,
                             muse_processing *aProcessing, unsigned char aIFU)
{
  cpl_ensure_code(aList && aProcessing, CPL_ERROR_NULL_INPUT);

  if (muse_processing_check_intags(aProcessing, MUSE_TAG_LINEARITY_BIAS, 16) ||
      muse_processing_check_intags(aProcessing, MUSE_TAG_LINEARITY_FLAT, 15))
  {
    return CPL_ERROR_NONE;
  }

  cpl_frame *fnonlingain = muse_frameset_find_master(aProcessing->inframes,
                                                     MUSE_TAG_NONLINGAIN, aIFU);
  if (!fnonlingain) {
    /* it's OK if there is no nonlinearity frame, it's optional... */
    return CPL_ERROR_NONE;
  }
  if (getenv("MUSE_BASICPROC_SKIP_GAIN_OVERRIDE") &&
      atoi(getenv("MUSE_BASICPROC_SKIP_GAIN_OVERRIDE")) > 0) {
    cpl_msg_info(__func__, "Skipping gain override, although %s is given",
                 MUSE_TAG_NONLINGAIN);
    return CPL_ERROR_NONE;
  }
  cpl_errorstate state = cpl_errorstate_get();
  const char *fn = cpl_frame_get_filename(fnonlingain);
  /* immediately try to load from channel extension */
  char *extname = cpl_sprintf("CHAN%02hhu", aIFU);
  int extension = cpl_fits_find_extension(fn, extname);
  cpl_propertylist *header = cpl_propertylist_load(fn, extension);
  cpl_msg_info(__func__, "Overriding gain in IFU %hhu using \"%s[%s]\"",
               aIFU, fn, extname);
  cpl_free(extname);
  muse_processing_append_used(aProcessing, fnonlingain, CPL_FRAME_GROUP_CALIB, 0);

  unsigned int k;
  for (k = 0; k < aList->size; k++) {
    muse_image *image = muse_imagelist_get(aList, k);
    /* XXX need to verify the detector/chip properties here, too! */
    unsigned char n;
    for (n = 1; n <= 4; n++) {
      /* transfer the GAIN of this quadrant into the image */
      double gain = muse_pfits_get_gain(header, n);
      char *kw = cpl_sprintf("ESO DET OUT%d GAIN", n);
      cpl_propertylist_update_double(image->header, kw, gain);
      cpl_free(kw);
    } /* for n (quadrants) */
  } /* for k (all images) */
  cpl_propertylist_delete(header);
  return cpl_errorstate_is_equal(state) ? CPL_ERROR_NONE : cpl_error_get_code();
} /* muse_basicproc_gain_override() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Use the gain to convert images in the list from counts to electrons.
  @param  aList         the image list
  @param  aProcessing   the processing structure
  @return CPL_ERROR_NONE on success or a CPL error code on failure

  This function should be called on data that was bias corrected.
  It will give bad results if the GAIN values in the FITS headers are wrong.

  @error{set and return CPL_ERROR_NULL_INPUT, aList or aProcessing are NULL}
  @error{return CPL_ERROR_NONE, intags of aProcessing is not MUSE_TAG_BIAS}
 */
/*---------------------------------------------------------------------------*/
static cpl_error_code
muse_basicproc_adu_to_count(muse_imagelist *aList, muse_processing *aProcessing)
{
  cpl_ensure_code(aList && aProcessing, CPL_ERROR_NULL_INPUT);
  /* this only makes sense to do for inputs that are not bias, and not *
   * inputs used to measure gain and linearity                         */
  if (muse_processing_check_intags(aProcessing, MUSE_TAG_BIAS, 5) ||
      muse_processing_check_intags(aProcessing, MUSE_TAG_LINEARITY_BIAS, 16) ||
      muse_processing_check_intags(aProcessing, MUSE_TAG_LINEARITY_FLAT, 15)) {
    return CPL_ERROR_NONE;
  }

  unsigned char ifu = muse_utils_get_ifu(muse_imagelist_get(aList, 0)->header);
  cpl_msg_info(__func__, "Converting %u exposures from adu to count (= electron)"
               " units in IFU %hhu", aList->size, ifu);
  cpl_error_code rc = CPL_ERROR_NONE;
  unsigned int k;
  for (k = 0; k < aList->size; k++) {
    rc = muse_image_adu_to_count(muse_imagelist_get(aList, k));
  } /* for k (all images) */
  return rc;
} /* muse_basicproc_adu_to_count() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Correct nonlinearity (of illuminated exposures).
  @param  aList         the image list
  @param  aProcessing   the processing structure
  @param  aIFU          the IFU/channel number
  @return CPL_ERROR_NONE on success or a CPL error code on failure

  This function uses the per-quadrant polynomials from the headers of the
  file tagged MUSE_TAG_NONLINGAIN to correct data and variance.

  For files tagged MUSE_TAG_BIAS or MUSE_TAG_DARK this function does not do
  anything.

  If the environment variable MUSE_BASICPROC_SKIP_NONLIN_CORR is set to a
  positive integer, this function returns without applying the correction.

  @error{set and return CPL_ERROR_NULL_INPUT, aList or aProcessing are NULL}
  @error{return CPL_ERROR_NONE,
         no nonlinearity file found in inframes of aProcessing}
 */
/*---------------------------------------------------------------------------*/
static cpl_error_code
muse_basicproc_corr_nonlinearity(muse_imagelist *aList,
                                 muse_processing *aProcessing,
                                 unsigned char aIFU)
{
  cpl_ensure_code(aList && aProcessing, CPL_ERROR_NULL_INPUT);
  /* this only makes sense to do for inputs that are illuminated and  *
   * inputs which are not used to measure gain and linearity          */
  if (muse_processing_check_intags(aProcessing, MUSE_TAG_BIAS, 5) ||
      muse_processing_check_intags(aProcessing, MUSE_TAG_DARK, 5) ||
      muse_processing_check_intags(aProcessing, MUSE_TAG_LINEARITY_BIAS, 16) ||
      muse_processing_check_intags(aProcessing, MUSE_TAG_LINEARITY_FLAT, 15)) {
    return CPL_ERROR_NONE;
  }

  cpl_frame *fnonlingain = muse_frameset_find_master(aProcessing->inframes,
                                                     MUSE_TAG_NONLINGAIN, aIFU);
  if (!fnonlingain) {
    /* it's OK if there is no nonlinearity frame, it's optional... */
    return CPL_ERROR_NONE;
  }
  if (getenv("MUSE_BASICPROC_SKIP_NONLIN_CORR") &&
      atoi(getenv("MUSE_BASICPROC_SKIP_NONLIN_CORR")) > 0) {
    cpl_msg_debug(__func__, "Skipping nonlinearity correction, although %s is "
                  "given", MUSE_TAG_NONLINGAIN);
    return CPL_ERROR_NONE;
  }
  const char *fn = cpl_frame_get_filename(fnonlingain);
  /* immediately try to load from channel extension */
  char *extname = cpl_sprintf("CHAN%02hhu", aIFU);
  int extension = cpl_fits_find_extension(fn, extname);
  cpl_propertylist *header = cpl_propertylist_load(fn, extension);
  cpl_msg_info(__func__, "Correcting nonlinearity in IFU %hhu using \"%s[%s]\"",
               aIFU, fn, extname);
  cpl_free(extname);
  muse_processing_append_used(aProcessing, fnonlingain, CPL_FRAME_GROUP_CALIB, 0);

  cpl_error_code rc = CPL_ERROR_NONE;
  unsigned int k;
  for (k = 0; k < aList->size; k++) {
    muse_image *image = muse_imagelist_get(aList, k);
    int nx = cpl_image_get_size_x(image->data);
    float *data = cpl_image_get_data_float(image->data),
          *stat = cpl_image_get_data_float(image->stat);

    /* XXX need to verify the detector/chip properties here, too! */

    unsigned char n;
    for (n = 1; n <= 4; n++) {
      /* create the 1D nonlinearity polynomial for this quadrant, then read   *
       * all the parameters from the header set the polynomial up accordingly */
      cpl_polynomial *poly = cpl_polynomial_new(1);
      char *kw = cpl_sprintf(MUSE_HDR_NONLINn_ORDER, n);
      unsigned char o, order = cpl_propertylist_get_int(header, kw);
      cpl_free(kw);
      for (o = 0; o <= order; o++) {
        kw = cpl_sprintf(MUSE_HDR_NONLINn_COEFFo, n, o);
        cpl_size pows = o;
        cpl_polynomial_set_coeff(poly, &pows,
                                 cpl_propertylist_get_double(header, kw));
        cpl_free(kw);
      } /* for i (polynomial orders) */
      /* compute linear extrapolations */
      kw = cpl_sprintf(MUSE_HDR_NONLINn_LLO, n);
      double lolim = cpl_propertylist_get_double(header, kw);
      cpl_free(kw);
      kw = cpl_sprintf(MUSE_HDR_NONLINn_LHI, n);
      double hilim = cpl_propertylist_get_double(header, kw);
      cpl_free(kw);
      /* coefficients for linear extrapolation beyond low and high limits */
      double p1lo, p0lo = cpl_polynomial_eval_1d(poly, lolim, &p1lo),
             p1hi, p0hi = cpl_polynomial_eval_1d(poly, hilim, &p1hi);
      /* convert limits from log10(adu) to adu */
      lolim = pow(10, lolim);
      hilim = pow(10, hilim);
#if 0
      double values[] = { 1., 20., 200., 2000., 20000., 65000.,
                          lolim, hilim, -1. };
      int idx = 0;
      while (values[idx] > 0) {
        double v = values[idx],
               logv = log10(v);
        cpl_msg_debug(__func__, "%f adu -> %f log10(adu) ==> poly = %f ==> x %f",
                      v, logv, cpl_polynomial_eval_1d(poly, logv, NULL),
                      1. / (1. + cpl_polynomial_eval_1d(poly, logv, NULL)));
        idx++;
      } /* while */
      cpl_polynomial_dump(poly, stdout);
      fflush(stdout);
      cpl_msg_debug(__func__, "beyond low limit (%f adu):  %f + (%f) * (log10(adu) - %f)",
                    lolim, p0lo, p1lo, log10(lolim));
      cpl_msg_debug(__func__, "beyond high limit (%f adu): %f + (%f) * (log10(adu) - %f)",
                    hilim, p0hi, p1hi, log10(hilim));
#endif

      /* now loop through the full quadrant and correct the scaling */
      cpl_size *w = muse_quadrants_get_window(image, n);
      int i;
      for (i = w[0] - 1; i < w[1]; i++) {
        int j;
        for (j = w[2] - 1; j < w[3]; j++) {
          if (data[i + j*nx] <= 0) {
            continue; /* don't do anything for non-positive datapoints */
          }
          /* compute the percent-level deviation (applying to log10(adu)) */
          double pcor;
          if (data[i + j*nx] < lolim) {
            pcor = p0lo + p1lo * (log10(data[i + j*nx]) - log10(lolim));
          } else if (data[i + j*nx] > hilim) {
            pcor = p0hi + p1hi * (log10(data[i + j*nx]) - log10(hilim));
          } else {
            pcor = cpl_polynomial_eval_1d(poly, log10(data[i + j*nx]), NULL);
          }
          /* then derive the multiplicative correction factor */
          double fcor = 1. / (1. + pcor);
          data[i + j*nx] *= fcor;
          stat[i + j*nx] *= fcor*fcor;
        } /* for j (vertical pixels) */
      } /* for i (horizontal pixels) */
      cpl_free(w);
      cpl_polynomial_delete(poly);
    } /* for n (quadrants) */
  } /* for k (all images) */
  cpl_propertylist_delete(header);

  return rc;
} /* muse_basicproc_corr_nonlinearity() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Apply the dark to the muse_image list.
  @param  aList         the image list
  @param  aProcessing   the processing structure
  @param  aIFU          the IFU/channel number
  @return CPL_ERROR_NONE on success or a CPL error code on failure

  @error{set and return CPL_ERROR_NULL_INPUT, aList or aProcessing are NULL}
  @error{return CPL_ERROR_NONE, no dark found in inframes of aProcessing}
  @error{propagate error from muse_image_load()/muse_image_load_from_extensions(),
         dark image could not be loaded}
  @error{propagate error from muse_image_scale() and/or muse_image_subtract(),
         dark scaling or subtraction fails}
 */
/*---------------------------------------------------------------------------*/
static cpl_error_code
muse_basicproc_apply_dark(muse_imagelist *aList, muse_processing *aProcessing,
                          unsigned char aIFU)
{
  cpl_ensure_code(aList && aProcessing, CPL_ERROR_NULL_INPUT);

  cpl_frame *darkframe = muse_frameset_find_master(aProcessing->inframes,
                                                   MUSE_TAG_MASTER_DARK, aIFU);
  cpl_errorstate prestate = cpl_errorstate_get();
  if (!darkframe) return CPL_ERROR_NONE;
  cpl_errorstate_set(prestate);
  const char *darkname = cpl_frame_get_filename(darkframe);
  muse_image *darkimage = muse_image_load(darkname);
  if (!darkimage) {
    char *errmsg = cpl_strdup(cpl_error_get_message());
    cpl_errorstate_set(prestate);
    darkimage = muse_image_load_from_extensions(darkname, aIFU);
    if (!darkimage) {
      cpl_msg_error(__func__, "%s", errmsg);
      cpl_msg_error(__func__, "%s", cpl_error_get_message());
      cpl_free(errmsg);
      cpl_frame_delete(darkframe);
      return cpl_error_get_code();
    } /* if 2nd load failure */
    cpl_free(errmsg);
    cpl_msg_info(__func__, "Dark correction in IFU %hhu using \"%s[CHAN%02hhu."
                 "DATA]\"", aIFU, darkname, aIFU);
  } else {
    cpl_msg_info(__func__, "Dark correction in IFU %hhu using \"%s[DATA]\"",
                 aIFU, darkname);
  }
  cpl_propertylist_append_string(darkimage->header, MUSE_HDR_TMP_FN, darkname);

  muse_processing_append_used(aProcessing, darkframe, CPL_FRAME_GROUP_CALIB, 0);
  cpl_error_code rc = CPL_ERROR_NONE;
  unsigned int k;
  for (k = 0; k < aList->size; k++) {
    muse_image *image = muse_imagelist_get(aList, k);
    rc = muse_basicproc_verify_setup(image, darkimage);
    if (rc != CPL_ERROR_NONE) {
      break;
    }

    /* duplicate the dark, because the scaling will destroy its normalization */
    muse_image *dark = muse_image_duplicate(darkimage);
    /* scale and subtract the dark by comparing exposure time *
     * of the dark and the other image                        */
    double scale = muse_pfits_get_exptime(image->header);
    if (muse_pfits_get_exptime(dark->header) > 0) {
      scale /= muse_pfits_get_exptime(dark->header);
    }
    rc = muse_image_scale(dark, scale);
    rc = muse_image_subtract(image, dark);

    muse_image_delete(dark);
  } /* for k (all images) */
  muse_image_delete(darkimage);

  return rc;
} /* muse_basicproc_apply_dark() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Apply cosmic ray finder to bias/dark corrected image.
  @param  aList    the image list
  @param  aBPars   basic processing parameters
  @return CPL_ERROR_NONE on success or a CPL error code on failure

  Only the DCR method is currently supported.

  @error{set and return CPL_ERROR_NULL_INPUT, aList is NULL}
  @error{return CPL_ERROR_NONE,
         aBPars is NULL or the crmethod of aBPars is not "dcr"}
  @error{propagate error from muse_cosmics_dcr(), cosmic ray rejection fails}
 */
/*---------------------------------------------------------------------------*/
static cpl_error_code
muse_basicproc_apply_cr(muse_imagelist *aList,
                        muse_basicproc_params *aBPars)
{
  cpl_ensure_code(aList, CPL_ERROR_NULL_INPUT);

  /* if we didn't get cr parameters we were not supposed *
   * to find cosmic rays, so just return without error   */
  cpl_boolean isdcr = aBPars && aBPars->crmethod
                    && !strncmp(aBPars->crmethod, "dcr", 4);
  if (!isdcr) {
    return CPL_ERROR_NONE;
  }

  cpl_error_code rc = CPL_ERROR_NONE;
  unsigned int k;
  for (k = 0; k < aList->size; k++) {
    muse_image *image = muse_imagelist_get(aList, k);
    unsigned char ifu = muse_utils_get_ifu(image->header);
    int ncr = 0;
    ncr = muse_cosmics_dcr(image, aBPars->dcrxbox, aBPars->dcrybox,
                           aBPars->dcrpasses, aBPars->dcrthres);
    if (ncr <= 0) {
      cpl_msg_error(__func__, "Cosmic ray rejection using DCR in IFU %hhu failed"
                    " for image %u (ncr = %d): %s", ifu, k+1, ncr,
                    cpl_error_get_message());
      rc = cpl_error_get_code();
    } else {
      cpl_msg_info(__func__, "Cosmic ray rejection using DCR in IFU %hhu found "
                   "%d affected pixels in image %u", ifu, ncr, k+1);
    }
  } /* for k (all images) */

  return rc;
} /* muse_basicproc_apply_cr() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Apply the flat to the muse_image list.
  @param  aList         the image list
  @param  aProcessing   the processing structure
  @param  aBPars        basic processing parameters
  @param  aIFU          the IFU/channel number
  @return CPL_ERROR_NONE on success or a CPL error code on failure

  @error{set and return CPL_ERROR_NULL_INPUT, aList or aProcessing are NULL}
  @error{return CPL_ERROR_NONE, no flat found in inframes of aProcessing}
  @error{propagate error from muse_image_load()/muse_image_load_from_extensions(),
         flat image could not be loaded}
  @error{propagate error from muse_image_divide(), flat division fails}
 */
/*---------------------------------------------------------------------------*/
static cpl_error_code
muse_basicproc_apply_flat(muse_imagelist *aList, muse_processing *aProcessing,
                          muse_basicproc_params *aBPars, unsigned char aIFU)
{
  cpl_ensure_code(aList && aProcessing, CPL_ERROR_NULL_INPUT);

  cpl_frame *flatframe = muse_frameset_find_master(aProcessing->inframes,
                                                   MUSE_TAG_MASTER_FLAT, aIFU);
  cpl_errorstate prestate = cpl_errorstate_get();
  if (!flatframe) return CPL_ERROR_NONE;
  cpl_errorstate_set(prestate);
  const char *flatname = cpl_frame_get_filename(flatframe);
  prestate = cpl_errorstate_get();
  muse_image *flatimage = muse_image_load(flatname);
  if (!flatimage) {
    char *errmsg = cpl_strdup(cpl_error_get_message());
    cpl_errorstate_set(prestate);
    flatimage = muse_image_load_from_extensions(flatname, aIFU);
    if (!flatimage) {
      cpl_msg_error(__func__, "%s", errmsg);
      cpl_msg_error(__func__, "%s", cpl_error_get_message());
      cpl_free(errmsg);
      cpl_frame_delete(flatframe);
      return cpl_error_get_code();
    } /* if 2nd load failure */
    cpl_free(errmsg);
    cpl_msg_info(__func__, "Flat-field correction in IFU %hhu using \"%s[CHAN%02hhu."
                 "DATA]\"", aIFU, flatname, aIFU);
  } else {
    cpl_msg_info(__func__, "Flat-field correction in IFU %hhu using \"%s[DATA]\"",
                 aIFU, flatname);
  }
  /* copy QC parameter containing integrated flux */
  double fflux = cpl_propertylist_get_double(flatimage->header,
                                             QC_FLAT_MASTER_INTFLUX);
  cpl_propertylist_append_string(flatimage->header, MUSE_HDR_TMP_FN, flatname);
  muse_processing_append_used(aProcessing, flatframe, CPL_FRAME_GROUP_CALIB, 0);

  cpl_error_code rc = CPL_ERROR_NONE;
  unsigned int k;
  for (k = 0; k < aList->size; k++) {
    cpl_msg_debug(__func__, "Flat-field division in IFU %hhu of image %u of %u",
                  aIFU, k+1, aList->size);
    muse_image *image = muse_imagelist_get(aList, k);
    rc = muse_basicproc_verify_setup(image, flatimage);
    if (rc != CPL_ERROR_NONE) {
      break;
    }
    rc = muse_image_divide(image, flatimage);
    cpl_propertylist_update_double(image->header, MUSE_HDR_FLAT_FLUX_LAMP, fflux);
  } /* for k (all images) */

  if (aBPars && aBPars->keepflat) {
    aBPars->flatimage = flatimage;
  } else {
    muse_image_delete(flatimage);
  }
  return rc;
} /* muse_basicproc_apply_flat() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  read the raw input files.
  @param  aProcessing   the processing structure
  @param  aIFU          the IFU/channel number
  @return the muse_image list

  Read the raw input files and store them as images in the images
  struct.  The files are taken from the inframes member of the
  processing struct.  On success, the used frames are added to the
  usedframes member of this struct.

  @error{set CPL_ERROR_NULL_INPUT\, return NULL, aProcessing is NULL}
  @error{set CPL_ERROR_NULL_INPUT\, return NULL,
         there were no raw frames matching the processing tag}
  @error{free image list\, set CPL_ERROR_ILLEGAL_OUTPUT\, return NULL,
         generated image list was empty or non uniform}
  @error{free image list\, propagate MUSE_ERROR_CHIP_NOT_LIVE\, return NULL,
         generated image list was empty but chip was known to be dead}
 */
/*---------------------------------------------------------------------------*/
static muse_imagelist *
muse_basicproc_load_raw(muse_processing *aProcessing, unsigned char aIFU)
{
  cpl_ensure(aProcessing, CPL_ERROR_NULL_INPUT, NULL);

  /* rawframes contains the list of input files that we want to load */
  cpl_frameset *rawframes = muse_frameset_check_raw(aProcessing->inframes,
                                                    aProcessing->intags, aIFU);
  cpl_ensure(rawframes, CPL_ERROR_NULL_INPUT, NULL);

  muse_imagelist *images = muse_imagelist_new();
  unsigned int k = 0;
  cpl_size iframe, nframes = cpl_frameset_get_size(rawframes);
  for (iframe = 0; iframe < nframes; iframe++) {
    cpl_frame *frame = cpl_frameset_get_position(rawframes, iframe);
    const char *fileName = cpl_frame_get_filename(frame);
    int extension = muse_utils_get_extension_for_ifu(fileName, aIFU);
    if (extension == -1) {
      cpl_msg_error(__func__, "\"%s\" does not contain data from IFU %hhu",
                    fileName, aIFU);
      break;
    }
    muse_image *raw = muse_image_load_from_raw(fileName, extension);
    if (!raw) {
      continue;
    }
    /* add temporary input filename for diagnostics */
    cpl_propertylist_append_string(raw->header, MUSE_HDR_TMP_FN, fileName);
    /* add temporary input tag, for use when deriving the output tag */
    cpl_propertylist_append_string(raw->header, MUSE_HDR_TMP_INTAG,
                                   cpl_frame_get_tag(frame));

    muse_imagelist_set(images, raw, k);
    muse_processing_append_used(aProcessing, frame, CPL_FRAME_GROUP_RAW, 1);
    k++;
  } /* for frame (all rawframes) */
  cpl_frameset_delete(rawframes);

  /* some checks if we ended up with a valid imagelist */
  if (!images->size) {
    if ((int)cpl_error_get_code() != MUSE_ERROR_CHIP_NOT_LIVE) {
      cpl_error_set(__func__, CPL_ERROR_ILLEGAL_OUTPUT);
      cpl_msg_error(__func__, "No raw images loaded for IFU %hhu", aIFU);
    }
    muse_imagelist_delete(images); /* still free the structure */
    return NULL;
  }
  if (muse_imagelist_is_uniform(images) != 0) {
    cpl_error_set(__func__, CPL_ERROR_ILLEGAL_OUTPUT);
    cpl_msg_error(__func__, "Non-uniform imagelist for IFU %hhu", aIFU);
    muse_imagelist_delete(images);
    return NULL;
  }

  return images;
} /* muse_basicproc_load_raw() */

/*---------------------------------------------------------------------------*/
/**
  @brief  Load the raw input files from disk and do basic processing.
  @param  aProcessing   the processing structure
  @param  aIFU          the IFU/channel number
  @param  aBPars        basic processing parameters
  @return the muse_imagelist * or NULL on error

  Read the raw input files from disk and store them as images in the images
  structure. Apply bad pixel table, bias, dark, and flat files to them,
  thereby adding two images with bad pixel information and variance to each
  of them, so that they are all of format muse_image. Saturated pixels are
  added to the bad pixel extension of every image.

  For raw bias exposures, the overscan levels are checked against the first
  input file. A mean overscan level is then added to the header of the first
  image to be propagated into the final combined master bias.

  For raw flat-field exposures, the gain is estimated from the data and compared
  to the one listed in the input FITS header.

  The necessary calibrations for each step are loaded as required.

  @error{propagate error from muse_processing_check_input()\, return NULL,
         input data is invalid somehow (NULL aProcessing pointer etc.)}
  @error{set an error code\, return NULL, raw images could not be loaded}
  @error{reset errors and continue, applying the bad pixel table failed}
  @error{set a cpl_error_code\, return NULL, saturation check failed}
  @error{set a cpl_error_code\, return NULL,
         correcting quadrants by overscan fitting failed}
  @error{set a cpl_error_code\, return NULL,
         computing overscan statistics and rejecting cosmic rays in overscans failed}
  @error{set a cpl_error_code\, return NULL, images could not be trimmed}
  @error{output warning and continue, overscan correction was requested but failed}
  @error{set a cpl_error_code\, return NULL, overscan check failed}
  @error{set a cpl_error_code\, return NULL, applying bias failed}
  @error{ignore errors and continue, checking gain failed}
  @error{set a cpl_error_code\, return NULL, coversion from adu to count failed}
  @error{set a cpl_error_code\, return NULL, applying dark failed}
  @error{set a cpl_error_code\, return NULL, carrying out cr rejection failed}
  @error{set a cpl_error_code\, return NULL, applying simple sky subtraction failed}
  @error{set a cpl_error_code\, return NULL, applying flat failed}
 */
/*---------------------------------------------------------------------------*/
muse_imagelist *
muse_basicproc_load(muse_processing *aProcessing, unsigned char aIFU,
                    muse_basicproc_params *aBPars)
{
  if (muse_processing_check_input(aProcessing, aIFU) != CPL_ERROR_NONE) {
    return NULL;
  }
  muse_imagelist *images = muse_basicproc_load_raw(aProcessing, aIFU);
  if (!images) {
    return NULL;
  }
  cpl_errorstate prestate = cpl_errorstate_get();
  if (muse_basicproc_apply_badpix(images, aProcessing, aIFU) != CPL_ERROR_NONE) {
    cpl_msg_warning(__func__, "Applying bad pixel table to IFU %hhu failed: %s",
                    aIFU, cpl_error_get_message());
    cpl_errorstate_set(prestate); /* this is not critical, continue */
  }
  if (muse_basicproc_check_saturation(images) != CPL_ERROR_NONE) {
    muse_imagelist_delete(images);
    return NULL;
  }
  muse_basicproc_quadrant_statistics(images, aProcessing); /* non-fatal */
  if (muse_basicproc_overscans_compute_stats(images, aBPars) /* incl. CR rejection */
      != CPL_ERROR_NONE) {
    muse_imagelist_delete(images); /* this failure should be fatal */
    return NULL;
  }
  if (muse_basicproc_correct_overscans_vpoly(images, aBPars)
      != CPL_ERROR_NONE) {
    muse_imagelist_delete(images); /* this failures should be fatal! */
    return NULL;
  }
  if (muse_basicproc_trim_images(images) != CPL_ERROR_NONE) {
    muse_imagelist_delete(images);
    return NULL;
  }
  if (muse_basicproc_correct_overscans_offset(images, aProcessing, aBPars)
      != CPL_ERROR_NONE) { /* for BIAS files only */
    cpl_msg_warning(__func__, "Bias-level correction using overscans failed in "
                    "IFU %hhu: %s", aIFU, cpl_error_get_message());
  }
  if (muse_basicproc_check_overscans(images, aProcessing, aBPars)
      != CPL_ERROR_NONE) { /* only for BIAS files and when overscan==none */
    muse_imagelist_delete(images);
    return NULL;
  }
  muse_basicproc_gain_override(images, aProcessing, aIFU);
  muse_basicproc_check_gain(images, aProcessing, aIFU); /* for flats */
  if (muse_basicproc_apply_bias(images, aProcessing, aIFU, aBPars)
      != CPL_ERROR_NONE) {
    muse_imagelist_delete(images);
    return NULL;
  }
  if (muse_basicproc_corr_nonlinearity(images, aProcessing, aIFU) != CPL_ERROR_NONE) {
    muse_imagelist_delete(images);
    return NULL;
  }
  if (muse_basicproc_adu_to_count(images, aProcessing) != CPL_ERROR_NONE) {
    muse_imagelist_delete(images);
    return NULL;
  }
  if (muse_basicproc_apply_dark(images, aProcessing, aIFU) != CPL_ERROR_NONE) {
    muse_imagelist_delete(images);
    return NULL;
  }
  if (muse_basicproc_apply_cr(images, aBPars) != CPL_ERROR_NONE) {
    muse_imagelist_delete(images);
    return NULL;
  }
  if (muse_basicproc_apply_flat(images, aProcessing, aBPars, aIFU) != CPL_ERROR_NONE) {
    muse_imagelist_delete(images);
    return NULL;
  }
  return images;
} /* muse_basicproc_load() */

/*---------------------------------------------------------------------------*/
/**
  @brief  Load reduced input files from disk.
  @param  aProcessing   the processing structure
  @param  aIFU          the IFU/channel number
  @return the muse_imagelist * or NULL on error

  Read the reduced input files from disk and store them as images in the images
  structure.

  @error{propagate error from muse_processing_check_input()\, return NULL,
         input data is invalid somehow (NULL aProcessing pointer etc.)}
 */
/*---------------------------------------------------------------------------*/
muse_imagelist *
muse_basicproc_load_reduced(muse_processing *aProcessing, unsigned char aIFU)
{
  muse_imagelist *images = muse_imagelist_new();
  cpl_frameset *redframes = muse_frameset_find_tags(aProcessing->inframes,
                                                    aProcessing->intags,
                                                    aIFU, CPL_FALSE);
  cpl_size iframe, nframes = cpl_frameset_get_size(redframes);
  for (iframe = 0; iframe < nframes; iframe++) {
    cpl_frame *frame = cpl_frameset_get_position(redframes, iframe);
    cpl_errorstate prestate = cpl_errorstate_get();
    const char *imagename = cpl_frame_get_filename(frame);
    muse_image *image = muse_image_load(imagename);
    if (!image) {
      cpl_errorstate_set(prestate);
      image = muse_image_load_from_extensions(imagename, aIFU);
    }
    muse_imagelist_set(images, image, iframe);
    muse_processing_append_used(aProcessing, frame, CPL_FRAME_GROUP_RAW, 1);
  } /* for iframe (all used frames) */
  cpl_frameset_delete(redframes);
  return images;
} /* muse_basicproc_load_reduced() */

/*---------------------------------------------------------------------------*/
/**
  @private
  @brief  Prepare an illum/attached flat-field to usable form from its pixel
          table.
  @param  aPT   pixel table of an attached flat-field exposure
  @return the prepared CPL table or NULL on failure

  This function converts the pixel table of the attached flat-field into the
  form expected by muse_basicproc_apply_attached(), with columns "slice" and
  "fflat".

  @error{set CPL_ERROR_NULL_INPUT\, return NULL,
         aPT or one of it's components are NULL}
 */
/*---------------------------------------------------------------------------*/
static cpl_table *
muse_basicproc_prepare_illum(muse_pixtable *aPT)
{
  cpl_ensure(aPT && aPT->header && aPT->table,
             CPL_ERROR_NULL_INPUT, NULL);

  /* crop pixel table to small wavelength range *
   * and separate into per-slice tables         */
  muse_pixtable_restrict_wavelength(aPT, 6500., 7500.);
  muse_pixtable **pts = muse_pixtable_extracted_get_slices(aPT);
  int ipt, npt = muse_pixtable_extracted_get_size(pts);

  unsigned char ifu = muse_utils_get_ifu(aPT->header);
  cpl_msg_info(__func__, "Preparing %s flat: %d slices in the data of IFU %hhu "
               "found.", MUSE_TAG_ILLUM, npt, ifu);
  cpl_table *tattached = cpl_table_new(npt);
  cpl_table_new_column(tattached, "slice", CPL_TYPE_INT);
  cpl_table_new_column(tattached, "fflat", CPL_TYPE_DOUBLE);
  for (ipt = 0; ipt < npt; ipt++) {
    uint32_t origin = cpl_table_get_int(pts[ipt]->table, MUSE_PIXTABLE_ORIGIN, 0, NULL);
    unsigned short slice = muse_pixtable_origin_get_slice(origin);
    double median = cpl_table_get_column_median(pts[ipt]->table,
                                                MUSE_PIXTABLE_DATA);
    cpl_msg_debug(__func__, "Found median of %f in slice %d of IFU %hhu "
                  "of illum flat.", median, slice, ifu);
    cpl_table_set_int(tattached, "slice", ipt, slice);
    cpl_table_set_double(tattached, "fflat", ipt, 1. / median);
  } /* for ipt */
  muse_pixtable_extracted_delete(pts);
  /* normalize the flat-field scales */
  double mean = cpl_table_get_column_mean(tattached, "fflat");
  cpl_msg_debug(__func__, "Normalizing whole illum-flat table if IFU %hhu to "
                "%e.", ifu, mean);
  cpl_table_multiply_scalar(tattached, "fflat", 1. / mean);
  cpl_table_set_column_format(tattached, "fflat", "%.6f");
  return tattached;
} /* muse_basicproc_prepare_illum() */

/*---------------------------------------------------------------------------*/
/**
  @brief  Get an illum/attached flat-field from an imagelist and prepare it for
          use.
  @param  aImages   image list to search for
  @param  aTrace    trace table
  @param  aWave     wavelength calibration table
  @param  aGeo      geometry table
  @return the prepared CPL table or NULL on failure

  This function finds and erases all attached flat-fields or illumination
  flat-fields from the input image lists. The first ILLUM found is then
  converted to a table that can be used by muse_basicproc_apply_illum() to
  correct a pixel table of a science exposure by the actual illumination.

  @error{set CPL_ERROR_NULL_INPUT\, return NULL,
         one of the input pointers is NULL}
 */
/*---------------------------------------------------------------------------*/
cpl_table *
muse_basicproc_get_illum(muse_imagelist *aImages, cpl_table *aTrace,
                         cpl_table *aWave, cpl_table *aGeo)
{
  cpl_ensure(aImages && aTrace && aWave && aGeo, CPL_ERROR_NULL_INPUT, NULL);

  cpl_table *tattached = NULL;
  unsigned int k, nimages = muse_imagelist_get_size(aImages);
  /* also create list of illum exposures to be able to erase them at the end */
  cpl_boolean *isillum = cpl_calloc(nimages, sizeof(cpl_boolean));
  for (k = 0; k < nimages; k++) {
    isillum[k] = CPL_FALSE;
    muse_image *image = muse_imagelist_get(aImages, k);
    const char *tag = cpl_propertylist_get_string(image->header,
                                                  MUSE_HDR_TMP_INTAG);
    if (tag && !strncmp(tag, MUSE_TAG_ILLUM, strlen(MUSE_TAG_ILLUM) + 1)) {
      isillum[k] = CPL_TRUE; /* is an attached FLAT,ILLUM */
      /* also verify the template ID */
      if (cpl_propertylist_has(image->header, "ESO TPL ID")) {
        const char *tplid = cpl_propertylist_get_string(image->header,
                                                        "ESO TPL ID"),
                   *fn = cpl_propertylist_get_string(image->header,
                                                     MUSE_HDR_TMP_FN),
                   *tplatt = "MUSE_wfm_cal_specflatatt",
                   *tplill = "MUSE_wfm_cal_illum",
                   *tpliln = "MUSE_nfm_cal_illum";
        if (strncmp(tplid, tplatt, strlen(tplatt) + 1) &&
            strncmp(tplid, tplill, strlen(tplill) + 1) &&
            strncmp(tplid, tpliln, strlen(tpliln) + 1)) {
          cpl_msg_warning(__func__, "%s input (\"%s\") was taken with neither"
                          " %s nor %s template, but %s!", MUSE_TAG_ILLUM, fn,
                          tplatt, tplill, tplid);
        } else {
          cpl_msg_debug(__func__, "%s input (\"%s\") was taken with template "
                        "%s", MUSE_TAG_ILLUM, fn, tplid);
        } /* else */
      } /* if TPL.ID present */
    } /* if */
    unsigned char ifu = muse_utils_get_ifu(image->header);
    if (isillum[k]) {
      if (tattached) {
        cpl_msg_warning(__func__, "Image %u of %u of IFU %hhu is illum flat, "
                        "but not the first; not using it!", k + 1, nimages, ifu);
        continue;
      }
      cpl_msg_debug(__func__, "Image %u of %u of IFU %hhu is illum flat.",
                    k + 1, nimages, ifu);
      muse_pixtable *pt = muse_pixtable_create(image, aTrace, aWave, aGeo);
      tattached = muse_basicproc_prepare_illum(pt);
      muse_pixtable_delete(pt);
    } else {
      cpl_msg_debug(__func__, "Image %u of %u of IFU %hhu is not an illum flat.",
                    k + 1, nimages, ifu);
    }
  } /* for k */

  /* remove the ILLUM image(s) from the image list */
  unsigned int k2;
  for (k = 0, k2 = 0; k < nimages; k++) {
    if (isillum[k]) {
      muse_image *image = muse_imagelist_unset(aImages, k2);
      muse_image_delete(image);
    } else {
      k2++; /* only step, if the image was not removed */
    }
  } /* for k, k2 */
  cpl_free(isillum);

  return tattached;
} /* muse_basicproc_get_illum() */

/*---------------------------------------------------------------------------*/
/**
  @brief  Apply an illum/attached flat-field to a pixel table.
  @param  aPT         the pixel table
  @param  aAttached   the table computed from the attached flat
  @return CPL_ERROR_NONE on success or a CPL error code on failure

  This function expects the attached flat-field to be in form of a table, in
  the format created by muse_basicproc_prepare_attached().

  @error{set and return CPL_ERROR_NULL_INPUT,
         aPT\, one of it's components\, or aAttached are NULL}
 */
/*---------------------------------------------------------------------------*/
cpl_error_code
muse_basicproc_apply_illum(muse_pixtable *aPT, cpl_table *aAttached)
{
  cpl_ensure_code(aPT && aPT->header && aPT->table && aAttached,
                  CPL_ERROR_NULL_INPUT);

  unsigned char ifu = muse_utils_get_ifu(aPT->header);

  muse_pixtable **pts = muse_pixtable_extracted_get_slices(aPT);
  int ipt, npt = muse_pixtable_extracted_get_size(pts);
  cpl_msg_info(__func__, "Applying %s flat-field in IFU %hhu (%d slices)",
               MUSE_TAG_ILLUM, ifu, npt);
  cpl_array *afactors = cpl_array_new(npt, CPL_TYPE_DOUBLE);
  for (ipt = 0; ipt < npt; ipt++) {
    uint32_t origin = cpl_table_get_int(pts[ipt]->table, MUSE_PIXTABLE_ORIGIN, 0, NULL);
    unsigned short slice = muse_pixtable_origin_get_slice(origin),
                   fslice = cpl_table_get_int(aAttached, "slice", ipt, NULL);
    int err;
    double fflat = cpl_table_get_double(aAttached, "fflat", ipt, &err);
    if (!err && slice == fslice) {
      cpl_table_multiply_scalar(pts[ipt]->table, MUSE_PIXTABLE_DATA, fflat);
      cpl_table_multiply_scalar(pts[ipt]->table, MUSE_PIXTABLE_STAT, fflat*fflat);
      cpl_array_set_double(afactors, ipt, fflat);
      char *kw = cpl_sprintf(MUSE_HDR_PT_ILLUMi, slice);
      cpl_propertylist_update_double(aPT->header, kw, fflat);
      cpl_free(kw);
    } else {
      cpl_msg_warning(__func__, "some error (%d) occurred when applying illum-"
                      "flat correction to slice %hu/%hu of IFU %hhu: %s", err,
                      slice, fslice, ifu, cpl_error_get_message());
    } /* else */
  } /* for ipt */
  muse_pixtable_extracted_delete(pts);
  cpl_propertylist_update_double(aPT->header, MUSE_HDR_PT_ILLUM_MEAN,
                                 cpl_array_get_mean(afactors));
  cpl_propertylist_update_double(aPT->header, MUSE_HDR_PT_ILLUM_STDEV,
                                 cpl_array_get_stdev(afactors));
  cpl_array_delete(afactors);
  return CPL_ERROR_NONE;
} /* muse_basicproc_apply_illum() */

/*---------------------------------------------------------------------------*/
/**
  @brief  Apply an attached flat-field to a pixel table.
  @param  aPT         the pixel table
  @param  aTwilight   the cube of the twilight skyflat exposure
  @return CPL_ERROR_NONE on success or a CPL error code on failure

  This function applies the twilight-sky based illumination correction using a
  cube of the skyflat prepared by the muse_twilight recipe.

  Spatially, the geometry-table based spatial pixel coordinates are used to find
  the nearest neighbor. In dispersion direction, the (up to two pixels) in the
  nearest wavelength plane(s) are interpolated linearly.

  @error{set and return CPL_ERROR_NULL_INPUT,
         aPT\, one of it's components\, or aAttached are NULL}
 */
/*---------------------------------------------------------------------------*/
cpl_error_code
muse_basicproc_apply_twilight(muse_pixtable *aPT, muse_datacube *aTwilight)
{
  cpl_ensure_code(aPT && aPT->header && aPT->table && aTwilight,
                  CPL_ERROR_NULL_INPUT);

  // XXX should check, if the twilight cube is of the same INS.MODE as the data

  unsigned char ifu = muse_utils_get_ifu(aPT->header);

  /* transfer the sky flat flux value */
  char *kw = cpl_sprintf(MUSE_HDR_FLAT_FLUX_SKY"%hhu", ifu);
  double flux = cpl_propertylist_get_double(aTwilight->header, kw);
  cpl_free(kw);
  cpl_propertylist_update_double(aPT->header, MUSE_HDR_FLAT_FLUX_SKY, flux);

  /* get the WCS from the twilight cube */
  int nx = cpl_image_get_size_x(cpl_imagelist_get(aTwilight->data, 0)),
      ny = cpl_image_get_size_y(cpl_imagelist_get(aTwilight->data, 0)),
      nz = cpl_imagelist_get_size(aTwilight->data);
  cpl_msg_debug(__func__, "Working with %d planes in twilight cube", nz);
  double cd12 = muse_pfits_get_cd(aTwilight->header, 1, 2),
         cd21 = muse_pfits_get_cd(aTwilight->header, 2, 1);
  if (cd12 > DBL_EPSILON || cd21 > DBL_EPSILON) {
    cpl_msg_warning(__func__, "Twilight cube contains WCS cross-terms (CD1_2"
                    "=%e, CD2_1=%e), results will be inaccurate!", cd12, cd21);
  }
  // XXX guard against more rubbish (CTYPEi, CUNITi, cross-terms with axis3
  double crval1 = muse_pfits_get_crval(aTwilight->header, 1),
         crpix1 = muse_pfits_get_crpix(aTwilight->header, 1),
         cd11 = muse_pfits_get_cd(aTwilight->header, 1, 1),
         crval2 = muse_pfits_get_crval(aTwilight->header, 2),
         crpix2 = muse_pfits_get_crpix(aTwilight->header, 2),
         cd22 = muse_pfits_get_cd(aTwilight->header, 2, 2),
         crval3 = muse_pfits_get_crval(aTwilight->header, 3),
         crpix3 = muse_pfits_get_crpix(aTwilight->header, 3),
         cd33 = muse_pfits_get_cd(aTwilight->header, 3, 3);

  /* loop through the pixel table and apply the correction */
  float *data = cpl_table_get_data_float(aPT->table, MUSE_PIXTABLE_DATA),
        *stat = cpl_table_get_data_float(aPT->table, MUSE_PIXTABLE_STAT),
        *xpos = cpl_table_get_data_float(aPT->table, MUSE_PIXTABLE_XPOS),
        *ypos = cpl_table_get_data_float(aPT->table, MUSE_PIXTABLE_YPOS),
        *lbda = cpl_table_get_data_float(aPT->table, MUSE_PIXTABLE_LAMBDA);
  cpl_size irow, nrow = muse_pixtable_get_nrow(aPT),
           nfailed = 0;
  for (irow = 0; irow < nrow; irow++) {
    /* find closest spatial pixel in twilight cube *
     * using the coordinates from the pixel table      */
    int x = lround((xpos[irow] - crval1) / cd11 + crpix1), /* nearest neighbor */
        y = lround((ypos[irow] - crval2) / cd22 + crpix2); /* nearest neighbor */
    /* boundary conditions */
    if (x < 1) {
      x = 1;
    }
    if (x > nx) {
      x = nx;
    }
    if (y < 1) {
      y = 1;
    }
    if (y > ny) {
      y = ny;
    }
    double z = (lbda[irow] - crval3) / cd33 + crpix3; /* plane */
#if 0
    cpl_msg_debug(__func__, "%"CPL_SIZE_FORMAT": %.3f %.3f %.3f => %d %d %.3f",
                  irow, xpos[irow], ypos[irow], lbda[irow], x, y, z);
#endif
    /* get indices of the two closest image planes for linear interpolation */
    int zp1 = floor(z) - 1,
        zp2 = ceil(z) - 1;
    /* boundary conditions */
    if (zp1 < 1) {
      zp1 = 0;
    }
    if (zp1 >= nz) {
      zp1 = nz - 1;
    }
    if (zp2 < 1) {
      zp2 = 0;
    }
    if (zp2 >= nz) {
      zp2 = nz - 1;
    }
    int err1, err2;
    double v1 = cpl_image_get(cpl_imagelist_get(aTwilight->data, zp1),
                              x, y, &err1),
           v2 = cpl_image_get(cpl_imagelist_get(aTwilight->data, zp2),
                              x, y, &err2);
    double villum = 1.;
    double f = 1.;
    /* linearly interpolate the given twilight factors from both planes */
    if (err1 && err2) {
      nfailed++;
      continue;
    }
    if (zp1 == zp2) {
      villum = v1; /* same planes, just take the first value */
    } else if (err1) {
      villum = v2;
    } else if (err2) {
      villum = v1;
    } else { /* real interpolation */
      f = fabs(z - 1 - zp1);
      villum = v1 * (1. - f) + v2 * f;
    }
    double fillum = 1. / villum; /* the inverse as correction factor */
#if 0
    cpl_msg_debug(__func__, "%d/%d, %f/%f -> %f -> %f => %f", zp1, zp2, v1, v2,
                  f, villum, fillum);
#endif

    /* multiply the data value by the inverse twilight factor */
    data[irow] *= fillum;
    /* multiply the stat value by the squared inverse twilight factor */
    stat[irow] *= fillum*fillum;
  } /* for irow (all pixel table rows) */
  if (nfailed) {
    cpl_msg_warning(__func__, "Failed to correct twilight in %"CPL_SIZE_FORMAT
                    " of %"CPL_SIZE_FORMAT", pixels in IFU %hhu!", nfailed,
                    nrow, ifu);
  } else {
    cpl_msg_debug(__func__, "No failures during twilight correction of %"
                  CPL_SIZE_FORMAT" pixels in IFU %hhu", nrow, ifu);
  }

  return CPL_ERROR_NONE;
} /* muse_basicproc_apply_twilight() */

/*----------------------------------------------------------------------------*/
/**
  @brief  Mask the range of the NaD notch filter in the given pixel table.
  @param  aPT   the pixel table to mask the range in
  @param  aIFU   the IFU number (for output)
  @return CPL_ERROR_NONE on success, and another CPL error code on failure.

  This function selects the range of wavelengths inside the notch filter (given
  by kMuseNa2LambdaMin and kMuseNa2LambdaMax or related constants, depending on
  the exact instrument mode) and sets the DQ entries to EURO3D_MISSDATA.

  The function returns without clearing the table row selection.

  @error{set and return CPL_ERROR_NULL_INPUT,
         aPT or one of its components is NULL}
  @error{set and return CPL_ERROR_ILLEGAL_INPUT,
         instrument mode without notch filter found}
 */
/*----------------------------------------------------------------------------*/
cpl_error_code
muse_basicproc_mask_notch_filter(muse_pixtable *aPT, unsigned char aIFU)
{
  cpl_ensure_code(aPT && aPT->header && aPT->table, CPL_ERROR_NULL_INPUT);

  /* for the moment set values for [WFM-AO-N] (Na2 filter) upfront */
  float lmin = kMuseNa2LambdaMin,
        lmax = kMuseNa2LambdaMax;
  muse_ins_mode insmode = muse_pfits_get_mode(aPT->header);
  const char *mode = muse_pfits_get_insmode(aPT->header);
  if (insmode == MUSE_MODE_WFM_AO_N) {
    /* nothing to be done */
  } else if (insmode == MUSE_MODE_WFM_AO_E) {
    lmin = kMuseNaLambdaMin;
    lmax = kMuseNaLambdaMax;
  } else if (insmode == MUSE_MODE_NFM_AO_N) {
    lmin = kMuseNaGLambdaMin;
    lmax = kMuseNaGLambdaMax;
  } else {
    cpl_msg_warning(__func__, "No notch filter for mode %s!", mode);
    return CPL_ERROR_ILLEGAL_INPUT;
  }
  cpl_msg_info(__func__, "%s mode: marking NaD region (%.1f..%.1f Angstrom) of "
               "IFU %d as 0x%08lx", mode, lmin, lmax, aIFU, EURO3D_NOTCH_NAD);
  cpl_table_unselect_all(aPT->table);
  cpl_table_or_selected_float(aPT->table, MUSE_PIXTABLE_LAMBDA,
                              CPL_GREATER_THAN, lmin);
  cpl_table_and_selected_float(aPT->table, MUSE_PIXTABLE_LAMBDA,
                               CPL_LESS_THAN, lmax);
  cpl_array *asel = cpl_table_where_selected(aPT->table);
  cpl_size isel, nsel = cpl_array_get_size(asel);
  const cpl_size *sel = cpl_array_get_data_cplsize_const(asel);
  int *dq = cpl_table_get_data_int(aPT->table, MUSE_PIXTABLE_DQ);
  for (isel = 0; isel < nsel; isel++) {
    dq[sel[isel]] = EURO3D_NOTCH_NAD;
  } /* for isel (all selected table rows) */
  cpl_array_delete(asel);

  return CPL_ERROR_NONE;
} /* muse_basicproc_mask_notch_filter() */

/*----------------------------------------------------------------------------*/
/**
  @private
  @brief  Check if two frames are equal in terms of lamps being switched on.
  @param  aFrame1   the first frame to compare
  @param  aFrame2   the second frame to compare
  @return 1 if the frames are identical, 0 if they are different, and -1 on
          error.

  This function is to be used as the comparison function in conjunction with
  cpl_frameset_labelise(). It checks the headers of two frames to see if the
  same lamps are switched on and the same lamp shutters are open.

  In principle, one needs to check that the shutters have the same names as the
  lamps with the same number, but due to the software setup this is guaranteed
  to be true. So this function does not make any effort to parse the keyword
  strings to verify this.

  @error{set CPL_ERROR_NULL_INPUT\, return -1,
         one of the input arguments is NULL}
  @error{propagate error code\, return -1, loading one of the headers fails}
  @error{set CPL_ERROR_INCOMPATIBLE_INPUT\, return -1,
         same lamp number has different name in the two exposures}
  @error{set CPL_ERROR_INCOMPATIBLE_INPUT\, return -1,
         same shutter number has different name in the two exposures}
 */
/*----------------------------------------------------------------------------*/
static int
muse_basicproc_combine_compare_lamp(const cpl_frame *aFrame1, const cpl_frame *aFrame2)
{
  cpl_ensure(aFrame1 && aFrame2, CPL_ERROR_NULL_INPUT, -1);
  const char *fn1 = cpl_frame_get_filename(aFrame1),
             *fn2 = cpl_frame_get_filename(aFrame2);
  cpl_propertylist *head1 = cpl_propertylist_load(fn1, 0),
                   *head2 = cpl_propertylist_load(fn2, 0);
  if (!head1 || !head2) {
    cpl_propertylist_delete(head1); /* in case one was loaded... */
    cpl_propertylist_delete(head2);
    return -1;
  }

  /* Loop through all lamps in the header and find their status. The first     *
   * missing shutter entry will cause the FITS query to throw an error, that's *
   * when we stop. Otherwise, we are done when the lamp status is not equal.   */
  int nlamp = 1, status1, status2;
  cpl_errorstate prestate = cpl_errorstate_get();
  do {
    /* ensure that we are dealing with the same lamps! */
    const char *name1 = muse_pfits_get_lamp_name(head1, nlamp),
               *name2 = muse_pfits_get_lamp_name(head2, nlamp);
    cpl_errorstate_set(prestate); /* lamps may be missing */
    if (name1 && name2 && strncmp(name1, name2, strlen(name1) + 1)) {
      cpl_error_set_message(__func__, CPL_ERROR_INCOMPATIBLE_INPUT,
                            "Files \"%s\" and \"%s\" have incompatible lamp "
                            "setups", fn1, fn2);
      cpl_propertylist_delete(head1);
      cpl_propertylist_delete(head2);
      return -1;
    }
    name1 = muse_pfits_get_shut_name(head1, nlamp);
    name2 = muse_pfits_get_shut_name(head2, nlamp);
    if (name1 && name2 && strncmp(name1, name2, strlen(name1) + 1)) {
      cpl_error_set_message(__func__, CPL_ERROR_INCOMPATIBLE_INPUT,
                            "Files \"%s\" and \"%s\" have incompatible shutter "
                            "setups", fn1, fn2);
      cpl_propertylist_delete(head1);
      cpl_propertylist_delete(head2);
      return -1;
    }

    status1 = muse_pfits_get_lamp_status(head1, nlamp);
    status2 = muse_pfits_get_lamp_status(head2, nlamp);
    cpl_errorstate_set(prestate); /* lamps may be missing */
    if (status1 != status2) {
      break;
    }
    status1 = muse_pfits_get_shut_status(head1, nlamp);
    status2 = muse_pfits_get_shut_status(head2, nlamp);
    if (status1 != status2) {
      break;
    }
    nlamp++;
  } while (cpl_errorstate_is_equal(prestate));
  cpl_errorstate_set(prestate);

  cpl_propertylist_delete(head1);
  cpl_propertylist_delete(head2);
  return status1 == status2;
} /* muse_basicproc_combine_compare_lamp() */

/*---------------------------------------------------------------------------*/
/**
  @brief  Combine several images into a lampwise image list.
  @param  aProcessing      the processing structure
  @param  aIFU             the IFU/channel number
  @param  aBPars           basic processing parameters
  @param  aLabeledFrames   optional output array of framesets containing used
                           frames for each lamp
  @return the list as a muse_imagelist * or NULL on error

  Combine several images into one, using method and parameters specified by
  aProcessing->parameters. Note that this function changes the input images if
  the "scale" option is used.  If all images are found to be of the same lamp,
  they are all combined into a single output image, which is the only entry in
  the returned image list.

  If aLabeledFrames is not NULL, it will contain a list of framesets or the
  same size as the returned imagelist. Each frameset contains the frames that
  were used to preprocess the image at the same index in the image list, and
  has to be deallocated using cpl_frameset_delete(). The returned pointer has
  to be deallocated with cpl_free(). The returned pointer will be NULL on error.

  @error{set CPL_ERROR_NULL_INPUT\, return NULL, aProcessing is NULL}
  @error{call muse_basicproc_load() with NULL parameter, aBPars is NULL}
  @error{propagate error code\, return NULL, image combination failed}
 */
/*---------------------------------------------------------------------------*/
muse_imagelist *
muse_basicproc_combine_images_lampwise(muse_processing *aProcessing,
                                       unsigned char aIFU,
                                       muse_basicproc_params *aBPars,
                                       cpl_frameset ***aLabeledFrames)
{
  if (aLabeledFrames) { /* NULL out return pointer, in case it's given... */
    *aLabeledFrames = NULL; /* ... so it's always NULL in case of problems */
  }
  cpl_ensure(aProcessing, CPL_ERROR_NULL_INPUT, NULL);

  /* find out different lamps in input */
  cpl_frameset *rawframes = muse_frameset_find_tags(aProcessing->inframes,
                                                    aProcessing->intags, aIFU,
                                                    CPL_FALSE);
  char *prefix = cpl_sprintf("muse.%s", aProcessing->name);
  muse_combinepar *cpars = muse_combinepar_new(aProcessing->parameters,
                                               prefix);
  cpl_free(prefix);
#if 0
  printf("rawframes\n");
  cpl_frameset_dump(rawframes, stdout);
  fflush(stdout);
#endif
  cpl_size nlabels,
           *labels = cpl_frameset_labelise(rawframes,
                                           muse_basicproc_combine_compare_lamp,
                                           &nlabels);
  if (!labels || nlabels <= 1) {
    /* if labeling didn't work, process the list of one lamp and return */
    cpl_free(labels);
    cpl_frameset_delete(rawframes);
    muse_imagelist *list = muse_basicproc_load(aProcessing, aIFU, aBPars),
                   *images = NULL;
    if (nlabels == 1) {
      muse_image *image = muse_combine_images(cpars, list);
      images = muse_imagelist_new();
      muse_imagelist_set(images, image, 0);
      if (aLabeledFrames) {
        *aLabeledFrames = cpl_calloc(1, sizeof(cpl_frameset *));
        (*aLabeledFrames)[0] = cpl_frameset_duplicate(aProcessing->usedframes);
      } /* if */
    } /* if only one label */
    muse_imagelist_delete(list);
    muse_combinepar_delete(cpars);
    return images;
  } /* if one of no labels */

#if 0
  cpl_array *alabels = cpl_array_wrap_int(labels, cpl_frameset_get_size(rawframes));
  cpl_array_dump(alabels, 0, 1000, stdout);
  fflush(stdout);
  cpl_array_unwrap(alabels);
#endif
  /* output list of lampwise combined images */
  muse_imagelist *images = muse_imagelist_new();
  if (aLabeledFrames) {
    *aLabeledFrames = cpl_calloc(nlabels, sizeof(cpl_frameset *));
  }

  /* duplicate aProcessing into a local structure the contents of which  *
   * we can manipulate here (don't change it directly for threadsafety!) */
  muse_processing *proc = cpl_malloc(sizeof(muse_processing));
  memcpy(proc, aProcessing, sizeof(muse_processing));
  cpl_frameset *inframes = proc->inframes;
  /* copy frames with all extra frames somewhere else */
  cpl_frameset *auxframes = muse_frameset_find_tags(inframes, aProcessing->intags,
                                                    aIFU, CPL_TRUE);
  /* loop through labels for all lamps */
  int ilabel, ipos = 0;
  for (ilabel = 0; ilabel < nlabels; ilabel++) {
    /* create new sub-frameset for this lamp */
    cpl_frameset *frames = cpl_frameset_extract(rawframes, labels, ilabel);
    /* append other files for the initial processing */
    cpl_frameset_join(frames, auxframes);
    /* substitute aProcessing->inframes */
    proc->inframes = frames;
    /* load and combine frames for this sub-frameset */
    muse_imagelist *list = muse_basicproc_load(proc, aIFU, aBPars);
    /* reinstate original aProcessing->inframes */
    proc->inframes = inframes;
    if (!list) { /* break this loop to fail the function below */
      muse_imagelist_delete(images);
      cpl_frameset_delete(frames);
      images = NULL;
      /* if muse_basicproc_load() fails then because of some missing *
       * calibration, and then it will already fail for the first    *
       * ilabel; it is therefore enough to free the pointer          */
      if (aLabeledFrames) {
        cpl_free(*aLabeledFrames);
        *aLabeledFrames = NULL;
      }
      break;
    }

    muse_image *lampimage = muse_combine_images(cpars, list);
    if (!lampimage) {
      cpl_msg_error(__func__, "Image combination failed for IFU %hhu for lamp "
                    "with label %d of %"CPL_SIZE_FORMAT, aIFU, ilabel + 1, nlabels);
      muse_imagelist_delete(list);
      cpl_frameset_delete(frames);
      continue;
    }

    if (aLabeledFrames) {
      /* if a given frame was used now, copy its group */
      cpl_size iframe, nframes = cpl_frameset_get_size(frames);
      for (iframe = 0; iframe < nframes; iframe++) {
        cpl_frame *frame = cpl_frameset_get_position(frames, iframe);
        const char *fn = cpl_frame_get_filename(frame),
                   *tag = cpl_frame_get_tag(frame);
        cpl_size iuframe, nuframes = cpl_frameset_get_size(aProcessing->usedframes);
        for (iuframe = 0;
             (iuframe < nuframes) && fn && tag; /* only check with valid info */
             iuframe++) {
          cpl_frame *uframe = cpl_frameset_get_position(aProcessing->usedframes,
                                                        iuframe);
          const char *fnu = cpl_frame_get_filename(uframe),
                     *tagu = cpl_frame_get_tag(uframe);
          if (fnu && !strncmp(fn, fnu, strlen(fn) + 1) &&
              tagu && !strncmp(tag, tagu, strlen(tag) + 1)) {
            cpl_frame_set_group(frame, cpl_frame_get_group(uframe));
            break;
          }
        } /* for uframe (all usedframes) */
      } /* for frame */
      (*aLabeledFrames)[ipos] = frames;
    } else {
      cpl_frameset_delete(frames);
    }

    /* transfer NSATURATION headers from invidual images to lamp-combined image */
    unsigned int k;
    for (k = 0; k < muse_imagelist_get_size(list); k++) {
      char *keyword = cpl_sprintf(QC_WAVECAL_PREFIXi" "QC_BASIC_NSATURATED, k+1);
      int nsaturated = cpl_propertylist_get_int(muse_imagelist_get(list, k)->header,
                                                MUSE_HDR_TMP_NSAT);
      cpl_propertylist_update_int(lampimage->header, keyword, nsaturated);
      cpl_free(keyword);
    }
    muse_imagelist_delete(list);
    /* append to imagelist */
    muse_imagelist_set(images, lampimage, ipos++);
  } /* for ilabel (labels) */
  cpl_free(labels);
  cpl_free(proc);
  muse_combinepar_delete(cpars);
  cpl_frameset_delete(rawframes);
  cpl_frameset_delete(auxframes);

  if (images && muse_imagelist_get_size(images) == 0) {
    muse_imagelist_delete(images);
    images = NULL;
    if (aLabeledFrames) {
      cpl_free(*aLabeledFrames);
      *aLabeledFrames = NULL;
    }
  }

  return images;
} /* muse_basicproc_combine_images_lampwise() */

/*---------------------------------------------------------------------------*/
/**
  @brief Compute wavelength corrections for science data based on reference sky
         lines.
  @param aPt          the pixel table with the science data
  @param aLines       array with the reference wavelengths of sky emission lines
  @param aHalfWidth   half-width of the wavelength region around each sky line
  @param aBinWidth    pixel size in Angstrom of the intermediate spectrum
  @param aLo          low sigma-clipping limit for the intermediate spectrum
  @param aHi          high sigma-clipping limit for the intermediate spectrum
  @param aIter        number of iterations for sigma-clipping the spectrum
  @return CPL_ERROR_NONE on success or another CPL error code on failure

  This function subsequently calls @ref muse_utils_pixtable_fit_line_gaussian()
  for all sky lines given in aLines. It computes a weighted mean offset and
  applies this shift in wavelength to the data in aPt.

  The input array aLines needs to be of simple floating-point type.

  @error{return CPL_ERROR_NULL_INPUT, aPt and/or aLines are NULL}
  @error{return CPL_ERROR_ILLEGAL_INPUT, aLines is not of floating-point type}
 */
/*---------------------------------------------------------------------------*/
cpl_error_code
muse_basicproc_shift_pixtable(muse_pixtable *aPt, cpl_array *aLines,
                              double aHalfWidth, double aBinWidth,
                              float aLo, float aHi, unsigned char aIter)
{
  cpl_ensure_code(aPt && aLines, CPL_ERROR_NULL_INPUT);
  cpl_ensure_code(cpl_array_get_type(aLines) == CPL_TYPE_DOUBLE ||
                  cpl_array_get_type(aLines) == CPL_TYPE_FLOAT,
                  CPL_ERROR_ILLEGAL_INPUT);

  double lmin = cpl_propertylist_get_float(aPt->header, MUSE_HDR_PT_LLO),
         lmax = cpl_propertylist_get_float(aPt->header, MUSE_HDR_PT_LHI);
  double shift = 0., wshift = 0.;
  cpl_array *errors = cpl_array_new(4, CPL_TYPE_DOUBLE);
  int i, n = cpl_array_get_size(aLines), nvalid = 0;
  for (i = 0; i < n; i++) {
    int err;
    double lambdasign = cpl_array_get(aLines, i, &err),
           lambda = fabs(lambdasign);
    if (err || lambda >= lmax || lambda <= lmin) {
      cpl_msg_debug(__func__, "Invalid wavelength at position %d of %d in "
                    "skylines", i + 1, n);
      continue;
    }
    nvalid++;
    double center = muse_utils_pixtable_fit_line_gaussian(aPt, lambdasign, aHalfWidth,
                                                          aBinWidth, aLo, aHi, aIter,
                                                          NULL, errors),
           cerr = cpl_array_get_double(errors, 0, NULL);
    shift += (lambda - center) / cerr;
    wshift += 1. / cerr;
    cpl_msg_debug(__func__, "dlambda = %.4f +/- %.4f (for skyline at %.4f "
                  "Angstrom)", lambda - center, cerr, lambda);
  } /* for i (all entries in aLines) */
  cpl_array_delete(errors);
  shift /= wshift;
  if (nvalid > 0 && isfinite(shift)) {
    cpl_msg_info(__func__, "Skyline correction (%d lines): shifting data of IFU "
                 "%hhu by %.4f Angstrom", nvalid, muse_utils_get_ifu(aPt->header),
                 shift);
    cpl_table_add_scalar(aPt->table, MUSE_PIXTABLE_LAMBDA, shift);
    cpl_propertylist_update_float(aPt->header, QC_SCIBASIC_SHIFT, shift);
  } else {
    cpl_propertylist_update_float(aPt->header, QC_SCIBASIC_SHIFT, 0.);
  }
  return CPL_ERROR_NONE;
} /* muse_basicproc_shift_pixtable() */

/*---------------------------------------------------------------------------*/
/**
  @brief Compute image statistics of an image and add them to a header
  @param aImage    the image to derive statistics for and add them to
  @param aHeader   the header to add them to
  @param aPrefix   the prefix of the output FITS keywords
  @param aStats    the statistics properties to compute
  @return a CPL error code on failure or CPL_ERROR_NONE on success

  This function only supports a subset of the properties supported by
  cpl_stats_new_from_image(), see muse_basicproc_stats_append_header_window().
 */
/*---------------------------------------------------------------------------*/
cpl_error_code
muse_basicproc_stats_append_header(cpl_image *aImage, cpl_propertylist *aHeader,
                                   const char *aPrefix, unsigned aStats)
{
  cpl_ensure_code(aImage, CPL_ERROR_NULL_INPUT);

  int nx = cpl_image_get_size_x(aImage),
      ny = cpl_image_get_size_y(aImage);
  return muse_basicproc_stats_append_header_window(aImage, aHeader, aPrefix,
                                                   aStats, 1, 1, nx, ny);
} /* muse_basicproc_stats_append_header() */

/*---------------------------------------------------------------------------*/
/**
  @brief Compute image statistics of an image window and add them to a header
  @param aImage    the image to derive statistics for and add them to
  @param aHeader   the header to add them to
  @param aPrefix   the prefix of the output FITS keywords
  @param aStats    the statistics properties to compute
  @param aX1       the lower left x-coordinate of the window
  @param aY1       the lower left y-coordinate of the window
  @param aX2       the upper right x-coordinate of the window
  @param aY2       the upper right y-coordinate of the window
  @return a CPL error code on failure or CPL_ERROR_NONE on success

  This function only supports a subset of the properties supported by
  cpl_stats_new_from_image_window(), namely median, mean, stdev, min, max,
  and flux,
 */
/*---------------------------------------------------------------------------*/
cpl_error_code
muse_basicproc_stats_append_header_window(cpl_image *aImage,
                                          cpl_propertylist *aHeader,
                                          const char *aPrefix, unsigned aStats,
                                          int aX1, int aY1, int aX2, int aY2)
{
  cpl_ensure_code(aImage && aHeader, CPL_ERROR_NULL_INPUT);
  cpl_ensure_code(aPrefix, CPL_ERROR_ILLEGAL_INPUT);
  cpl_ensure_code(aPrefix, CPL_ERROR_ILLEGAL_INPUT);

  cpl_stats *stats = cpl_stats_new_from_image_window(aImage, aStats,
                                                     aX1, aY1, aX2, aY2);
  if (!stats) {
    return cpl_error_get_code();
  }

  char keyword[KEYWORD_LENGTH];
  if (aStats & CPL_STATS_MEDIAN) {
    snprintf(keyword, KEYWORD_LENGTH, "%s MEDIAN", aPrefix);
    cpl_propertylist_append_float(aHeader, keyword,
                                  cpl_stats_get_median(stats));
  }
  if (aStats & CPL_STATS_MEAN) {
    snprintf(keyword, KEYWORD_LENGTH, "%s MEAN", aPrefix);
    cpl_propertylist_append_float(aHeader, keyword, cpl_stats_get_mean(stats));
  }
  if (aStats & CPL_STATS_STDEV) {
    snprintf(keyword, KEYWORD_LENGTH, "%s STDEV", aPrefix);
    cpl_propertylist_append_float(aHeader, keyword, cpl_stats_get_stdev(stats));
  }
  if (aStats & CPL_STATS_MIN) {
    snprintf(keyword, KEYWORD_LENGTH, "%s MIN", aPrefix);
    cpl_propertylist_append_float(aHeader, keyword, cpl_stats_get_min(stats));
  }
  if (aStats & CPL_STATS_MAX) {
    snprintf(keyword, KEYWORD_LENGTH, "%s MAX", aPrefix);
    cpl_propertylist_append_float(aHeader, keyword, cpl_stats_get_max(stats));
  }
  if (aStats & CPL_STATS_FLUX) {
    snprintf(keyword, KEYWORD_LENGTH, "%s INTFLUX", aPrefix);
    cpl_propertylist_append_float(aHeader, keyword, cpl_stats_get_flux(stats));
  }

  cpl_stats_delete(stats);

  return CPL_ERROR_NONE;
} /* muse_basicproc_stats_append_header_window() */

/*---------------------------------------------------------------------------*/
/**
  @brief  Add QC parameter about saturated pixels to a muse_image.
  @param  aImage    the image to count on and modify
  @param  aPrefix   prefix of the QC header
  @return CPL_ERROR_NONE on success, another CPL error code on failure.

  Count EURO3D_SATURATED pixels in the dq component of aImage and add a
  aPrefix+QC_BASIC_NSATURATED (if aPrefix does not contain a trailing space,
  one is added in between) keyword to its header component.

  @error{return CPL_ERROR_NULL_INPUT,
         aImage or its dq or header components\, or aPrefix are NULL}
  @error{propagate error code, keyword updating failed}
 */
/*---------------------------------------------------------------------------*/
cpl_error_code
muse_basicproc_qc_saturated(muse_image *aImage, const char *aPrefix)
{
  cpl_ensure_code(aImage && aImage->dq && aImage->header && aPrefix,
                  CPL_ERROR_NULL_INPUT);

  cpl_mask *mask = cpl_mask_threshold_image_create(aImage->dq,
                                                   EURO3D_SATURATED - 0.1,
                                                   EURO3D_SATURATED + 0.1);
  int nsaturated = cpl_mask_count(mask);
  cpl_mask_delete(mask);
  /* check if the prefix has a trailing space, add it if not */
  char *keyword = NULL;
  if (aPrefix[strlen(aPrefix)-1] == ' ') {
    keyword = cpl_sprintf("%s%s", aPrefix, QC_BASIC_NSATURATED);
  } else {
    keyword = cpl_sprintf("%s %s", aPrefix, QC_BASIC_NSATURATED);
  }
  cpl_error_code rc = cpl_propertylist_update_int(aImage->header, keyword,
                                                  nsaturated);
  cpl_free(keyword);
  return rc;
} /* muse_basicproc_qc_saturated() */

/**@}*/