/*                                                                              *
 *   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: 2010-09-24 09:32:09 $
 * $Revision: 1.23 $
 * $Name: not supported by cvs2svn $
 * $Log: not supported by cvs2svn $
 * Revision 1.21  2007/08/21 13:08:26  jmlarsen
 * Removed irplib_access module, largely deprecated by CPL-4
 *
 * Revision 1.20  2007/06/06 08:17:33  amodigli
 * replace tab with 4 spaces
 *
 * Revision 1.19  2007/05/25 07:05:21  jmlarsen
 * Decreased  warning verbosity
 *
 * Revision 1.18  2007/04/26 13:21:04  jmlarsen
 * Made more robust against inaccurate abs_order polynomial
 *
 * Revision 1.17  2007/04/10 07:11:56  jmlarsen
 * Changed interface of polynomial_regression_2d()
 *
 * Revision 1.16  2007/03/05 10:22:24  jmlarsen
 * Fixed bug in computation of max/min physical order number
 *
 * Revision 1.15  2007/01/15 08:58:20  jmlarsen
 * More robust polynomial fitting
 *
 * Revision 1.14  2007/01/10 12:40:12  jmlarsen
 * Removed unused parameter
 *
 * Revision 1.13  2006/12/07 08:29:58  jmlarsen
 * Compute correct Ynew column for FLAMES
 *
 * Revision 1.12  2006/11/24 16:24:32  jmlarsen
 * Added check of abs order polynomial
 *
 * Revision 1.11  2006/11/15 15:02:15  jmlarsen
 * Implemented const safe workarounds for CPL functions
 *
 * Revision 1.9  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.8  2006/11/06 15:19:42  jmlarsen
 * Removed unused include directives
 *
 * Revision 1.7  2006/08/18 07:07:43  jmlarsen
 * Switched order of cpl_calloc arguments
 *
 * Revision 1.6  2006/07/14 12:43:47  jmlarsen
 * Documentation update
 *
 * Revision 1.5  2006/07/03 13:29:45  jmlarsen
 * Reduced max line length
 *
 * Revision 1.4  2006/03/03 13:54:11  jmlarsen
 * Changed syntax of check macro
 *
 * Revision 1.3  2006/02/15 13:19:15  jmlarsen
 * Reduced source code max. line length
 *
 * Revision 1.2  2006/02/08 09:25:05  jmlarsen
 * Fixed bug caused by == comparison of doubles
 *
 * Revision 1.1  2006/02/03 07:46:30  jmlarsen
 * Moved recipe implementations to ./uves directory
 *
 * Revision 1.27  2005/12/20 08:11:44  jmlarsen
 * Added CVS  entry
 *
 */
/*----------------------------------------------------------------------------*/
/*
 * @addtogroup uves_wavecal
 */
/*----------------------------------------------------------------------------*/
/**@{*/

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

#include <uves_wavecal_firstsolution.h>

#include <uves_utils.h>
#include <uves_utils_wrappers.h>
#include <uves_dump.h>
#include <uves_error.h>
#include <uves_msg.h>

#include <cpl.h>

#include <math.h>

static int *
write_physical_order(cpl_table *linetable,
             const polynomial *absolute_order, 
                     const cpl_table *ordertable,
             const polynomial *order_locations,
             int *first_abs_order, int *last_abs_order);

static double
calculate_shift(const cpl_table *linetable, const cpl_table *previous, 
        const char *column, const char *reference_column, 
        double range, double step, double tolerance);

static double
cross_correlation(double shift, 
          const cpl_table *t1, const cpl_table *t2,
          const char *column, const char* reference_column, 
          int minref, int maxref, double tolerance);

static polynomial *apply_shift(const cpl_table *previous, 
                   const double shift, const int degree, double *mse);

/*----------------------------------------------------------------------------*/
/**
  @brief    Obtain initial dispersion relation
  @param    linetable         The line table containing the already detected 
                              emission lines
  @param    guess             A previously computed line table
  @param    absolute_order    (input/output) A polynomial mapping from pixel
                              coordinates to absolute order number, 
                              @em m = @em f(@em x, @em y). On input the
                              polynomial must contain the initial (guess) mapping,
                              and it will (on output) be updated once the initial
                              dispersion relation is obtained.
  @param    order_locations   A polynomial defining the relative order locations 
                              y = f(x, m)
  @param    flames            FLAMES reduction?
  @param    offset            y-offset (in pixels) of current window wrt to the position
                              defined by @em order_locations
  @param    relative_order    (output) Pointer to an array containing the map from 
                              absolute order number to relative order number:
                              m_relative = f[m_absolute]
  @param    DEGREE            Degree of both independent variables of dispersion 
                              polynomial
  @param    CORREL_RANGE      The maximum shift (in pixels) between current and 
                              guess solution
  @param    CORREL_STEP       The step size (in pixels) when searching for 
                              optimum shift
  @param    CORREL_TOLERANCE  The tolerance (in pixels) when matching shifted lines
  @param    MAXERROR          If the RMS of the initial fit is larger than 
                              this @em MAXERROR, a warning is displayed.
  @param    first_abs_order   The absolute order number of the minimum relative order
  @param    last_abs_order    The absolute order number of the maximum relative order
  @return   The obtained initial dispersion solution @em lambda * @em m = @em f(@em x, @em m), 
            or NULL on error

  This function uses a "guess" line table that contains a map 
  from pixel-order space to wavelengths to obtain an initial dispersion relation,
  @em lambda * @em m = @em f(@em x, @em m), for the provided @em linetable.
  This is achieved by:

  - Calculating the absolute order number, @em m, using the parameter @em absolute_order.
    See @c write_physical_order() .
  - Estimating the x-shift (in pixels) between the "guess" and current line table.
    See @c calculate_shift() .
  - Applying this shift to the "guess" line table, thereby obtaining an initial
    dispersion relation. See @c apply_shift() .

 */
/*----------------------------------------------------------------------------*/
polynomial *
uves_wavecal_firstsolution(cpl_table *linetable,
               const cpl_table *guess, 
               polynomial **absolute_order, 
                           const cpl_table *ordertable,
                           const polynomial *order_locations,
               bool flames,
               double offset,
               int **relative_order, 
               int DEGREE, double CORREL_RANGE, double CORREL_STEP,
               double CORREL_TOLERANCE, double MAXERROR, 
               int *first_abs_order, int *last_abs_order)
{
    polynomial *initial_dispersion = NULL;
    polynomial *new_absorder = NULL;
    const char *er_msg = NULL;
    double shift;
    double mse;

    /* Get physical order numbering */
    check( *relative_order =   write_physical_order(linetable, *absolute_order, 
                                                    ordertable,
                            order_locations,
                            first_abs_order,
                            last_abs_order),
       "Could not calculate absolute order numbers");

    /* Update the 'absolute_order' map */
    {
    int row;

    /* Create column for Y-location (in pixels) of order */
    cpl_table_new_column(linetable, "Ynew", CPL_TYPE_DOUBLE);
    for (row = 0; row < cpl_table_get_nrow(linetable); row++)
        {
        /* For historical reasons, the column 'Y' contains the
           (relative) order number while 'Ynew' contains 
           the y-coordinate (in pixels) of the emission line. */
        int order = cpl_table_get_int   (linetable, "Y", row, NULL);
        double x  = cpl_table_get_double(linetable, "X", row, NULL);
        
        cpl_table_set_double(
            linetable, "Ynew", row, 
            uves_polynomial_evaluate_2d(order_locations, x, order));
        }

    assure_nomsg( cpl_error_get_code() == CPL_ERROR_NONE,
              cpl_error_get_code() );

    new_absorder =
        uves_polynomial_regression_2d(linetable, "X", "Ynew", "Order",
                      NULL,              /* uncertainty of order number */
                      DEGREE, DEGREE,
                      NULL, NULL, NULL,  /* New columns */
                      NULL, NULL,        /* mse, chi^2 */
                      NULL,              /* variance pol. */
                      -1, -1);           /* kappa */
    cpl_table_set_column_unit(linetable,"X","pix");
    cpl_table_set_column_unit(linetable,"Ynew","pix");
    cpl_table_set_column_unit(linetable,"Order"," ");

    if (cpl_error_get_code() != CPL_ERROR_NONE) /* Singular matrix, or too few points */
        {
        er_msg = uves_sprintf("%s", cpl_error_get_message());
        
        uves_error_reset();
        uves_msg_warning("Could not make global fit of absolute order number (%s). "
                 "Polynomial is not updated",
                 er_msg);
        }
    else
        {
        uves_polynomial_delete(absolute_order);
        *absolute_order = uves_polynomial_duplicate(new_absorder);
        }

    /* Calculate absolute_order wrt center of orders, but add offset to Ynew column */
    if (flames)
        {
        cpl_table_add_scalar(linetable, "Ynew", + offset);
        }
    }

    /* Sort linetable by 'Order' (ascending), then 'X' (ascending) */
    uves_sort_table_2(linetable, "Order", "X", false, false);

    /* Cross correlation of guess (linetable) and linetable */
    /* Step size should not be less than 2*tolerance */
    check( shift = calculate_shift(guess, linetable, "X", "Order", 
                   CORREL_RANGE, CORREL_STEP, CORREL_TOLERANCE),
       "Could not calculate shift of position w.r.t. guess solution");

    /* Apply shift to guess solution
     * Note that it doesn't help to simply call uves_polynomial_shift()
     * on the guess solution
     * because the requested 'DEGREE' might be different from
     * the degree used in the guess solution
     */
    
    check( initial_dispersion = apply_shift(guess, shift, DEGREE, &mse),
       "Could not calculate initial dispersion relation");
    /* This fit may fail if the input guess table has too few or badly
       distributed points, but there is not much to do about that */

    /* Check quality of initial solution */
    if(mse > MAXERROR*MAXERROR) 
    {
        uves_msg_warning("RMS of initial fit (%f pixels) is greater "
                 "than tolerance (%f pixels)", sqrt(mse), MAXERROR);
    }
    
  cleanup:
    uves_free_string_const(&er_msg);
    uves_polynomial_delete(&new_absorder);
    if (cpl_error_get_code() != CPL_ERROR_NONE)
    {
        uves_polynomial_delete(&initial_dispersion);
    }
    
    return initial_dispersion;
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Apply x-shift to line table and obtain dispersion relation
  @param    guess             The line table to shift
  @param    shift             The amount (in pixels) to shift the line positions
                              before fitting
  @param    degree            Degree of both independent variables of dispersion polynomial
  @param    mse               (output) Mean squared error of the fit
  @return   The obtained dispersion solution 
            @em lambda * @em m = @em f(@em x, @em m), or NULL on error

  This function shifts all line positions by the amount @em shift, then makes the
  fit of wavelengths.

 */
/*----------------------------------------------------------------------------*/
static polynomial *
apply_shift(const cpl_table *guess, double shift, int degree, double *mse)
{
    polynomial *result = NULL;
    cpl_table *t = NULL;
    
    /* Copy guess table */
    check( t = cpl_table_duplicate(guess),
       "Error duplicating table");
    
    /* Create auxillary column  Ident*Order  */
    check(( cpl_table_duplicate_column(t, "ident_order", t, "Ident"),
        cpl_table_multiply_columns(t, "ident_order", "Order")),
      /* ident_order = Ident * Order */
      "Error creating auxillary column");
    
    /* Shift x values */
    check( cpl_table_add_scalar(t, "X", shift), "Error shifting column 'X'");

    /* Fit lambda*m = f(x, m) */
    /* Don't use uncertainties because they might not exist in guess solution */
    result = uves_polynomial_regression_2d(t, "X", "Order", "ident_order", NULL,
                       degree, degree,
                       NULL, NULL, NULL,
                       mse, NULL,
                       NULL, -1, -1);

    /* If failed, set error to SINGULAR_MATRIX */
    if (cpl_error_get_code() != CPL_ERROR_NONE) /* Singular matrix or too few points */
    {
        uves_error_reset();

        assure( false, CPL_ERROR_SINGULAR_MATRIX,
            "Polynomial fitting failed");
    }

  cleanup:
    uves_free_table(&t);
    return result;
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Calculate shift compared to guess solution
  @param    linetable         Line table containing the detected emission lines
  @param    guess             The "guess" line table
  @param    column            Calculate shift of this numerical column
  @param    reference_column  Values in the numerical @em column are compared
                              when the values in the reference column (which 
                  must by of integer type) are identical.
  @param    CORREL_TOLERANCE  

  @param    range             The maximum shift
  @param    step              The step size
  @param    tolerance         The tolerance (in pixels) when matching lines
  @return   The shift

  The function calculates the shift of the specified @em column (which must
  be double type) between the two tables. Shifts in the interval [-@em range; +@em range]
  are considered using the resolutiong @em step. Two entries match by definition iff
  |@em column_1-@em column_2 - @em shift| < @em tolerance and the values in the
  reference columns are identical.

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

static double
calculate_shift(const cpl_table *linetable, const cpl_table *guess, const char *column,
        const char *reference_column, double range, double step, double tolerance)
{
    cpl_type t;
    int minorder, maxorder;
    int N, i;
    double shift, max_corr, median_corr, maxpos = 0;
    cpl_table *temp = NULL;

    assure( cpl_table_has_column(linetable, column), 
        CPL_ERROR_ILLEGAL_INPUT, "Table has no '%s' column", column);
    assure( cpl_table_has_column(guess , column), 
        CPL_ERROR_ILLEGAL_INPUT, "Table has no '%s' column", column);
    assure( cpl_table_has_column(linetable, reference_column),
        CPL_ERROR_ILLEGAL_INPUT, "Table has no '%s' column", reference_column);
    assure( cpl_table_has_column(guess , reference_column), 
        CPL_ERROR_ILLEGAL_INPUT, "Table has no '%s' column", reference_column);
    assure( range > 0, CPL_ERROR_ILLEGAL_INPUT, "Range = %f", range);

    t = cpl_table_get_column_type(linetable, column);
    assure( t == CPL_TYPE_DOUBLE, CPL_ERROR_TYPE_MISMATCH,
        "Column '%s' has type '%s'. Double expected", column, uves_tostring_cpl_type(t));

    t = cpl_table_get_column_type(guess, column);
    assure( t == CPL_TYPE_DOUBLE, CPL_ERROR_TYPE_MISMATCH,
        "Column '%s' has type '%s'. Double expected", column, uves_tostring_cpl_type(t));

    t = cpl_table_get_column_type(linetable, reference_column);
    assure( t == CPL_TYPE_INT, CPL_ERROR_TYPE_MISMATCH,
        "Ref. column '%s' has type '%s'. Integer expected", 
        reference_column, uves_tostring_cpl_type(t));
    
    t = cpl_table_get_column_type(guess, reference_column);
    assure( t == CPL_TYPE_INT, CPL_ERROR_TYPE_MISMATCH,
        "Ref. column '%s' has type '%s'. Integer expected",
        reference_column, uves_tostring_cpl_type(t));

    /* Identify common orders    */
    check(( minorder = 
        uves_max_int(cpl_table_get_column_min(guess, reference_column), 
             cpl_table_get_column_min(linetable, reference_column)),
        maxorder = 
        uves_min_int(cpl_table_get_column_max(guess, reference_column), 
             cpl_table_get_column_max(linetable, reference_column))),
      "Error reading column '%s'", reference_column);
    
    assure(maxorder >= minorder, CPL_ERROR_ILLEGAL_INPUT, "No common orders found");
    
    uves_msg("Min/max common absolute orders = %d - %d", minorder, maxorder);
    
    /* Find maximum of cross correlation function 
       for shifts in [-range ; range]
    */

    /* Count number of candidates,
       so we can create a table of the correct size
       which is used to get the median of
       all cross-correlation values */
    N = 0;
    for (shift = -range; shift <= range; shift += step) 
    {
        N += 1;
    }

    temp = cpl_table_new(N);
    cpl_table_new_column(temp, "Corr", CPL_TYPE_DOUBLE);

    max_corr = -1;
    maxpos = 0;
    for (shift = -range, i = 0;
     i < N;
     shift += step , i++) 
    {
        double corr;
        check( corr = cross_correlation(shift, linetable, guess, column, 
                        reference_column, minorder, maxorder, tolerance),
           "Error calculating spectrum cross correlation for shift = %f pixel(s)", 
           shift);
        
        /* Update table */
        check( cpl_table_set_double(temp, "Corr", i, corr),
           "Error updating table");
        
        uves_msg_debug("Correlation(shift=%f) = %f", shift, corr);
        
        if (corr > max_corr) 
        {
            max_corr = corr;
            maxpos = shift;
        }
    }

    /* To estimate significance,
       compare the detected max cross-correlation 
       value to "no correlation" estimated as the
       median of all cross-corr. values */

    median_corr = cpl_table_get_column_median(temp, "Corr");
    
    /* Correlation value is integer ; don't divide by zero */
    if (median_corr < 0.5)
    {
        median_corr = 1;
    }

    uves_msg("Estimated shift compared to guess solution is %f pixels (%.2f sigma detection)",
         maxpos, max_corr / median_corr);

    /* The correlation peak is usually 
       ~30 or more times the background,
       so warn if peak value is less than, say,
       10 times background. */
    if (max_corr / median_corr < 10)
    {
        uves_msg_warning("Cross-correlation with guess solution is "
                 "only %f times no correlation (usually >30). "
                 "Make sure that the guess solution is within ~10 pixels "
                 "of the real dispersion relation; otherwise the following "
                 "wavelength calibration is likely to fail or converge "
                 "to a wrong solution",
                 max_corr / median_corr);
    }
    
  cleanup:
    uves_free_table(&temp);
    return maxpos;
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Calculate cross correlation between two table columns
  @param    shift             The cross correlation function is evaluated here
  @param    t1                First table to compare
  @param    t2                Second table to compare
  @param    column            Calculate cross correlation of this numerical column
  @param    reference_column  Values in the numerical @em column are compared
                              only when the values in the reference column (which 
                  must by of integer type) are identical.
  @param    minref            Minimum reference value
  @param    maxref            Maximum reference value
  @param    tolerance         The tolerance when matching values in the 
                              specified @em column
  @return   The number of matching values

  This function returns the number of matching @em column values 
  in @em t1 and @em t2 when the values in @em t1 are shifted by the
  amount @em shift. Two values match iff their difference
  is less than the specified @em tolerance.  

 */
/*----------------------------------------------------------------------------*/
static double
cross_correlation(double shift,
          const cpl_table *t1, const cpl_table *t2,
          const char *column, const char* reference_column, 
          int minref, int maxref, double tolerance)
{
    double result = 0;  /* The result */
    int i1 = 0;         /* Pointers to table rows */
    int i2 = 0;

    /* For efficiency reasons, retrieve the pointers to the columns */
    const double *col1 = cpl_table_get_data_double_const(t1, column);
    const double *col2 = cpl_table_get_data_double_const(t2, column);
    const int *ref1 = cpl_table_get_data_int_const(t1, reference_column);
    const int *ref2 = cpl_table_get_data_int_const(t2, reference_column);

    int N1 = cpl_table_get_nrow(t1);
    int N2 = cpl_table_get_nrow(t2);

    assure( cpl_error_get_code() == CPL_ERROR_NONE, cpl_error_get_code(),
        "Error reading input table");
    
    /* Search for matching rows */
    while (i1 < N1 && ref1[i1] <= maxref && 
       i2 < N2 && ref2[i2] <= maxref) {
    if      (i1 < minref || ref1[i1] < ref2[i2])
        i1++;
    else if (i2 < minref || ref1[i1] > ref2[i2])
        i2++;
    else {
        /* Reference values match */
        double difference = col2[i2] - (col1[i1] + shift);
        
        if      (difference > tolerance)
        {
            i1++;
        }
        else if (difference < -tolerance)
        {
            i2++;
        }
        else {
        /* Matching rows found: |col2-col1-shift| <= tolerance.
           Update result and continue search */
        result += 1.0;
        i2++;
        }
    }
    }


  cleanup:
    return result;
}


/*----------------------------------------------------------------------------*/
/**
  @brief    Get absolute order numbering
  @param    linetable         Line table containing relative order numbers
  @param    absolute_order    Polynomial mapping from pixel coordinates to 
                              absolute order number, @em m = @em f(@em x, @em y)
  @param    order_locations   A polynomial defining the relative order
                              locations y = f(x, m)
  @param    first_abs_order   The absolute order number of the minimum relative order
  @param    last_abs_order    The absolute order number of the maximum relative order

  @return   The array mapping from absolute order number to relative order number

  This function creates a column named 'Order' in the line table which contains
  the absolute order number and is obtained by evaluating the polynomial
  @em absolute_order (x,y) at points along each order. @em x is read from the
  column @em "X" in the line table. The y-position is calculated from @em x and
  the relative order number which is also stored in the line table.   
 */
/*----------------------------------------------------------------------------*/
static int *
write_physical_order(cpl_table *linetable, const polynomial *absolute_order,
                     const cpl_table *ordertable,
             const polynomial *order_locations,
             int *first_abs_order, int *last_abs_order)
{
    int *relative_order = NULL; /* Result */
    int *physical_order = NULL;
    int minorder, maxorder;
    int maxphysical;
    cpl_table *temp = NULL;
    const polynomial *map = NULL;

    double *sum = NULL;   /* Auxillary variables used to calculate the average */
    int      *N = NULL;
    
    int i;

    check( cpl_table_new_column(linetable, "Order", CPL_TYPE_INT),
       "Error creating column");

    check( cpl_table_new_column(linetable, "AbsOrder", CPL_TYPE_DOUBLE),
       "Error creating column");
    
    check(( minorder = cpl_table_get_column_min(ordertable, "Order"),
        maxorder = cpl_table_get_column_max(ordertable, "Order")),
      "Could not read min. and max. order numbers");

    assure( minorder > 0, CPL_ERROR_ILLEGAL_INPUT,
        "Non-positive order number (%d) in linetable", minorder);
    
    physical_order = cpl_calloc(maxorder + 1, sizeof(int));
    assure_mem( physical_order );
    
    /* First calculate the estimation of the
       absolute order number at each line position */
    for (i = 0; i < cpl_table_get_nrow(linetable); i++) {
    double x, y;
    double absorder;
    int order;
    
    order = cpl_table_get_int   (linetable, "Y", i, NULL); 
        /* The column 'Y' contains the (relative) order number */

    x     = cpl_table_get_double(linetable, "X", i, NULL);

    y = uves_polynomial_evaluate_2d(order_locations, x, order);

        absorder = uves_polynomial_evaluate_2d(absolute_order, x, y);

        uves_msg_debug("Order #%d: Absolute order = %f at x = %f",
               order, absorder, x);

        cpl_table_set_double(linetable, "AbsOrder", i, absorder);
    }
 
    {
        int degree = 1;
        int coeff1, coeff2;  /* absorder = coeff1 + coeff2 * relative_order */
        int order;
        int relorder_median;
        int absorder_median;

        check_nomsg( map = 
                     uves_polynomial_regression_1d(linetable,
                                                   "Y", "AbsOrder", NULL,
                                                   degree, 
                                                   NULL, NULL, NULL, -1));
        
        relorder_median = uves_round_double(cpl_table_get_column_median(linetable, "Y"));
        absorder_median = uves_round_double(uves_polynomial_evaluate_1d(map, relorder_median));
            
        if (uves_polynomial_derivative_1d(map, relorder_median) > 0) {
            coeff2 = 1;
        }
        else {
            coeff2 = -1;
        }
        
        coeff1 = absorder_median - coeff2 * relorder_median;

    uves_msg_debug("Assuming relation: abs.order = %d + (%d) * rel.order",
                       coeff1, coeff2);
        
        maxphysical = -1;
        for (order = minorder; order <= maxorder; order++) {
            physical_order[order] = coeff1 + coeff2 * order;
            
            assure(physical_order[order] > 0, CPL_ERROR_ILLEGAL_OUTPUT,
                   "Estimated physical order number is non-positive (%d)", 
                   physical_order[order]);
            
            if (physical_order[order] > maxphysical) 
                {
                    maxphysical = physical_order[order];
                }

            uves_msg_debug("Mapping relative order #%d to absolute order #%d", 
                           order, physical_order[order]);
        }
        
        /* Get first and last physical orders */
        *first_abs_order = physical_order[minorder];
        *last_abs_order  = physical_order[maxorder];
        
        passure( *first_abs_order - *last_abs_order == coeff2*(minorder - maxorder),
                 "%d %d %d %d %d",
                 *first_abs_order, *last_abs_order, coeff2, minorder, maxorder);
        
    }

    /* Then write this rounded mean value to every row of the table */
    for (i = 0; i < cpl_table_get_nrow(linetable); i++) {
    int order;
    order = cpl_table_get_int (linetable, "Y", i, NULL);
    cpl_table_set_int(linetable, "Order", i, physical_order[order]);
    }

    /* Calculate the inverse of 'physical_order' */
    relative_order = cpl_calloc(maxphysical + 1, sizeof(int));
    for (i = 0; i <= maxorder; i++)
    {
        relative_order[physical_order[i]] = i;
    }
    
  cleanup:
    uves_free_table(&temp);
    uves_polynomial_delete_const(&map);
    cpl_free(sum);
    cpl_free(physical_order);
    cpl_free(N);

    return relative_order;
}
/**@}*/