/*                                                                              *
 *   This file is part of the ESO UVES Pipeline                                 *
 *   Copyright (C) 2004,2005 European Southern Observatory                      *
 *                                                                              *
 *   This library 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, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA       *
 *                                                                              */


/*
 * $Author: amodigli $
 * $Date: 2013-04-16 15:46:57 $
 * $Revision: 1.65 $
 * $Name: not supported by cvs2svn $
 * $Log: not supported by cvs2svn $
 * Revision 1.64  2013/04/05 08:02:31  amodigli
 * changed WAVELENGTH to AWAV to be FITS compliant
 *
 * Revision 1.63  2012/05/02 06:08:32  amodigli
 * replace Ang by Angstrom
 *
 * Revision 1.62  2012/03/02 16:53:31  amodigli
 * fixed warning related to upgrade to CPL6
 *
 * Revision 1.61  2011/12/08 14:03:09  amodigli
 * Fix warnings with CPL6
 *
 * Revision 1.60  2010/12/16 16:57:40  amodigli
 * fixed compiler warnings
 *
 * Revision 1.59  2010/12/08 11:07:59  amodigli
 * added chip parameter to uves_merge_orders() to have proper filenames. Fixed content err data for noappend case
 *
 * Revision 1.58  2010/09/27 06:32:22  amodigli
 * fixed mem leaks
 *
 * Revision 1.57  2010/09/24 09:32:04  amodigli
 * put back QFITS dependency to fix problem spot by NRI on FIBER mode (with MIDAS calibs) data
 *
 * Revision 1.55  2010/06/07 09:48:00  amodigli
 * changed units: A-->Ang, FLUX-->ADU
 *
 * Revision 1.54  2010/06/01 16:06:55  amodigli
 * added unit [A] to reduced product
 *
 * Revision 1.53  2008/09/29 06:57:41  amodigli
 * add #include <string.h>
 *
 * Revision 1.52  2008/09/27 16:05:12  amodigli
 * fixed bug using delta parmeters
 *
 * Revision 1.51  2008/09/23 11:32:39  amodigli
 * added check on array upper bound access to fix DFS05803
 *
 * Revision 1.50  2008/08/29 09:52:39  amodigli
 * fixed compiler warning
 *
 * Revision 1.49  2008/06/26 08:30:38  amodigli
 * fixed bug in setting delta
 *
 * Revision 1.48  2008/06/11 14:42:50  amodigli
 * fixed seg fault
 *
 * Revision 1.45  2008/03/04 07:35:31  amodigli
 * generate spectra of each order only if NOAPPEND
 *
 * Revision 1.44  2008/02/21 07:50:38  amodigli
 * added method NOAPPEND
 *
 * Revision 1.43  2007/08/21 13:08:26  jmlarsen
 * Removed irplib_access module, largely deprecated by CPL-4
 *
 * Revision 1.42  2007/06/21 11:28:57  jmlarsen
 * Added support for type float (for FLAMES)
 *
 * Revision 1.41  2007/06/06 08:17:33  amodigli
 * replace tab with 4 spaces
 *
 * Revision 1.40  2007/05/03 15:21:13  jmlarsen
 * Decreased output message verbosity
 *
 * Revision 1.39  2007/04/24 12:50:29  jmlarsen
 * Replaced cpl_propertylist -> uves_propertylist which is much faster
 *
 * Revision 1.38  2007/02/09 08:57:54  jmlarsen
 * Include <float.h>
 *
 * Revision 1.37  2007/02/09 08:14:16  jmlarsen
 * Do not use CPL_PIXEL_MAXVAL which works only for integer images
 *
 * Revision 1.36  2007/02/08 07:34:28  jmlarsen
 * Minor doc change
 *
 * Revision 1.35  2006/11/15 15:02:14  jmlarsen
 * Implemented const safe workarounds for CPL functions
 *
 * Revision 1.33  2006/11/15 14:04:08  jmlarsen
 * Removed non-const version of parameterlist_get_first/last/next which is already
 * in CPL, added const-safe wrapper, unwrapper and deallocator functions
 *
 * Revision 1.32  2006/11/13 14:23:55  jmlarsen
 * Removed workarounds for CPL const bugs
 *
 * Revision 1.31  2006/11/06 15:19:41  jmlarsen
 * Removed unused include directives
 *
 * Revision 1.30  2006/09/11 14:00:11  jmlarsen
 * Minor documentation change
 *
 * Revision 1.29  2006/08/17 13:56:53  jmlarsen
 * Reduced max line length
 *
 * Revision 1.28  2006/08/17 09:16:47  jmlarsen
 * Removed CPL2 code
 *
 * Revision 1.27  2006/08/11 14:37:59  jmlarsen
 * Added input validation
 *
 * Revision 1.26  2006/08/10 10:50:12  jmlarsen
 * Removed workaround for cpl_image_get_bpm
 *
 * Revision 1.25  2006/07/03 13:16:42  jmlarsen
 * Reduced number of significant digits in message
 *
 * Revision 1.24  2006/04/06 08:38:56  jmlarsen
 * Changed char* -> const char* for static string
 *
 * Revision 1.23  2006/03/03 13:54:11  jmlarsen
 * Changed syntax of check macro
 *
 * Revision 1.22  2006/02/03 07:46:30  jmlarsen
 * Moved recipe implementations to ./uves directory
 *
 * Revision 1.21  2006/01/31 08:24:29  jmlarsen
 * Wrapper for cpl_image_get_bpm
 *
 * Revision 1.20  2006/01/25 16:13:20  jmlarsen
 * Changed interface of gauss.fitting routine
 *
 * Revision 1.19  2005/12/19 16:17:56  jmlarsen
 * Replaced bool -> int
 *
 * Revision 1.18  2005/12/16 14:22:23  jmlarsen
 * Removed midas test data; Added sof files
 *
 * Revision 1.17  2005/11/24 15:09:06  jmlarsen
 * Implemented 2d extraction/rebinning/merging
 *
 * Revision 1.16  2005/11/24 11:54:46  jmlarsen
 * Added support for CPL 3 interface
 *
 * Revision 1.15  2005/11/18 10:52:06  jmlarsen
 * Split into optimal/sum merge methods
 *
 * Revision 1.14  2005/11/11 13:18:54  jmlarsen
 * Reorganized code, renamed source files
 *
 */

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

/*----------------------------------------------------------------------------*/
/**
 * @defgroup uves_merge   Order merging
 *
 * Merge echelle orders to obtain a 1D spectrum.
 */
/*----------------------------------------------------------------------------*/
/**@{*/


/*-----------------------------------------------------------------------------
                                Includes
 -----------------------------------------------------------------------------*/

#include <uves_merge.h>

#include <uves_pfits.h>
#include <uves_utils.h>
#include <uves_utils_wrappers.h>
#include <uves_dump.h>
#include <uves_dfs.h>
#include <uves_error.h>

#include <cpl.h>
#include <float.h>
#include <string.h>
/*-----------------------------------------------------------------------------
                            Functions prototypes
 ----------------------------------------------------------------------------*/

/*-----------------------------------------------------------------------------
                            Implementation
 ----------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------*/
/**
  @brief    Merge orders
  @param    spectrum        The 2D spectrum in (wavelength, order)-space
  @param    spectrum_noise  The spectrum noise
  @param    spectrum_header Header describing the wavelength offset for each row
                            in the spectrum image.
  @param    m_method        The merge method
  @param    n_traces        Number of traces per order (equal to 1, or more if 2d extraction)
  @param    merged_header  (out) Newly allocated header of merged image
  @param    merged_noise   (out) Newly allocated merged 1D noise spectrum
  @return   The merged spectrum as an image of height @em n_traces.
  
  If @em m_method is MERGE_OPTIMAL, the resulting spectrum is at each
  wavelength optimally merged from
  the (one, two or more) single order spectra:

  flux      = sum [1/sigma_i^2  flux_i] / sum [1/sigma_i^2]
  sigma^2   = 1 / sum [1/sigma_i^2]

  If @em m_method is MERGE_SUM, the sum (not average as above) of input
  fluxes is computed:

  flux    = sum [flux_i]
  sigma^2 = sum [sigma_i^2]

 */
/*---------------------------------------------------------------------------*/

cpl_image *
uves_merge_orders(const cpl_image *spectrum, 
                  const cpl_image *spectrum_noise,
                  const uves_propertylist *spectrum_header,
                  merge_method m_method,
                  int n_traces,
                  uves_propertylist **merged_header,
                  const double delt1, 
                  const double delt2,
                  enum uves_chip chip,
                  cpl_image **merged_noise)
{
    cpl_image *merged = NULL;    /* Result */

    const double *spectrum_data_double = NULL;
    const float  *spectrum_data_float  = NULL;
    const cpl_mask *spectrum_badmap    = NULL;
    const cpl_binary *spectrum_bad     = NULL;

    const double *noise_data_double = NULL;
    const float *noise_data_float   = NULL;
    const cpl_mask *noise_badmap    = NULL;
    const cpl_binary *noise_bad     = NULL;
    cpl_type type;               /* input/output images type */

    int nbins, ny, norders;          /* Input image size. ny = norders*n_traces */
    double wavestep;
    int bin_min = 0, bin_max = 0;    /* wavelength of min/max bin in units of 'wavestep' */
    int total_bins;
    int order, trace;
    int spectrum_sx=0;
    int spectrum_sy=0;
    double delt1_bin=0;
    double delt2_bin=0;
    char* filename=NULL;

    cpl_vector* image_1d=NULL;
    uves_propertylist* hext=NULL;
    //double med_noise=0;
    passure( spectrum != NULL, " ");
    passure( spectrum_noise != NULL, " ");
    passure( spectrum_header != NULL, " ");
    passure( merged_header != NULL, " ");
    passure( merged_noise != NULL, " ");

    assure( m_method == MERGE_OPTIMAL || 
            m_method == MERGE_SUM || 
            m_method == MERGE_NOAPPEND, 
            CPL_ERROR_ILLEGAL_INPUT,
            "Unknown merge method: %d", m_method);
    
    assure( cpl_image_get_type(spectrum) == CPL_TYPE_DOUBLE ||
            cpl_image_get_type(spectrum) == CPL_TYPE_FLOAT,
            CPL_ERROR_TYPE_MISMATCH,
            "Spectrum must have type double or float. It is '%s'",
            uves_tostring_cpl_type(cpl_image_get_type(spectrum)));
    
    assure( cpl_image_get_type(spectrum_noise) == CPL_TYPE_DOUBLE ||
            cpl_image_get_type(spectrum_noise) == CPL_TYPE_FLOAT,
            CPL_ERROR_TYPE_MISMATCH,
            "Spectrum noise must have type double. It is '%s'",
            uves_tostring_cpl_type(cpl_image_get_type(spectrum_noise)));

    assure( cpl_image_get_type(spectrum) ==
            cpl_image_get_type(spectrum_noise),
            CPL_ERROR_TYPE_MISMATCH,
            "Spectrum and spectrum noise must have same type. They are "
            "%s and %s, respectively",
            uves_tostring_cpl_type(cpl_image_get_type(spectrum)),
            uves_tostring_cpl_type(cpl_image_get_type(spectrum_noise)) );

    type = cpl_image_get_type(spectrum);

    /* Read input spectrum geometry */
    nbins           = cpl_image_get_size_x(spectrum);
    ny              = cpl_image_get_size_y(spectrum);

    assure( cpl_image_get_size_x(spectrum_noise) == nbins &&
            cpl_image_get_size_y(spectrum_noise) == ny,
            CPL_ERROR_INCOMPATIBLE_INPUT,
            "Incompatible spectrum/noise image sizes: %dx%d vs. %" CPL_SIZE_FORMAT "x%" CPL_SIZE_FORMAT "",
            nbins, ny,
            cpl_image_get_size_x(spectrum_noise),
            cpl_image_get_size_y(spectrum_noise) );
    
    assure( ny % n_traces == 0, CPL_ERROR_INCOMPATIBLE_INPUT,
            "Spectrum image height (%d) is not a multiple of "
            "the number of traces (%d). Confused, bailing out",
            ny, n_traces);

    norders         = ny / n_traces;
    
    check( wavestep = uves_pfits_get_cdelt1(spectrum_header),
       "Error reading bin width");

    /* Get data pointers (for efficiency) */
   
    if (type == CPL_TYPE_DOUBLE) {
        spectrum_data_double   = cpl_image_get_data_double_const(spectrum);
    }
    else {
        spectrum_data_float   = cpl_image_get_data_float_const(spectrum);
    }

    spectrum_sx=cpl_image_get_size_x(spectrum);
    spectrum_sy=cpl_image_get_size_y(spectrum);

    spectrum_badmap = cpl_image_get_bpm_const(spectrum);
    spectrum_bad    = cpl_mask_get_data_const(spectrum_badmap);
    
    if (type == CPL_TYPE_DOUBLE) {
        noise_data_double = cpl_image_get_data_double_const(spectrum_noise);        
    }
    else {
        noise_data_float  = cpl_image_get_data_float_const(spectrum_noise);
    }
    noise_badmap = cpl_image_get_bpm_const(spectrum_noise);
    noise_bad    = cpl_mask_get_data_const(noise_badmap);
    uves_msg("delt1=%f delt2=%f",delt1,delt2);
    /* Read max/min lambda */
    for (order = 1; order <= norders; order++)
    {
        double wstart, wend;
        check( wstart = uves_pfits_get_wstart(spectrum_header, order),
           "Error reading start wavelength for order #%d", order);

        check( wend = uves_pfits_get_wend(spectrum_header, order),
           "Error reading end wavelength for order #%d", order);
        
        uves_msg_debug("Order #%d: wstart - wend = %f - %f wlu", order, wstart, wend);

	//wstart+=delt1;
	//wend-=delt2;

        if (order == 1)
        {
            bin_min = uves_round_double(wstart/wavestep);
            bin_max = uves_round_double(wend  /wavestep);
        }
        
        bin_min = uves_min_int(bin_min, uves_round_double(wstart/wavestep));
        bin_max = uves_max_int(bin_max, uves_round_double(wend  /wavestep));
    }
    total_bins = (bin_max - bin_min) + 1;

    uves_msg_debug("Merging orders into %d bins covering wavelengths %.3f - %.3f wlu", 
            total_bins, bin_min * wavestep, bin_max * wavestep);

    /* Initialize spectrum to zero and noise to negative */
    if(m_method == MERGE_NOAPPEND) {
      merged        = cpl_image_new(total_bins, n_traces*norders, type);
      *merged_noise = cpl_image_new(total_bins, n_traces*norders, type);
    } else {
      merged        = cpl_image_new(total_bins, n_traces, type);
      *merged_noise = cpl_image_new(total_bins, n_traces, type);
    }
    cpl_image_add_scalar(*merged_noise, -1.0);
    //cpl_image_power(spectrum_noise,0.5);
    //cpl_image_multiply_scalar(spectrum_noise,0.25);
    /* Distribute input in output bins */
    for (order = 1; order <= norders; order++)
    {
        double wstart, wend;
        int wstart_bin, wend_bin;      /* In 1d space */



        check( wstart = uves_pfits_get_wstart(spectrum_header, order),
           "Error reading start wavelength for order #%d", order);
        
        check( wend = uves_pfits_get_wend(spectrum_header, order),
           "Error reading end wavelength for order #%d", order);

	//wstart+=delt1;
	//wend-=delt2;
        
        wstart_bin = uves_round_double(wstart/wavestep);
        wend_bin   = uves_round_double(wend/wavestep);
        delt1_bin  = uves_round_double(delt1/wavestep);
        delt2_bin  = uves_round_double(delt2/wavestep);

        
	int bin_min_ord = uves_round_double(wstart/wavestep);
	int bin_max_ord = uves_round_double(wend  /wavestep);
        int nbins_ord = (bin_max_ord - bin_min_ord) + 1;
	  cpl_image* merged_ord=NULL;
	  cpl_image * noise_ord=NULL;

	if(m_method == MERGE_NOAPPEND) {
	  merged_ord = cpl_image_new(nbins_ord, n_traces, type);
	  noise_ord = cpl_image_new(nbins_ord, n_traces, type);
	}

        //if(order>1 && order<norders) {
	//  med_noise=cpl_image_get_median_window(spectrum_noise,
	//					1,order,ny,order);
	//} else {
	//  med_noise=cpl_image_get_median(spectrum_noise);
	//}
        /* Loop over spatial traces (only 1 trace, unless extraction was 2d) */
        for (trace = 1; trace <= n_traces; trace++)
        {
	  int merged_row = 0;
	  int spectrum_row = (order - 1)*n_traces + trace;
	  if(m_method == MERGE_NOAPPEND) {
	    merged_row = (order - 1)*n_traces + trace;
	  } else {
	    merged_row = trace;

	  }
            int rel_bin;                   /* Counting columns in input spectrum */
            
            for (rel_bin = 1+delt1_bin; 
                 (rel_bin <= wend_bin - wstart_bin + 1-delt2_bin) &&
                 (rel_bin <(spectrum_sx*spectrum_sy+1-(spectrum_row-1)*nbins));
                 rel_bin++)
            {
                double flux, noise;
                double current_flux, new_flux;
                double current_noise, new_noise;
                double weight;
                int pis_rejected, noise_rejected;
                
                /* merged_bin = (offset of order)  +  (offset inside order) */
                int merged_bin = (wstart_bin - bin_min) + rel_bin;
                int merged_bin_ord = rel_bin;
                
                passure(1 <= merged_bin && merged_bin <= total_bins,
                    "%d %d %d", rel_bin, merged_bin, total_bins);
                
                /* This is slow:
                   check( flux  = cpl_image_get(spectrum      , 
                          rel_bin, spectrum_row, &pis_rejected);
                   noise = cpl_image_get(spectrum_noise, rel_bin, 
                    spectrum_row, &noise_rejected),
                   "Error reading input spectrum");
                */
                
                if (type == CPL_TYPE_DOUBLE) {
                    flux  = spectrum_data_double[(rel_bin-1) + (spectrum_row-1) * nbins];
                    noise = noise_data_double   [(rel_bin-1) + (spectrum_row-1) * nbins];

                }
                else {
                    flux  = spectrum_data_float[(rel_bin-1) + (spectrum_row-1) * nbins];
                    noise = noise_data_float   [(rel_bin-1) + (spectrum_row-1) * nbins];
                }

                pis_rejected   = spectrum_bad[(rel_bin-1) + (spectrum_row-1) * nbins];
                noise_rejected = noise_bad   [(rel_bin-1) + (spectrum_row-1) * nbins];
                
                if (!pis_rejected && !noise_rejected)
                {

		  if(m_method == MERGE_NOAPPEND) {
                    check(( current_flux  = cpl_image_get(merged, 
                                                          merged_bin, 
                                                          merged_row, 
							  &pis_rejected),
			    current_noise = cpl_image_get(*merged_noise, 
							  merged_bin, 
                                                          merged_row, 
							  &pis_rejected)),
			  "Error reading merged spetrum");
		  } else {
                    check(( current_flux  = cpl_image_get(
                        merged      , merged_bin, trace, &pis_rejected),
                        current_noise = cpl_image_get(
                        *merged_noise, merged_bin, trace, &pis_rejected)),
                       "Error reading merged spetrum");
		  }
		  weight = 1/(noise*noise);
                    
                    /*
                     * Optimal formulas for Variance and Flux are
                     *
                     *   Vn = ( 1/sigma1^2 + ... +  1/sigmaN^2)^-1
                     *   Fn = (f1/sigma1^2 + ... + fN/sigmaN^2) * Vn
                     *
                     * Update by using these recurrence relations
                     *
                     *   Fn+1 = (Fn/Vn + fn+1/(sigma_{n+1})^2) * Vn+1  for n > 1
                     *   Vn+1 = (Vn^-1 + 1/(sigma_{n+1})^2)^-1  for n > 1
                     *
                     *
                     *  In the case of method = sum,
                     *
                     *  Vn = sigma1^2 + ... + sigmaN^2
                     *  Fn = f1 + ... + fN
                     *
                     */
                    
                    if (current_noise > 0)
                    {
                        if (m_method == MERGE_OPTIMAL)
                        {
                            new_noise  = 1/(current_noise*current_noise);
                            new_noise += weight;
                            new_noise  = 1/sqrt(new_noise);
                        }
                        else if (m_method == MERGE_SUM)
                        {
                            new_noise = sqrt(current_noise*current_noise
                                     + noise*noise);
                        }
                        else if (m_method == MERGE_NOAPPEND)
                        {
                            new_noise = current_noise;
                        }
                        else
                        {
                            /* Impossible */
                            passure( false, "%d", m_method);
                        }
                    }
                    else
                    {
                        /* First time in this bin */
                        new_noise = noise;
                    }
                    
                    if (current_noise > 0)
                    {
                        if (m_method == MERGE_OPTIMAL)
                        {
                            new_flux = (current_flux / 
                                (current_noise*current_noise)
                                + flux * weight) * 
                            (new_noise*new_noise);
                        }
                        else if (m_method == MERGE_SUM)
                        {
                            new_flux = current_flux + flux;
                        }
                        else if (m_method == MERGE_NOAPPEND)
			  {
                            new_flux = flux;
			  }
                        else
                        {
                            /* Impossible */
                            passure( false, "%d", m_method);
                        }
                    }
                    else
                    {
                        new_flux = flux;
                    }
                    
		  if(m_method == MERGE_NOAPPEND) {

		    /*
		      uves_msg_warning("flux[%d,%d]=%g noise=%g",
				       merged_bin,merged_row,
				       new_flux,current_noise);
		    */
                    check( cpl_image_set(
					 merged, 
                                         merged_bin, 
                                         merged_row, 
                                         new_flux),
                       "Error updating merged spectrum");

                    check( cpl_image_set(
					 *merged_noise, 
                                          merged_bin, 
					 merged_row, 
                                         new_noise),
                       "Error updating weights");


                    check( cpl_image_set(
					 merged_ord, 
                                         merged_bin_ord, 
                                         trace, 
                                         new_flux),
                       "Error updating merged spectrum");

                    check( cpl_image_set(
					 noise_ord, 
                                         merged_bin_ord, 
                                         trace, 
                                         new_noise),
                       "Error updating merged spectrum");


		  } else {
                    check( cpl_image_set(
                           merged      , merged_bin, trace, new_flux),
                       "Error updating merged spectrum");
                    check( cpl_image_set(
                           *merged_noise, merged_bin, trace, new_noise),
                       "Error updating weights");
		  }
/*            uves_msg("Input flux = %e +- %e ;  
            Binned flux changed from %e +- %e to %e +- %e",
            flux, noise,
            current_flux, current_noise, new_flux, new_noise);*/
                    
                } /* If pixel is good ... */

            } /* For each input bin */

        }/* For trace ... */


	if (merged_header == NULL)
        {
           uves_free_propertylist(merged_header);
           check( *merged_header = uves_initialize_image_header(
                     "AWAV", (n_traces > 1) ? "PIXEL" : " ",
                     "Angstrom", (n_traces > 1) ? "PIXEL" : NULL, 
                     "ADU",0,
                     bin_min_ord * wavestep, 1.0,
                     1.0, 1.0,
                     wavestep, 1.0),
                  "Error initializing merged spectrum header");
        }

        if(m_method == MERGE_NOAPPEND) {
           filename=uves_sprintf("merged_data_noappend_%s.fits",
                           uves_chip_tostring_upper(chip));
           image_1d = cpl_vector_wrap(cpl_image_get_size_x(merged_ord),
                                      cpl_image_get_data_double(merged_ord));
           uves_free_propertylist(&hext);
           hext=uves_propertylist_new();
           uves_propertylist_append_double(hext,"CRVAL1",wstart-delt1);
           uves_propertylist_append_double(hext,"CDELT1",wavestep);
           uves_propertylist_append_double(hext,"CRPIX1",0);

           if(order==1) {
            uves_vector_save(image_1d,filename,CPL_BPP_IEEE_FLOAT,hext, 
                                    CPL_IO_DEFAULT);
           } else {
            uves_vector_save(image_1d,filename,CPL_BPP_IEEE_FLOAT,hext, 
                                    CPL_IO_EXTEND);
           }
           cpl_vector_unwrap(image_1d);

           image_1d = cpl_vector_wrap(cpl_image_get_size_x(noise_ord),
                                      cpl_image_get_data_double(noise_ord));

           uves_free(filename);
      
           filename=uves_sprintf("merged_sigma_noappend_%s.fits",
                                 uves_chip_tostring_upper(chip));
           if(order==1) {
            uves_vector_save(image_1d,filename,CPL_BPP_IEEE_FLOAT,hext, 
                                    CPL_IO_DEFAULT);
           } else {
            uves_vector_save(image_1d,filename,CPL_BPP_IEEE_FLOAT,hext, 
                                    CPL_IO_EXTEND);
           }
           cpl_vector_unwrap(image_1d);

           uves_free(filename);
           uves_free_image(&merged_ord);
           uves_free_image(&noise_ord);
           uves_free_propertylist(&hext);
           uves_free_propertylist(merged_header);
        }

    }/* For order ... */
    
    /* Undefined bins have  (flux, noise) = (0, -1)  (the initial values).
     * Set to (flux, noise) = (0, 1)
     */

/* Here commented out piece of code to be used eventually to have
   noappend files files as products
    if(m_method == MERGE_NOAPPEND) {
       cpl_frame* mdata_noappend=NULL;
       cpl_frame* msigma_noappend=NULL;
       mdata_noappend=cpl_frame_new();
       msigma_noappend=cpl_frame_new();
       cpl_frame_set_group(mdata_noappend,CPL_FRAME_GROUP_PRODUCT);
       cpl_frame_set_type(mdata_noappend,CPL_FRAME_TYPE_IMAGE);
       cpl_frame_set_level(mdata_noappend,CPL_FRAME_LEVEL_FINAL);

       tag=uves_sprintf("MERGED_SCI_NOAPPEND_%s",
                             uves_chip_tostring_upper(chip));
       filename=uves_sprintf("%s_%s.fits",tag,uves_chip_tostring_upper(chip));
       cpl_frame_set_tag(mdata_noappend,tag);
       cpl_frame_set_filename(mdata_noappend,filename);
       uves_free(filename);
       uves_free(tag);

       cpl_frame_set_group(msigma_noappend,CPL_FRAME_GROUP_PRODUCT);
       cpl_frame_set_type(msigma_noappend,CPL_FRAME_TYPE_IMAGE);
       cpl_frame_set_level(msigma_noappend,CPL_FRAME_LEVEL_FINAL);

       tag=uves_sprintf("MERGED_ERR_NOAPPEND_%s",
                             uves_chip_tostring_upper(chip));
       filename=uves_sprintf("%s_%s.fits",tag,uves_chip_tostring_upper(chip));
       cpl_frame_set_tag(msigma_noappend,tag);
       cpl_frame_set_filename(msigma_noappend,filename);

       uves_free_frame(&mdata_noappend);
       uves_free_frame(&msigma_noappend);
       uves_free(filename);
       uves_free(tag);

    }
*/
    check( cpl_image_threshold(*merged_noise, 
                   0, DBL_MAX,        /* Outside this interval */
                   1, 1),             /* Set to these values   */
       "Error setting undefined noise");
 
    
    if (merged_header != NULL)
    {
        check( *merged_header = uves_initialize_image_header(
               "AWAV", (n_traces > 1) ? "PIXEL" : " ", 
               "Angstrom", (n_traces > 1) ? "PIXEL" : NULL, 
               "ADU",0,
               bin_min * wavestep, 1.0,
               1.0, 1.0,
               wavestep, 1.0),
           "Error initializing merged spectrum header");
    }
 
    
  cleanup:
    return merged;
}

/*---------------------------------------------------------------------------*/
/**
  @brief    Read merging method from parameter list
  @param    parameters             The parameter list
  @param    context                Context of parameter (or NULL)
  @param    subcontext             Subcontext of parameter
                                           
  @return   The merging method as read from the parameter context.subcontext.merge
  
**/
/*---------------------------------------------------------------------------*/
merge_method
uves_get_merge_method(const cpl_parameterlist *parameters, 
                      const char *context, 
                      const char *subcontext)
{
    const char *mm = "";
    merge_method result = 0;
    
    check( uves_get_parameter(parameters, context, subcontext, "merge", CPL_TYPE_STRING, &mm),
       "Could not read parameter");
    
    if      (strcmp(mm, "optimal") == 0) result = MERGE_OPTIMAL;
    else if (strcmp(mm, "sum"    ) == 0) result = MERGE_SUM;
    else if (strcmp(mm, "noappend") == 0) result = MERGE_NOAPPEND;
    else
    {
        assure(false, CPL_ERROR_ILLEGAL_INPUT, "No such merging method: '%s'", mm);
    }
    
  cleanup:
    return result;
}
    

/**@}*/