File: uves_wavecal_firstsolution.c

package info (click to toggle)
cpl-plugin-uves 6.1.3+dfsg-2
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 23,128 kB
  • sloc: ansic: 171,056; sh: 4,359; python: 3,002; makefile: 1,322
file content (749 lines) | stat: -rw-r--r-- 29,373 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
/*                                                                              *
 *   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;
}
/**@}*/