/*                                                                              *
 *   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: 2012-05-02 06:11:40 $
 * $Revision: 1.38 $
 * $Name: not supported by cvs2svn $
 * $Log: not supported by cvs2svn $
 * Revision 1.37  2012/03/02 16:40:40  amodigli
 * fixed warning related to upgrade to CPL6
 *
 * Revision 1.36  2011/12/08 14:00:02  amodigli
 * Fox warnings with CPL6
 *
 * Revision 1.35  2011/04/14 11:25:40  amodigli
 * fixed typo QC key  in comments
 *
 * Revision 1.34  2011/04/11 09:07:41  amodigli
 * implemented QC comments corrections from DFO
 *
 * Revision 1.33  2011/04/11 07:53:12  amodigli
 * uniformed QC param key name
 *
 * Revision 1.32  2011/03/23 12:27:31  amodigli
 * changed QC key as user likes
 *
 * Revision 1.31  2011/03/23 10:08:47  amodigli
 * added QC to better characterize wave accuracy
 *
 * Revision 1.30  2010/09/24 09:32:09  amodigli
 * put back QFITS dependency to fix problem spot by NRI on FIBER mode (with MIDAS calibs) data
 *
 * Revision 1.28  2007/07/23 14:57:30  jmlarsen
 * Make workaround work
 *
 * Revision 1.27  2007/07/23 12:40:37  jmlarsen
 * Update to CPL4
 *
 * Revision 1.26  2007/06/06 08:17:33  amodigli
 * replace tab with 4 spaces
 *
 * Revision 1.25  2007/05/22 11:46:15  jmlarsen
 * Removed 1d wavecal mode which was not supported
 *
 * Revision 1.24  2007/05/16 16:33:42  amodigli
 * fixed leak
 *
 * Revision 1.23  2007/05/10 08:32:48  jmlarsen
 * Minor output message change
 *
 * Revision 1.22  2007/05/07 14:26:44  jmlarsen
 * Added QC.NLINSOL parameter
 *
 * Revision 1.21  2007/05/07 07:13:59  jmlarsen
 * Made resolution computation robust against negative dl/dx
 *
 * Revision 1.20  2007/04/27 07:22:57  jmlarsen
 * Implemented possibility to use automatic polynomial degree
 *
 * Revision 1.19  2007/04/13 07:34:54  jmlarsen
 * Removed dead code
 *
 * Revision 1.18  2007/04/10 07:12:09  jmlarsen
 * Changed interface of polynomial_regression_2d()
 *
 * Revision 1.17  2007/03/15 12:36:44  jmlarsen
 * Added experimental ppm code
 *
 * Revision 1.16  2007/03/05 10:24:14  jmlarsen
 * Do kappa-sigma rejection only in second loop
 *
 * Revision 1.15  2007/02/22 15:37:35  jmlarsen
 * Use kappa-sigma clipping when fitting dispersion
 *
 * Revision 1.14  2007/01/15 08:58:51  jmlarsen
 * Added text output
 *
 * Revision 1.13  2006/11/06 15:19:42  jmlarsen
 * Removed unused include directives
 *
 * Revision 1.12  2006/10/12 11:36:48  jmlarsen
 * Reduced max line length
 *
 * Revision 1.11  2006/10/10 11:20:11  jmlarsen
 * Renamed line table columns to match MIDAS
 *
 * Revision 1.10  2006/08/17 14:11:25  jmlarsen
 * Use assure_mem macro to check for memory allocation failure
 *
 * Revision 1.9  2006/08/17 13:56:53  jmlarsen
 * Reduced max line length
 *
 * Revision 1.8  2006/08/11 14:36:37  jmlarsen
 * Added profiling info
 *
 * Revision 1.7  2006/08/07 11:35:08  jmlarsen
 * Removed hardcoded constant
 *
 * Revision 1.6  2006/07/14 12:52:57  jmlarsen
 * Exported/renamed function find_nearest
 *
 * Revision 1.5  2006/07/14 12:44:26  jmlarsen
 * Use less significant digits
 *
 * Revision 1.4  2006/04/24 09:33:48  jmlarsen
 * Shortened max line length
 *
 * Revision 1.3  2006/03/03 13:54:11  jmlarsen
 * Changed syntax of check macro
 *
 * Revision 1.2  2006/02/15 13:19:15  jmlarsen
 * Reduced source code max. line length
 *
 * Revision 1.1  2006/02/03 07:46:30  jmlarsen
 * Moved recipe implementations to ./uves directory
 *
 * Revision 1.31  2005/12/20 08:11:44  jmlarsen
 * Added CVS  entry
 *
 */

/*----------------------------------------------------------------------------*/
/**
 * @addtogroup uves_wavecal
 */
/*----------------------------------------------------------------------------*/
/**@{*/

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

#include <uves_wavecal_identify.h>

#include <uves_wavecal_utils.h>
#include <uves_utils.h>
#include <uves_utils_wrappers.h>
#include <uves_error.h>
#include <uves_msg.h>
#include <cpl_ppm.h> /* missing from cpl.h */
#include <uves_qclog.h>
#include <cpl.h>

#include <math.h>
#include <float.h>

#define USE_PPM 0

static cpl_error_code verify_calibration(const cpl_table *selected,
                                         const cpl_table *linetable, 
                     double TOLERANCE,
                                         double red_chisq,cpl_table* qclog);
static cpl_error_code compute_lambda(cpl_table *linetable, 
                     const polynomial *dispersion_relation, 
                     const polynomial *dispersion_variance,
                                     bool verbose);

static int identify_lines(cpl_table *linetable, 
                          const cpl_table *line_refer, 
                          double ALPHA);

static polynomial *calibrate_global(const cpl_table *linetable,
                                    cpl_table **selected,
                    int degree, bool verbose,
                                    bool reject,
                    double TOLERANCE, 
                                    double kappa,
                    double *red_chisq, 
                    polynomial **dispersion_variance,
                    double *pixelsize,
                    double *rms_wlu,
                    double *rms_pixels);

/*----------------------------------------------------------------------------*/
/**
  @brief    Obtain final dispersion relation
  @param    linetable         The line table containing the already detected
                              emission lines
  @param    line_refer        The wavelength catalogue
  @param    guess_dispersion  The initial dispersion relation in the form 
                              @em lambda * @em m = @em f(@em x, @em m)
  @param    DEGREE            Degree of both independent variables of the 
                              dispersion polynomial.
  @param    TOLERANCE         When making the final fit, lines with residuals 
                              worse than @em TOLERANCE are excluded.
                  If positive, this tolerance is considered in pixel units,
                  otherwise in w.l.u.
  @param    ALPHA             Parameter that controls the next-neighbour distance
                              during line identification. See @c identify_lines()
                              for details.
  @param    MAXERROR          If the RMS of the fit is larger than this number
                              (in pixels), the calibration loop terminates with
                              an error. This is to ensure a graceful exit
                              (when incorrect identifications are made)
  @param    kappa             used in outlier rejectiong
  @return   The obtained initial dispersion solution of the form
            @em lambda * @em m = @em f(@em x, @em m), 
            or NULL on error

  This function performs a wavelength calibration of the detected lines
  listed in the provided @em linetable, starting from the @em guess_dispersion
  solution. Results of the calibration are written to the @em linetable.
  
  The function will iteratively identify as many of the detected lines as
  possible (using the specified wavelength catalogue, see also 
  @c identify_lines() ), then update the fit polynomial 
  (see also @c calibrate_global() ). This loop continues until no new line
  identifications can be made. After this first convergence all
  identifications are reset (to remove possible false identifications), and
  the loop repeats, but this time ignoring lines with residuals worse than
  @em TOLERANCE . The final solution is returned.

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

polynomial *
uves_wavecal_identify(cpl_table *linetable, 
              const cpl_table *line_refer, 
              const polynomial *guess_dispersion, 
              int DEGREE, double TOLERANCE, 
              double ALPHA, double MAXERROR,
                      double kappa,
                      const int trace,const int window,cpl_table* qclog)
{
    polynomial *dispersion_relation = NULL; /* Result */
    polynomial *dispersion_variance = NULL; /* Variance of result, 
                           written to line table */
    int current_id; /* Current and previous number of line identifications */
    int previous_id;
    int idloop;             /* Number of iterations of grand loop */
    int n;                  /* Number of iterations in ID loop */
    double pixelsize;       /* Average conversion factor between pixels and wlu */
    double red_chisq;       /* Reduced chi^2 of fit         */
    cpl_table *selected = NULL;  /* Lines used in final fit */
    char qc_key[40];

    passure( linetable        != NULL, " ");
    passure( line_refer       != NULL, " ");
    passure( guess_dispersion != NULL, " ");

    assure( 0 < ALPHA && ALPHA <= 1, CPL_ERROR_ILLEGAL_INPUT, 
        "Illegal alpha = %e", ALPHA);

    /* Calculate LambdaC from the initial dispersion relation */
    {
    cpl_table_new_column(linetable, LINETAB_LAMBDAC    , CPL_TYPE_DOUBLE);
    cpl_table_new_column(linetable, "dLambdaC"         , CPL_TYPE_DOUBLE);
    cpl_table_new_column(linetable, LINETAB_PIXELSIZE  , CPL_TYPE_DOUBLE);
    cpl_table_new_column(linetable, LINETAB_RESIDUAL   , CPL_TYPE_DOUBLE);
    cpl_table_new_column(linetable, "Residual_pix"     , CPL_TYPE_DOUBLE);
    cpl_table_new_column(linetable, "Lambda_candidate" , CPL_TYPE_DOUBLE);
    cpl_table_new_column(linetable, "dLambda_candidate", CPL_TYPE_DOUBLE);
    cpl_table_new_column(linetable, "dLambda_cat_sq"   , CPL_TYPE_DOUBLE);
    cpl_table_new_column(linetable, "dLambda_nn_sq"    , CPL_TYPE_DOUBLE);

    /* Create columns 'Ident' and 'dIdent' (uncertainty) and fill with
       invalid (no identification made) */
    cpl_table_new_column(linetable, "Ident", CPL_TYPE_DOUBLE);
    cpl_table_new_column(linetable, "dIdent",CPL_TYPE_DOUBLE);
    cpl_table_set_column_invalid(linetable, "Ident", 0, cpl_table_get_nrow(linetable));
    cpl_table_set_column_invalid(linetable, "dIdent",0, cpl_table_get_nrow(linetable));
    
    cpl_table_set_column_unit(linetable,LINETAB_LAMBDAC,"Angstrom" );
    cpl_table_set_column_unit(linetable,"dLambdaC","Angstrom" );
    cpl_table_set_column_unit(linetable,LINETAB_PIXELSIZE,"Angstrom" );
    cpl_table_set_column_unit(linetable,LINETAB_RESIDUAL,"Angstrom" );
    cpl_table_set_column_unit(linetable,"Residual_pix","Angstrom" );
    cpl_table_set_column_unit(linetable,"Lambda_candidate","Angstrom" );
    cpl_table_set_column_unit(linetable,"dLambda_candidate","Angstrom" );
    cpl_table_set_column_unit(linetable,"dLambda_cat_sq","Angstrom" );
    cpl_table_set_column_unit(linetable,"dLambda_nn_sq","Angstrom" );
    cpl_table_set_column_unit(linetable,"Ident","Angstrom" );
    cpl_table_set_column_unit(linetable,"dIdent","Angstrom" );
    /* Residuals are not calculated because 'Ident' is invalid */
    check( compute_lambda(linetable, guess_dispersion, NULL, false), 
           "Error applying dispersion relation");
    }


#if USE_PPM
    for (idloop = 2; idloop <= 2; idloop += 1)
#else
    for (idloop = 1; idloop <= 2; idloop += 1)
#endif
    {

        current_id = 0;
        n = 0;
        /* Iterate until no more identifications can be made */
        do {
        double rms_wlu;
        double rms_pixels;
                bool reject = (idloop == 2);
#if USE_PPM
                int nident_ppm;
#endif
        
        previous_id = current_id;
        n++;
        
        /* Identify lines */
        check( current_id = identify_lines(linetable, line_refer, ALPHA), 
               "Error identifying lines");


#if USE_PPM
                /* Try PPM */
                check( nident_ppm = uves_wavecal_identify_lines_ppm(linetable, line_refer),
                       "Error during point pattern matching");

                cpl_table_erase_column(linetable, "Ident");
                cpl_table_duplicate_column(linetable, "Ident", linetable, "Ident_ppm");
                current_id = nident_ppm;

                /* FIXME: This only works if 'dIdent' is constant.
                   We should propagate error bars during ppm matching */
                cpl_table_fill_column_window(linetable, "dIdent",
                                             0, cpl_table_get_nrow(linetable),
                                             cpl_table_get_column_mean(linetable, "dIdent"));
#endif

        /* Calibrate with 
         * 1st loop: tolerance=infinity (i.e. all identified lines are considered good). 
         * 2nd loop: use specified tolerance (ignore outliers)
         */
        uves_polynomial_delete(&dispersion_relation);
        uves_polynomial_delete(&dispersion_variance);

        check( dispersion_relation = calibrate_global(
               linetable, NULL,
                           DEGREE, false,
                           reject,
               TOLERANCE,
               kappa,
               &red_chisq,
               &dispersion_variance,
               &pixelsize,
               &rms_wlu,
               &rms_pixels),
               "Could not perform global calibration");

        uves_msg_debug("Average pixelsize = %f wlu", pixelsize);
                if (idloop == 1)
                    {
                        uves_msg("%d identifications made. RMS = %.5f wlu = %.3f "
                                 "pixels (no rejection)", 
                                 current_id, rms_wlu, rms_pixels);




                    }
                else
                    {
                        uves_msg("%d identifications made. RMS = %.5f wlu = %.3f "
                                 "pixels (%f %s rejection, kappa = %.1f)", 
                                 current_id, rms_wlu, rms_pixels,
                                 fabs(TOLERANCE), (TOLERANCE > 0) ? "pixels" : "wlu",
                                 kappa);
                    }

                sprintf(qc_key,"QC TRACE%d WIN%d NLINID%d",trace,window,idloop);
                ck0_nomsg(uves_qclog_add_int(qclog,qc_key,current_id,
                                             "ThAr lamp identified lines",
                                             "%d"));

#if USE_PPM
                uves_msg("%d identifications from point pattern matching",
                         nident_ppm);
#endif
        
        assure( rms_pixels < MAXERROR, CPL_ERROR_CONTINUE,
            "Wavelength calibration did not converge. "
            "After %d iterations the RMS was %f pixels. "
            "Try to improve on the initial solution", n, rms_pixels);
        

        /* Apply calibration result */
        check( compute_lambda(linetable, dispersion_relation, dispersion_variance,
                                      false),
               "Error applying dispersion relation");


        }
        while (current_id > previous_id) ;

        sprintf(qc_key,"QC TRACE%d WIN%d NLINID NITERS",trace,window);
        ck0_nomsg(uves_qclog_add_int(qclog,qc_key,idloop+1,
                                     "Number of iterations",
                                     "%d"));



        if (idloop == 1)
        {
            /* 
             * Remove all identifications and repeat
             */
            
            uves_msg("Identification loop converged. Resetting identifications");
            cpl_table_set_column_invalid(linetable, "Ident", 0, 
                         cpl_table_get_nrow(linetable));
        }
    }

    /* Calibrate again with a global polynomial, but this time don't
       use lines with residuals worse than TOLERANCE */
    uves_polynomial_delete(&dispersion_relation);
    uves_polynomial_delete(&dispersion_variance);
    uves_free_table(&selected);
    
    check( dispersion_relation = calibrate_global(linetable,
                                                  &selected,
                                                  DEGREE, true,
                                                  true,  /* do rejection? */
                                                  TOLERANCE,
                                                  kappa,
                                                  &red_chisq,
                                                  &dispersion_variance,
                                                  NULL, NULL, NULL),
           "Could not perform global calibration");
    
    /* Update the computed wavelengths */
    check( compute_lambda(linetable, dispersion_relation, dispersion_variance,
                          true), 
           "Error applying dispersion relation");
    
    /* Add columns 'Select' and 'NLinSol' to linetable.
       The columns defines which lines were identified,
       and which lines were used in the final fit */
    {
        int i, j;

        /* Tables are sorted by Order, X */

        cpl_table_new_column(linetable, "NLinSol", CPL_TYPE_INT);
        cpl_table_new_column(linetable, "Select", CPL_TYPE_INT);

        cpl_table_fill_column_window_int(linetable, "NLinSol", 
                                         0, cpl_table_get_nrow(linetable),
                                         0);
        cpl_table_fill_column_window_int(linetable, "Select", 
                                         0, cpl_table_get_nrow(linetable),
                                         0);

        j = 0;
        for (i = 0; i < cpl_table_get_nrow(selected); i++) {
            int order = cpl_table_get_int(selected, "Order", i, NULL);
            double  x = cpl_table_get_double(selected, "X", i, NULL);
            int order2;
            double x2;

            /* Find this line in the original linetable */
            passure( j < cpl_table_get_nrow(linetable), "%d %" CPL_SIZE_FORMAT "",
                     j, cpl_table_get_nrow(linetable));
            do {
                order2 = cpl_table_get_int(linetable, "Order", j, NULL);
                x2     = cpl_table_get_double(linetable, "X", j, NULL);
                if (cpl_table_is_valid(linetable, "Ident", j))
                    {
                        cpl_table_set_int(linetable, "Select", j, 1);
                    }
                j++;

            } while (order2 < order || x2 < x - 0.1);
            
            passure( order2 == order && fabs(x2 - x) < 0.1,
                     "%d %d %g %g", order2, order, x2, x);
            
            cpl_table_set_int(linetable, "NLinSol", j-1, 1);
        }
    }

    /* Display results */
    check( verify_calibration(selected, linetable, TOLERANCE, red_chisq,qclog), 
       "Error verifying calibration");
    
  cleanup:
    uves_free_table(&selected);
    uves_polynomial_delete(&dispersion_variance);
    return dispersion_relation;
}

/*----------------------------------------------------------------------------*/
/**
   @brief    Report quality of calibration
   @param    linetable     The line table
   @param    TOLERANCE     When reporting the RMS of the fit, exclude 
                           lines with residuals worse than @em TOLERANCE from 
               the computation. If positive, this tolerance is
               considered in pixel units, otherwise in w.l.u.
   @param    red_chisq     Reduced chi^2 of the calibration
   @return   CPL_ERROR_NONE iff okay

   This function reports the RMS of a wavelength calibration, as well as 
   the percentage of the brightest lines that were identified.
**/
/*----------------------------------------------------------------------------*/
static cpl_error_code
verify_calibration(const cpl_table *selected,
                   const cpl_table *linetable, double TOLERANCE,
                   double red_chisq, cpl_table* qclog)
{
    cpl_table *brightest  = NULL;
    double median_intensity;
    int ninvalid;    /* Number of unidentified lines among the brightest half */
    double ratio;
    double rms_wlu;
    double rms_pixels;
    double rms_speed;
    char qc_key[40];
    
    {
    double mean;
    double stdev;
    
    check(( mean = cpl_table_get_column_mean (selected, LINETAB_RESIDUAL),
        stdev= cpl_table_get_column_stdev(selected, LINETAB_RESIDUAL),
        rms_wlu = sqrt(mean*mean + stdev*stdev),
        
        mean = cpl_table_get_column_mean (selected, "Residual_pix"),
        stdev= cpl_table_get_column_stdev(selected, "Residual_pix"),
        rms_pixels = sqrt(mean*mean + stdev*stdev)),
          "Error reading RMS of fit");
    }
    rms_speed=rms_wlu * SPEED_OF_LIGHT/
       cpl_table_get_column_mean(selected,LINETAB_LAMBDAC);
    uves_msg("%" CPL_SIZE_FORMAT " lines accepted", cpl_table_get_nrow(selected));
    uves_msg("Average RMS of calibration (tolerance = %.3f %s) = %.5f wlu = %.4f pixels ~ %.1f m/s",
         fabs(TOLERANCE),
         (TOLERANCE > 0) ? "pixels" : "wlu",
         rms_wlu, rms_pixels, rms_speed);
    
   sprintf(qc_key,"QC LINE RESIDRMS WLU");
   ck0_nomsg(uves_qclog_add_double(qclog,qc_key,rms_wlu,
                                   "Line ID RMS TRACE0 WIN2 [Angstrom]",
                                   "%f"));
   sprintf(qc_key,"QC LINE RESIDRMS PIX");
   ck0_nomsg(uves_qclog_add_double(qclog,qc_key,rms_pixels,
                                   "Line ID RMS TRACE0 WIN2 [pix]",
                                   "%f"));
   sprintf(qc_key,"QC LINE RESIDRMS SPEED");
   ck0_nomsg(uves_qclog_add_double(qclog,qc_key,rms_speed,
                                   "Line ID RMS TRACE0 WIN2 [m/s]",
                                   "%f"));


    uves_msg("Reduced chi^2 of calibration = %f", red_chisq);
    sprintf(qc_key,"QC LINE IDCHI2");
    ck0_nomsg(uves_qclog_add_double(qclog,qc_key,red_chisq,
                                   "Reduced chi^2 of line ID TRACE0 WIN2",
                                   "%f"));

    if (red_chisq < .01)
    {
        uves_msg_warning("Reduced chi^2 of fit is less than 1/100: %f", 
                 red_chisq);
    }
    if (red_chisq > 100)
    {
        uves_msg_warning("Reduced chi^2 of fit is greater than 100: %f", 
                 red_chisq);
    }
    
    check(( median_intensity = cpl_table_get_column_median(linetable, "Peak"),
        brightest = uves_extract_table_rows(linetable, "Peak", 
                        CPL_GREATER_THAN, 
                        median_intensity),
        ninvalid = cpl_table_count_invalid(brightest, "Ident")),
      "Error counting identifications");

    ratio = 1 - ((double) ninvalid)/cpl_table_get_nrow(brightest);
    uves_msg("Percentage of identifications among the half brighter lines : %.2f %%",
         100*ratio);

    sprintf(qc_key,"QC LINE HALFBRIG");
    ck0_nomsg(uves_qclog_add_double(qclog,qc_key,100*ratio,
                                   "Half brighter lines frac TRACE0 WIN2",
                                   "%f"));

  cleanup:
    uves_free_table(&brightest);

    return cpl_error_get_code();
}

/*----------------------------------------------------------------------------*/
/**
   @brief    Apply dispersion relation to line table
   @param    linetable            The line table
   @param    dispersion_relation  The dispersion relation
   @param    dispersion_variance  The variance of the dispersion relation
   @param    verbose              print warning if dl/dx is negative?
   @return   CPL_ERROR_NONE iff okay

   This function applies the provided dispersion relation to the line table, 
   thereby calculating the predicted wavelengths for lines in the spectrum, 
   the residual of the fit, the local pixelsize in w.l.u., and if 
   @em dispersion_variance is non-NULL, the uncertainty of the fitted wavelength.
**/
/*----------------------------------------------------------------------------*/
static cpl_error_code
compute_lambda(cpl_table *linetable, 
           const polynomial *dispersion_relation, 
           const polynomial *dispersion_variance,
               bool verbose)
{
    int i;
    bool printed_warning = false;
    
    /* Check input */
    passure(linetable           != NULL, " ");
    passure(dispersion_relation != NULL, " ");
    /* 'dispersion_variance' may be NULL */
    
    passure( uves_polynomial_get_dimension(dispersion_relation) == 2, "%d", 
         uves_polynomial_get_dimension(dispersion_relation));
    
    /* Input columns */
    passure(cpl_table_has_column(linetable, "X")           , " ");
    passure(cpl_table_has_column(linetable, "Order")       , " ");
    passure(cpl_table_has_column(linetable, "Ident")       , " ");
    /* Output columns */
    passure(cpl_table_has_column(linetable, LINETAB_LAMBDAC)     , " ");
    /* The column 'dLambdaC' is set to invalid if 'dispersion_variance' is NULL */
    passure(cpl_table_has_column(linetable, "dLambdaC")    , " ");  
    passure(cpl_table_has_column(linetable, "dIdent")      , " ");
    passure(cpl_table_has_column(linetable, LINETAB_RESIDUAL), " ");
    passure(cpl_table_has_column(linetable, "Residual_pix"), " ");
    passure(cpl_table_has_column(linetable, LINETAB_PIXELSIZE)   , " ");
    
    /* The linetable is sorted w.r.t. order. 
       Move to the first order above minorder */
    for(i = 0; i < cpl_table_get_nrow(linetable); i++)
    {
        int order;
        double x, dfdx;
        double lambdac, dlambdac, pixelsize;
        order = cpl_table_get_int(linetable, "Order", i, NULL);
        
        x     = cpl_table_get_double(linetable, "X", i, NULL);
        
        /* Evaluate the dispersion relation
           m.lambda = f(x,m)  (2d global fit)  */
        
        lambdac =
        uves_polynomial_evaluate_2d(dispersion_relation, x, order) / order;
        
        /* Pixelsize = dl/dx = (df/dx)/m  (for fixed m) */
            dfdx = uves_polynomial_derivative_2d(dispersion_relation, x, order, 1);
            if (dfdx < 0) {
                if (!printed_warning && verbose) {
                    uves_msg_warning("Inferred dispersion (dlambda/dx) is negative at"
                                     "(x, order) = (%f, %d)", x, order);
                    printed_warning = true;  /* To avoid repeating the same warning */
                }
                else {
                    uves_msg_debug("Inferred dispersion (dlambda/dx) is negative at "
                                   "(x, order) = (%f, %d)", x, order);
                }
            }
            pixelsize = dfdx / order;
        
        check(( cpl_table_set_double(linetable, LINETAB_LAMBDAC , i, lambdac),
            cpl_table_set_double(linetable, LINETAB_PIXELSIZE, i, pixelsize)),
            "Error writing table");
        
        if (dispersion_variance != NULL)
        {
            /* d( lambda  (x, order) ) = 
               d( lambda*m(x, order) ) / m    */
            dlambdac = 
            sqrt(uves_polynomial_evaluate_2d(dispersion_variance, x, order))
            / order;
            
            cpl_table_set_double(linetable, "dLambdaC" , i, dlambdac);
        }
        else
        {
            /* Only the ratio of a line's "dLambdaC" to other
               lines' are used, so set "dLambdaC" to a constant value
               when the actual uncertainty is not known
            */
            cpl_table_set_double(linetable, "dLambdaC" , i, 1.0);
        }
        
        /* If line is identified, calculate residual */
        if (cpl_table_is_valid(linetable, "Ident", i)) 
        {
            double ident = cpl_table_get_double(linetable, "Ident", i, NULL);
            cpl_table_set_double(linetable, LINETAB_RESIDUAL, i,
                     ident - lambdac);
            cpl_table_set_double(linetable, "Residual_pix", i, 
                     (ident - lambdac)/pixelsize);
        }
        else
        {
            cpl_table_set_invalid(linetable, LINETAB_RESIDUAL, i);
            cpl_table_set_invalid(linetable, "Residual_pix", i);
        }
    }
    
    /* Sort by 'Order' (ascending), then 'X' (ascending) */
    check( uves_sort_table_2(linetable, "Order", "X", false, false), 
       "Error sorting table");
    
  cleanup:
    return cpl_error_get_code();
}


/*----------------------------------------------------------------------------*/
/**
   @brief    Identify lines by comparing to catalogue wavelengths
   @param    linetable     The line table containing the line positions and
                           predicted wavelengths
   @param    line_refer    The wavelength catalogue
   @param    ALPHA         Parameter < 1 to control distance to nearest neighbour
   @return   The total number of lines identified (including previous identifications),
             or undefined on error.

   This function identifies lines in the provided @em linetable by comparing
   to the reference table (the wavelength catalogue). The identified 
   wavelengths are written to the "Ident" column of the line table.

   An identification is made iff
   
   - The catalogue wavelength is within two linewidths of the computed
   (predicted) wavelength: 
   | @em lambda_cat - @em lambda_com | < 2 * @em sigma, 
   where @em sigma is the detected line width,
   
   - The nearest neighbour (in the spectrum and in the catalogue) is farther away
   than the candidate catalogue wavelength (after multiplying by the
   "safety parameter", @em ALPHA):
   | @em lambda_cat - @em lambda_nn| * ALPHA > | @em lambda_cat - @em lambda_com |,
   
   - The nearest neighbour (in the spectrum and in the catalogue) is farther away
   than the average tolerance distance, which measures the precision of the identifications:
   @em tolerance < ALPHA * | @em lambda_cat - @em lambda_nn| . See code for the exact
   definition of the @em tolerance .

   The purpose of the first criterion is to make the correct identifications.
   The purpose of the latter two criterions is to avoid making incorrect identifications.

   If a line was previously identified (implied by a valid @em "Ident" 
   column value) but now fails to meet the ID criterium, it is not deleted.
**/
/*----------------------------------------------------------------------------*/

static int
identify_lines(cpl_table *linetable, const cpl_table *line_refer, double ALPHA)
{
    int number_identifications = 0;      /* Result */
    int linetable_size;
    int linerefer_size;
    int row;
    int *histogram = NULL;
    const double minlog  = -5.0;         /* Histogram (it's sort of ugly
                        to hardcode these numbers, but
                        as long as it works, ...) */
    const double maxlog  = 15.0;
    const int nbins       = 400;
    double error = 0;                    /* Dimensionless factor
                        that controls IDs */
    double average_dlambda_com = 0;      /* Average of uncertainty of 
                        predicted wavelenghts */

    /* Check input */
    passure( linetable  != NULL, " ");
    /* Line table input columns */
    passure( cpl_table_has_column(linetable, LINETAB_LAMBDAC  ), " "); /* Predicted
                                      wavelength  */
    passure( cpl_table_has_column(linetable, "dLambdaC" ), " "); /* Predicted wavelength 
                                    uncertainty  */
    passure( cpl_table_has_column(linetable, "X"        ), " "); /* Line position, used
                                    only for messaging */
    passure( cpl_table_has_column(linetable, "Order"    ), " "); /* Absolute order number 
                                    of line */
    passure( cpl_table_has_column(linetable, "Xwidth"   ), " "); /* Line width (sigma) */
    passure( cpl_table_has_column(linetable, LINETAB_PIXELSIZE), " "); /* Pixelsize */

    /* Line table output columns */
    passure( cpl_table_has_column(linetable, "Ident"    ), " "); /* Identified catalogue 
                                    wavelength */
    passure( cpl_table_has_column(linetable, "dIdent"   ), " "); /* Uncertainty of IDed
                                    catalogue wavelength */

    /* Catalogue */
    passure( line_refer != NULL, " ");
    passure( cpl_table_has_column(line_refer, "Wave" ), " ");    /* Catalogue wavelength */
    passure( cpl_table_has_column(line_refer, "dWave"), " ");    /* Uncertainty of
                                    catalogue wavelength */
    
    linetable_size = cpl_table_get_nrow(linetable);
    linerefer_size = cpl_table_get_nrow(line_refer);
    assure(linerefer_size >= 1, CPL_ERROR_ILLEGAL_INPUT, "Empty line reference table");
    
    /* Parameter */
    passure( 0 < ALPHA && ALPHA <= 1, "%e", ALPHA);

    /* Get average uncertainty of predicted wavelength */
    average_dlambda_com = cpl_table_get_column_median(linetable, "dLambdaC");

    /* Initialize histogram to zero */
    histogram = cpl_calloc(nbins, sizeof(int));
    assure_mem( histogram );
    

    /* First: Find distance to closest catalogue match, 
       distance to nearest neighbour, 
       and calculate histogram (to get average of distances to nearest neighbour) */
    for (row = 0; row < linetable_size; row++) {
    double lambda_com;                 /* Computed (predicted) wavelength */
    double line_width;                 /* Line width (sigma) in wlu       */
    double line_fwhm;                  /* Line FWHM in wlu                */
    int order;                         /* (Absolute) order of detected wavelength */
    double lambda_cat;                 /* Catalogue wavelength */
    double lambda_cat_sigma;           /* Catalogue wavelength uncertainty */
    double distance_cat_sq;            /* Distance to catalogue wavelength (squared) */
    double nn_distance_sq;             /* Distance to nearest neighbour (squared) */
    int row_cat;                       /* Row number of best matching catalogue wavelength */
    
    /* Read line table */
    lambda_com  = cpl_table_get_double(linetable, LINETAB_LAMBDAC   , row, NULL);
    order       = cpl_table_get_int   (linetable, "Order"     , row, NULL);

    
    line_width = 
        cpl_table_get_double(linetable, "Xwidth"    , row, NULL) *
        fabs(cpl_table_get_double(linetable, LINETAB_PIXELSIZE , row, NULL)); 
    /* Convert pixel->wlu */

    line_fwhm = TWOSQRT2LN2 * line_width;
    
    /* Find closest match in catalogue */
    row_cat          = uves_wavecal_find_nearest(
        line_refer, lambda_com, 0, linerefer_size - 1);
    lambda_cat       = cpl_table_get_double(line_refer, "Wave", row_cat, NULL);
    lambda_cat_sigma = cpl_table_get_double(line_refer, "dWave",row_cat, NULL);

    /* Distance to closest match */
    distance_cat_sq = (lambda_com - lambda_cat)*(lambda_com - lambda_cat);
    
        /* Determine the distance to the next neighbour
     * There are (max) 4 candiates: 2 neigbours in spectrum (i.e. line table)
     *                          and 2 neigbours in line catalogue
     */
    {


        double lambda_cat_prev, lambda_cat_next;

        nn_distance_sq = DBL_MAX;

        /* Read previous and next rows of line table */
        if (row >= 1) 
        {
            int order_prev      = cpl_table_get_int   (
            linetable, "Order"  , row - 1, NULL);
            double lambda_com_prev = cpl_table_get_double(
            linetable, LINETAB_LAMBDAC, row - 1, NULL);
            
            if (order == order_prev) 
            {
                nn_distance_sq = uves_min_double(nn_distance_sq,
                                 (lambda_com_prev - lambda_com)*
                                 (lambda_com_prev - lambda_com)
                );
            }
        }

        if (row <= linetable_size - 2) 
        {
            int order_next      = cpl_table_get_int   (linetable, "Order",
                               row + 1, NULL);
            double lambda_com_next = cpl_table_get_double(linetable, LINETAB_LAMBDAC,
                               row + 1, NULL);
            
            if (order == order_next) 
            {
                nn_distance_sq = uves_min_double(nn_distance_sq,
                                 (lambda_com_next - lambda_com)*
                                 (lambda_com_next - lambda_com)
                );
            }
        }
        
        /* Read previous and next rows of catalogue */
        if (row_cat >= 1)
        {
            lambda_cat_prev = cpl_table_get_double(
            line_refer, "Wave", row_cat - 1, NULL);

            nn_distance_sq = uves_min_double(
            nn_distance_sq,
            (lambda_cat_prev - lambda_cat)*
            (lambda_cat_prev - lambda_cat)
            );
        }
        if (row_cat <= linerefer_size - 2) 
        {
            lambda_cat_next = cpl_table_get_double(
            line_refer, "Wave", row_cat + 1, NULL);

            nn_distance_sq = uves_min_double(
            nn_distance_sq,
            (lambda_cat_next - lambda_cat)*
            (lambda_cat_next - lambda_cat)
            );
        }

        /* Update distance to nearest neighbour with a 
           safety margin (determined by parameter ALPHA < 1) */
        if (nn_distance_sq < DBL_MAX)
        {
            nn_distance_sq *= ALPHA*ALPHA;
        }
        
    }/* Find next neighbour */
    
    /* Update line table */
    cpl_table_set_double(linetable, "Lambda_candidate", row, lambda_cat);
    cpl_table_set_double(linetable, "dLambda_candidate",row, lambda_cat_sigma);
    cpl_table_set_double(linetable, "dLambda_cat_sq", row, distance_cat_sq);
    cpl_table_set_double(linetable, "dLambda_nn_sq", row, nn_distance_sq);

    /* Update histogram with the interval
       [distance_cat_sq ; nn_distance_sq]  (in units of line_fwhm) */
    {
        int ilow  = uves_round_double((0.5*log(distance_cat_sq/(line_fwhm*line_fwhm))
                       - minlog)/(maxlog - minlog) * nbins);
        int ihigh = uves_round_double((0.5*log(nn_distance_sq /(line_fwhm*line_fwhm))
                       - minlog)/(maxlog - minlog) * nbins);
        int i;
        
        for (i = uves_max_int(ilow, 0); i < uves_min_int(ihigh, nbins); i++) 
        {
            histogram[i] += 1;
        }
    }
    }/* ... finding neighbours */
    
    /* Determine error as peak of histogram */
    {
    int i;
    int maxfreq = -1;
    for (i = 0; i < nbins; i++) 
        {
        uves_msg_debug("histogram[%d] = %d", i, histogram[i]);
        if (histogram[i] > maxfreq) 
            {
            maxfreq = histogram[i];
            error   = exp( i / ((double)nbins) * (maxlog - minlog) + minlog ) ;
            /* == the dimensionless factor to be multiplied by Xwidth */
            }
        }
    uves_msg_debug("Dimensionless error factor is %f", error);
    }
    
    /* Sketch of situation:
       
  lambda_com                  Nearest neighbour

      |                            |
      |    |                       |
      |    |                       |
      |    |                       |
           |

    lambda_cat


     The 'average' (as inferred from the histogram)
     midpoint between 'lambda_cat' and 'nearest neighbour'
     is at   'error' * 'line_fwhm' .
    */
    
    /* Make the identification if
       
    1) the catalogue candidate is within two sigma:
         | lambda_cat - lambda_com | < 2 * dlambda_com

    and

    2) after multiplying the distance to the nearest neighbour by ALPHA < 1,
    the nearest neighbour is farther away than the catalogue wavelength 
         distance_nn  >  distance_cat
    and farther away than the tolerance
         distance_nn  >  line_fwhm * error
     
    */
    for (row = 0; row < linetable_size; row++)
    {
        double distance_cat_sq;              /* Distance to catalogue wavelength (squared) */
        double nn_distance_sq;               /* Distance to nearest neighbour (squared) */
        double tolerance_sq;
        double dlambda_com;
        double line_width;                   /* Line width (1 sigma) */
        double line_fwhm;
        double lambda_cat;
        double lambda_cat_sigma;             /* Uncertainty of lambda_cat */
        
        lambda_cat       = cpl_table_get_double(linetable,  "Lambda_candidate", row, NULL);
        lambda_cat_sigma = cpl_table_get_double(linetable, "dLambda_candidate", row, NULL);
        
        
        /* Sigma less than 1 pixel is usually not
           justified by the data (which obviously 
           has a resolution of only 1 pixel). Such
           an underenstimation of the uncertainty
           leads to wrong identifications.
           Therefore use a width of at least 1 pixel */
        line_width =
        uves_max_double(1, cpl_table_get_double(linetable, "Xwidth"    , row, NULL)) *
        fabs(cpl_table_get_double(linetable, LINETAB_PIXELSIZE , row, NULL));
        /* convert to wlu */
        
        line_fwhm = TWOSQRT2LN2 * line_width;

        /* As the uncertainty of the computed wavelength is used
         *  line_fwhm (in w.l.u.)
         * To take into account the fact that lines near the edge of
         * the chip have larger error of the computed wavelength,
         * this is also scaled according to the accuracy of the dispersion
         * relation, i.e. multiplied by  dl/<dl>,
         * where <dl> is an average, say the median, of uncertainties of
         * all predicted wavelengths.
         */
        
        dlambda_com = line_fwhm 
        * cpl_table_get_double(linetable, "dLambdaC"  , row, NULL)
        / average_dlambda_com;
        
        tolerance_sq = line_fwhm*line_fwhm * error*error;
        
        distance_cat_sq = cpl_table_get_double(linetable, "dLambda_cat_sq", row, NULL);
        nn_distance_sq  = cpl_table_get_double(linetable, "dLambda_nn_sq" , row, NULL);
        
#if WANT_BIG_LOGFILE
        uves_msg_debug("(order,x) = (%d,%f) lcom = %f+-%f lcat = %f "
               "dist_cat = %f (%f pixels) tolerance = %.3f error = %f "
               "nn = %f (%f pixels)", 
               cpl_table_get_int   (linetable, "Order"  , row, NULL),
               cpl_table_get_double(linetable, "X"      , row, NULL),
               cpl_table_get_double(linetable, LINETAB_LAMBDAC, row, NULL),
               dlambda_com,
               lambda_cat,
               sqrt(distance_cat_sq),
               sqrt(distance_cat_sq)
               /cpl_table_get_double(linetable, LINETAB_PIXELSIZE, row, NULL),
               sqrt(tolerance_sq),
               error,
               sqrt(nn_distance_sq),
               sqrt(nn_distance_sq)
               /cpl_table_get_double(linetable, LINETAB_PIXELSIZE, row, NULL));
#endif
        
        /* Make the ID? */
        if (distance_cat_sq < (dlambda_com)*(dlambda_com)
        && tolerance_sq < nn_distance_sq
        && distance_cat_sq < nn_distance_sq)
        {
            number_identifications++;
            cpl_table_set_double(linetable, "Ident", row, lambda_cat);
            cpl_table_set_double(linetable, "dIdent",row, lambda_cat_sigma);
#if WANT_BIG_LOGFILE
            uves_msg_debug("ID made");
#endif
        }
        else 
        {
            if (cpl_table_is_valid(linetable, "Ident", row)) {
            number_identifications++;                      
            /* Also count lines that were already identified */
            uves_msg_debug("Line at (%d,%f) does not match ID criterion anymore",
                       cpl_table_get_int   (linetable, "Order", row, NULL),
                       cpl_table_get_double(linetable, "X", row, NULL)
            );
        }
        }
    }

  cleanup:
    cpl_free(histogram);
    return number_identifications;
}

/*----------------------------------------------------------------------------*/
/**
   @brief    Create a fit of all orders
   @param    linetable     The line table
   @param    selected      (output) if non-NULL, subset of linetable containing
                           the lines which were used in the final fit
   @param    degree        Degree of both independent variables of polynomial fit
   @param    verbose       Be verbose about autodegree fitting?
   @param    reject        Do rejection?
   @param    TOLERANCE     Before fitting, exclude lines with residuals worse than
                           @em TOLERANCE. If positive, this tolerance
                           is considered in pixel units, otherwise in w.l.u.

   @param    kappa         used for removing outliers
   @param    red_chisq            If non-NULL, the reduced chi square of the fit.
   @param    dispersion_variance  If non-NULL, the variance of the fit returned polynomial.
   @param    pixelsize            If non-NULL, the average of d(lambda)/dx
   @param    rms_wlu              If non-NULL, the root-mean-square residual (w.l.u)
   @param    rms_pixels           If non-NULL, the root-mean-square residual (pixels)
   @return   The obtained dispersion relation in the form
             @em lambda * @em m = @em f(@em x, @em m), 
         or NULL on error

   @note Un-identified lines and lines with residuals larger then @em TOLERANCE
   (from the previous fit) are excluded from the fit.
**/
/*----------------------------------------------------------------------------*/
static polynomial *
calibrate_global(const cpl_table *linetable,
                 cpl_table **selected,
         int degree, bool verbose,
                 bool reject,
         double TOLERANCE,
                 double kappa,
         double *red_chisq, polynomial **dispersion_variance,
         double *pixelsize,
         double *rms_wlu,
         double *rms_pixels)
{
    polynomial *dispersion_relation = NULL; /* Result */
    cpl_table *identified = NULL;
    int valid_ids = 
    cpl_table_get_nrow(linetable) - 
    cpl_table_count_invalid(linetable, "Ident");
    int rejected;
    
    passure( (pixelsize == NULL) == (rms_wlu    == NULL) &&
         (pixelsize == NULL) == (rms_pixels == NULL), " ");

    assure( degree < 0 ||
            valid_ids >= (degree + 1)*(degree + 1), CPL_ERROR_ILLEGAL_INPUT,
        "There are not enough identifications to create a %d.-degree global fit. "
        "%d needed. %d found", degree, (degree + 1)*(degree + 1), valid_ids);
    
    identified = cpl_table_duplicate(linetable);
    assure_mem(identified);

    /* Delete rows with invalid 'Ident' and large residuals */
    if (reject)
        {
            check_nomsg( rejected = uves_delete_bad_lines(identified, TOLERANCE, kappa) );
            uves_msg_debug("%d lines rejected %f %f", rejected, TOLERANCE, kappa);
        }
    else
        {
            check( uves_erase_invalid_table_rows(identified, "Ident"),
                   "Error erasing un-identified lines");
        }

    
    /* Create column 'Aux' = 'Order' * 'Ident' */
    check((  cpl_table_duplicate_column(identified, "Aux", identified, "Ident"),
             cpl_table_multiply_columns(identified, "Aux", "Order"),
             
             /* Create column 'dAux' = 'Order' * 'dIdent' */
             cpl_table_duplicate_column(identified, "dAux", identified, "dIdent"),
             cpl_table_multiply_columns(identified, "dAux", "Order")),
          "Error setting up temporary table");

    /* Fit */
    
    if (degree >= 0) {
        check( dispersion_relation =
               uves_polynomial_regression_2d(identified, 
                                             "X", "Order", "Aux", 
                                             "dAux", /* Use "dAux" for weighting,
                                                        to be able to compute an uncertainty
                                                        of WAVEC.
                                                        
                                                        It would probably make more sense
                                                        to use the uncertainty of 'dX' for
                                                        weighting. */
                                             degree, degree,
                                             NULL, NULL, NULL,     /* Don't add extra columns */
                                             NULL,                 /* mse */
                                             red_chisq,
                                             dispersion_variance, 
                                             reject ? kappa : -1, -1),
               "Error fitting polynomial. Possible cause: too few (%d) "
               "line identifications", valid_ids);
    } 
    else {
        int max_degree = 8;
        double min_rms = -1; /* disabled */
        double min_reject = -1; /* disabled */
        check( dispersion_relation =
               uves_polynomial_regression_2d_autodegree(identified,
                                                        "X", "Order", "Aux", 
                                                        "dAux", 
                                                        NULL, NULL, NULL,  
                                                        NULL, 
                                                        red_chisq,
                                                        dispersion_variance,
                                                        reject ? kappa : -1,
                                                        max_degree, max_degree, 
                                                        min_rms, min_reject,
                                                        verbose,
                                                        NULL, NULL, 0, NULL),
               "Error fitting polynomial. Possible cause: too few (%d) "
               "line identifications", valid_ids);
    }

    if (pixelsize != NULL)
    {
        /* Compute parameters if requested */

        check( compute_lambda(identified, dispersion_relation, NULL,
                                  false),
           "Error applying dispersion relation");
        
        *pixelsize = cpl_table_get_column_median(identified, LINETAB_PIXELSIZE);
        *rms_wlu   = cpl_table_get_column_stdev (identified, LINETAB_RESIDUAL);
        *rms_pixels= cpl_table_get_column_stdev (identified, "Residual_pix");
    }

    if (selected != NULL) {
        *selected = cpl_table_duplicate(identified);
    }

  cleanup:
    uves_free_table(&identified);
    if (cpl_error_get_code() != CPL_ERROR_NONE)
    {
        uves_polynomial_delete(&dispersion_relation);
    }
    
    return dispersion_relation;
}



/*----------------------------------------------------------------------------*/
/**
   @brief    Identify lines using point pattern matching
   @param    linetable     The line table containing the line positions 
   @param    line_refer    The wavelength catalogue
   @return   number of identifications
**/
/*----------------------------------------------------------------------------*/

int
uves_wavecal_identify_lines_ppm(cpl_table *linetable, const cpl_table *line_refer)
{
    int result = 0;
    int minorder, maxorder;
    int order;
    cpl_table *lt_order = NULL;
    cpl_table *refer_order = NULL;
    cpl_vector *peaks = NULL;
    cpl_vector *lines = NULL;
    cpl_bivector *ids = NULL;

    assure( cpl_table_has_column(linetable, LINETAB_LAMBDAC), CPL_ERROR_DATA_NOT_FOUND,
            "Missing column %s", LINETAB_LAMBDAC);

    assure( cpl_table_has_column(linetable, LINETAB_PIXELSIZE), CPL_ERROR_DATA_NOT_FOUND,
            "Missing column %s", LINETAB_PIXELSIZE);

    assure( cpl_table_has_column(linetable, "Order"), CPL_ERROR_DATA_NOT_FOUND,
            "Missing column %s", "Order");

    minorder = uves_round_double( cpl_table_get_column_min(linetable, "Order"));
    maxorder = uves_round_double( cpl_table_get_column_max(linetable, "Order"));

    /* Reset identifications */
    if (cpl_table_has_column(linetable, "Ident_ppm"))
        {
            cpl_table_erase_column(linetable, "Ident_ppm");
        }

    cpl_table_new_column(linetable, "Ident_ppm", CPL_TYPE_DOUBLE);
    
    for (order = minorder; order <= maxorder; order++)
        {
            const double tolerance = 0.05; /* relative tolerance on interval ratios */
            double min_lambda, max_lambda;
            double min_disp, max_disp;

            /* Extract current order */
           
            uves_free_table(&lt_order);
            lt_order = uves_extract_table_rows(linetable, "Order",
                                               CPL_EQUAL_TO, order); /* Uses integer comparison */

            check_nomsg((min_lambda = cpl_table_get_column_min(lt_order, LINETAB_LAMBDAC),
                         max_lambda = cpl_table_get_column_max(lt_order, LINETAB_LAMBDAC),
                         min_disp   = cpl_table_get_column_min(lt_order, LINETAB_PIXELSIZE)*0.99,
                         max_disp   = cpl_table_get_column_max(lt_order, LINETAB_PIXELSIZE)*1.01));
                        
            uves_free_table(&refer_order);
            refer_order = uves_extract_table_rows(line_refer, "Wave", CPL_GREATER_THAN,
                                                  min_lambda);
            uves_extract_table_rows_local(refer_order, "Wave", CPL_LESS_THAN,
                                          max_lambda);

            /* Convert to vectors */
            {
                int i;
                uves_free_vector(&peaks);
                peaks = cpl_vector_new(cpl_table_get_nrow(lt_order));
                for (i = 0; i < cpl_vector_get_size(peaks); i++)
                    {
                        cpl_vector_set(peaks, i, cpl_table_get_double(lt_order, "X", i, NULL));
                    }
                
                uves_free_vector(&lines);
                lines = cpl_vector_new(cpl_table_get_nrow(refer_order));
                for (i = 0; i < cpl_vector_get_size(lines); i++)
                    {
                        cpl_vector_set(lines, i, cpl_table_get_double(refer_order, "Wave", i, NULL));
                    }
            }
            
            /* Not sure if this is necessary for the PPM algorithm */
            cpl_vector_sort(peaks, 1);
            cpl_vector_sort(lines, 1);

            uves_msg_debug("Call ppm with %" CPL_SIZE_FORMAT " peaks, %" CPL_SIZE_FORMAT " lines, dispersion range = %f - %f A/pixel",
                           cpl_vector_get_size(peaks), 
                           cpl_vector_get_size(lines),
                           min_disp, max_disp);

            uves_free_bivector(&ids);

            ids = cpl_ppm_match_positions(peaks, lines,
                                          min_disp, max_disp,
                                          tolerance, 
                                          NULL, NULL);


            if (ids == NULL)
                {
                    uves_msg_warning("Order %d: Point pattern matching failed", order);
                    if (cpl_error_get_code() != CPL_ERROR_NONE)
                        {
                            uves_msg_debug("%s at %s", cpl_error_get_message(),
                                           cpl_error_get_where());
                            uves_error_reset();
                        }
                }
            else
                {
                    int i, j;

                    uves_msg_debug("%" CPL_SIZE_FORMAT " identifications from point pattern matching (order %d)",
                                   cpl_bivector_get_size(ids), order);

                    result += cpl_bivector_get_size(ids);

                    for (i = 0; i < cpl_table_get_nrow(linetable); i++) {

                        if (cpl_table_get_int(linetable, "Order", i, NULL) == order)
                            for (j = 0; j < cpl_bivector_get_size(ids); j++)
                                {
                                    if (fabs(cpl_table_get_double(linetable, "X", i, NULL) -
                                             cpl_bivector_get_x_data(ids)[j]) < 0.001)
                                        cpl_table_set_double(linetable, "Ident_ppm", i,
                                                             cpl_bivector_get_y_data(ids)[j]);
                                }
                    }
                }
        }
    
  cleanup:
    uves_free_table(&lt_order);
    uves_free_table(&refer_order);
    uves_free_vector(&peaks);
    uves_free_vector(&lines);
    uves_free_bivector(&ids);

    return result;
}
/**@}*/