File: uves_extract_profile.c

package info (click to toggle)
cpl-plugin-uves 6.1.3+dfsg-2
  • links: PTS, VCS
  • area: main
  • in suites:
  • size: 23,128 kB
  • sloc: ansic: 171,056; sh: 4,359; python: 3,002; makefile: 1,322
file content (834 lines) | stat: -rw-r--r-- 30,443 bytes parent folder | download | duplicates (4)
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
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
/*                                                                              *
 *   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:03 $
 * $Revision: 1.6 $
 * $Name: not supported by cvs2svn $
 * $Log: not supported by cvs2svn $
 * Revision 1.4  2010/02/13 12:22:31  amodigli
 * removed inlines (let's do work to compiler)
 *
 * Revision 1.3  2007/08/30 07:56:54  amodigli
 * fixed some doxygen warnings
 *
 * Revision 1.2  2007/06/06 08:17:33  amodigli
 * replace tab with 4 spaces
 *
 * Revision 1.1  2007/05/02 13:43:46  jmlarsen
 * Added source
 *
 * Revision 1.141  2007/04/26 06:55:35  amodigli
 * fixed mem leak adding uves_free_image(&spectrum_order)
 *
 * Revision 1.140  2007/04/24 12:50:29  jmlarsen
 * Replaced cpl_propertylist -> uves_propertylist which is much faster
 *
 * Revision 1.139  2007/04/24 09:40:37  jmlarsen
 * Removed deprecated irplib_string_concatenate_all
 *
 * Revision 1.138  2007/04/20 14:44:20  jmlarsen
 * Implemented QC parameter to measure small scale ripples
 *
 * Revision 1.137  2007/04/12 12:00:35  jmlarsen
 * Added testing code
 *
 * Revision 1.136  2007/04/10 11:34:14  jmlarsen
 * Removed debug message
 *
 * Revision 1.135  2007/04/10 08:05:49  jmlarsen
 * Disabled optimization (reduced kappa-sigma iterations, caught by unit test)
 *
 * Revision 1.134  2007/04/10 07:23:20  jmlarsen
 * Added commented out code to spline interpolate virtually resampled profile
 *
 * Revision 1.133  2007/03/28 11:38:38  jmlarsen
 * Removed dead code
 *
 * Revision 1.132  2007/03/19 15:12:14  jmlarsen
 * Optimization: use doubles rather than zero deg. poly.
 *
 * Revision 1.131  2007/03/19 13:50:18  jmlarsen
 * Fixed serious bug happening when object is at +-15 pixels
 *
 * Revision 1.130  2007/03/15 12:33:37  jmlarsen
 * Minor message change
 *
 * Revision 1.129  2007/03/13 15:33:30  jmlarsen
 * Use autodegree polynomials for virtual profile, not zero degree
 *
 * Revision 1.128  2007/03/05 10:16:37  jmlarsen
 * Support slope parameter in 1d fitting
 *
 * Revision 1.127  2007/02/26 13:29:40  jmlarsen
 * Don't use Gauss-Legendre 3 point interpolation, for efficiency
 *
 * Revision 1.126  2007/02/26 11:55:47  jmlarsen
 * Renamed and generalized function uves_raise_to_median() -> uves_raise_to_median_frac()
 *
 * Revision 1.125  2007/02/22 15:33:56  jmlarsen
 * Optimization: use double's rather than constant 2d polynomials
 *
 * Revision 1.124  2007/02/09 13:37:06  jmlarsen
 * Added bug in 2d extraction mode
 *
 * Revision 1.123  2007/02/09 08:14:16  jmlarsen
 * Do not use CPL_PIXEL_MAXVAL which works only for integer images
 *
 * Revision 1.122  2007/02/08 07:33:56  jmlarsen
 * Added doc
 *
 * Revision 1.121  2007/01/31 13:10:33  jmlarsen
 * Changed message
 *
 * Revision 1.120  2007/01/29 12:09:42  jmlarsen
 * Compute QC parameters (pos, fwhm, s/n) also for simple extraction
 *
 * Revision 1.119  2007/01/26 13:49:43  jmlarsen
 * Fixed sky subtraction residuals for optimal sky subtraction
 *
 * Revision 1.118  2007/01/15 08:46:01  jmlarsen
 * Made more robust against extended objects
 *
 * Revision 1.117  2007/01/05 07:22:07  jmlarsen
 * Eliminated compiler warnings
 *
 * Revision 1.116  2007/01/04 13:55:21  jmlarsen
 * Implemented order-by-order object tracing (disabled)
 *
 * Revision 1.115  2006/12/08 07:41:43  jmlarsen
 * Minor doc. change
 *
 * Revision 1.114  2006/11/16 09:48:30  jmlarsen
 * Renamed data type position -> uves_iterate_position, for namespace reasons
 *
 * Revision 1.113  2006/11/15 15:02:14  jmlarsen
 * Implemented const safe workarounds for CPL functions
 *
 * Revision 1.111  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.110  2006/11/08 14:04:34  jmlarsen
 * Implemented flag to select sky subtraction method
 *
 * Revision 1.109  2006/11/06 15:19:41  jmlarsen
 * Removed unused include directives
 *
 * Revision 1.108  2006/10/31 09:14:58  jmlarsen
 * Man page doc fix
 *
 * Revision 1.107  2006/10/02 08:34:40  jmlarsen
 * Do not recompute variance in last iteration
 *
 * Revision 1.106  2006/09/27 15:08:45  jmlarsen
 * Fixed doc. bug
 *
 * Revision 1.105  2006/09/27 13:08:49  jmlarsen
 * Use dynamic memory allocation to store bad pixels
 *
 * Revision 1.104  2006/09/20 12:53:57  jmlarsen
 * Replaced stringcat functions with uves_sprintf()
 *
 * Revision 1.103  2006/09/20 07:25:30  jmlarsen
 * Doc. bug fix
 *
 * Revision 1.102  2006/09/19 14:29:05  jmlarsen
 * Measure object position QC parameter from bottom of slit
 *
 * Revision 1.101  2006/09/19 07:15:35  jmlarsen
 * Added chip to argument list of uves_extract()
 *
 * Revision 1.100  2006/09/11 14:19:28  jmlarsen
 * Updated documentation
 *
 * Revision 1.99  2006/09/11 13:57:46  jmlarsen
 * Remove usage of cpl_image_set after getting bpm pointer
 *
 * Revision 1.98  2006/09/08 14:02:34  jmlarsen
 * Simplified code by using iterators, sky subtraction much optimized
 *
 * Revision 1.97  2006/09/06 15:35:51  jmlarsen
 * Changed indentations
 *
 * Revision 1.96  2006/09/06 14:50:23  jmlarsen
 * Worked on code to globally measure spatial profile
 *
 * Revision 1.95  2006/09/01 13:56:46  jmlarsen
 * Added commented out code (alternative way of measuring spatial profile)
 *
 * Revision 1.94  2006/08/23 15:08:56  jmlarsen
 * Improved plot of spatial profile
 *
 * Revision 1.93  2006/08/23 09:33:03  jmlarsen
 * Renamed local variables shadowing POSIX reserved names
 *
 * Revision 1.92  2006/08/22 15:35:48  jmlarsen
 * Auto-select profile method based on S/N estimate
 *
 * Revision 1.91  2006/08/22 14:20:56  jmlarsen
 * Implemented simultaneous optimal extraction of obj+sky
 *
 * Revision 1.90  2006/08/17 14:40:06  jmlarsen
 * Added missing documentation
 *
 * Revision 1.89  2006/08/17 14:11:25  jmlarsen
 * Use assure_mem macro to check for memory allocation failure
 *
 * Revision 1.88  2006/08/17 13:59:11  jmlarsen
 * Removed CPL2 const bug workaround
 *
 * Revision 1.87  2006/08/17 13:56:52  jmlarsen
 * Reduced max line length
 *
 * Revision 1.86  2006/08/17 09:17:42  jmlarsen
 * Removed CPL2 code
 *
 * Revision 1.85  2006/08/14 12:16:31  jmlarsen
 * Moved defines to top of file
 *
 * Revision 1.84  2006/08/11 14:56:05  amodigli
 * removed Doxygen warnings
 *
 * Revision 1.83  2006/08/11 09:20:06  jmlarsen
 * Implemented workaround for slow cpl_image_set
 *
 * Revision 1.82  2006/08/10 10:49:28  jmlarsen
 * Removed workaround for cpl_image_get_bpm
 *
 * Revision 1.81  2006/08/08 11:02:43  jmlarsen
 * Make temporary copy of image bad pixel map
 *
 * Revision 1.80  2006/08/08 08:19:17  amodigli
 * update to CPL3
 *
 * Revision 1.79  2006/08/07 11:35:35  jmlarsen
 * Disabled parameter environment variable mode
 *
 * Revision 1.78  2006/07/14 12:21:36  jmlarsen
 * Take bad pixels into account in sky subtraction
 *
 * Revision 1.77  2006/07/03 13:01:22  jmlarsen
 * Use analytical-fit sky subtraction method to improve S/N, use a
 * global model of chi square
 *
 * Revision 1.76  2006/06/16 08:23:04  jmlarsen
 * Added comment
 *
 * Revision 1.75  2006/06/05 08:51:55  amodigli
 * cleaned some warnings from static checks
 *
 * Revision 1.74  2006/06/02 06:41:59  jmlarsen
 * Added missing error code
 *
 * Revision 1.73  2006/06/01 14:43:17  jmlarsen
 * Added missing documentation
 *
 * Revision 1.72  2006/05/16 12:13:07  amodigli
 * added QC log
 *
 * Revision 1.71  2006/05/15 08:15:52  jmlarsen
 * Changed default kappa to 10.0
 *
 * Revision 1.70  2006/05/15 07:21:50  jmlarsen
 * Changed default kappa 3.5 -> 5.0
 *
 * Revision 1.69  2006/05/12 15:04:09  jmlarsen
 * Changed gauss/moffat/virtual profile measuring methods to use
 * global polynomials (rather than one polynomial per order)
 *
 * Revision 1.68  2006/04/24 09:21:18  jmlarsen
 * Implemented virtual resampling algorithm
 *
 * Revision 1.67  2006/04/10 12:36:35  jmlarsen
 * Fixed bug that caused extraction to halt if an order is completely 
 * outside an image
 *
 * Revision 1.66  2006/04/07 12:29:21  jmlarsen
 * Bugfix: in opt_evaluate_profile
 *
 * Revision 1.65  2006/04/07 07:10:12  jmlarsen
 * Use Gauss-Legendre rather than Simpson for profile integration
 *
 * Revision 1.64  2006/04/06 11:49:24  jmlarsen
 * Minor msg change
 *
 * Revision 1.63  2006/04/06 08:36:40  jmlarsen
 * Re-factored optimal extraction, added loop to measure 
 * profile until high statistics is achieved
 *
 * Revision 1.62  2006/03/24 14:46:39  jmlarsen
 * Doc. bugfix
 *
 * Revision 1.61  2006/03/24 14:17:37  jmlarsen
 * Mirror input image before/after extraction
 *
 * Revision 1.60  2006/03/03 13:54:11  jmlarsen
 * Changed syntax of check macro
 *
 * Revision 1.59  2006/02/28 09:15:22  jmlarsen
 * Minor update
 *
 * Revision 1.58  2006/02/15 13:19:15  jmlarsen
 * Reduced source code max. line length
 *
 * Revision 1.57  2006/01/25 16:13:20  jmlarsen
 * Changed interface of gauss.fitting routine
 *
 * Revision 1.56  2006/01/12 15:41:14  jmlarsen
 * Moved gauss. fitting to irplib
 *
 * Revision 1.55  2005/12/20 16:10:32  jmlarsen
 * Added some documentation
 *
 * Revision 1.54  2005/12/19 16:17:56  jmlarsen
 * Replaced bool -> int
 *
 */

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

/*----------------------------------------------------------------------------*/
/**
 * @addtogroup uves_extract    Substep: Profile Extraction
 *
 * This module implements simple (i.e. linear, average, weighted)
 * and optimal extraction of echelle spectra.
 */
/*----------------------------------------------------------------------------*/

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

#include <uves_extract_profile.h>

#include <uves_extract_iterate.h>
#include <uves_error.h>

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

uves_extract_profile *
uves_extract_profile_new_constant(double slit_length)
{
    uves_extract_profile *p = NULL;  /* Result */
    
    p = cpl_malloc(sizeof(uves_extract_profile));

    p->constant = true;
    p->slit_length = slit_length;

    /* remaining members not used */

    return p;
}

/*----------------------------------------------------------------------------*/
/**
   @brief    Allocate a spatial profile
   @param    f                       (zero resampling) The profile function.
                                     If NULL is specified, the virtual resampling
                                     algorithm is assumed.
   @param    dfda                    (zero resampling) The derivative of the
                                     profile function
   @param    M                       (zero resampling) Number of parameters for
                                     the profile function
   @param    slit_length             (virtual resampling) extraction slit
   @param    sampling_factor         (virtual resampling) oversampling factor

   @return   The profile which must be deallocated with
             uves_extract_profile_delete

   

*/
/*----------------------------------------------------------------------------*/
uves_extract_profile *
uves_extract_profile_new(int (*f)   (const double x[], const double a[], double *result),
                         int (*dfda)(const double x[], const double a[], double result[]),
                         int M,
                         double slit_length,
                         int sampling_factor)
{
    uves_extract_profile *p = NULL;  /* Result */
    
    p = cpl_malloc(sizeof(uves_extract_profile));

    p->constant = false;
    p->f = f;

    if (f != NULL)
        {
            /* Zero resampling */
            p->dfda      = dfda;
            p->M         = M;
#if ORDER_PER_ORDER
            p->y0        = cpl_calloc(sizeof(polynomial *), 100);  /* 1 poly. per order */
            p->sigma     = cpl_calloc(sizeof(polynomial *), 100);
            p->red_chisq = cpl_calloc(sizeof(polynomial *), 100);
#else
            p->y0        = NULL;
            p->sigma     = NULL;
            p->red_chisq = NULL;
#endif            
            /* Not used */
            p->spatial_bins    = 0;
            p->slit_length     = 0;
            p->sampling_factor = 0;
            p->is_zero_degree  = NULL;
            p->dy_poly         = NULL;
            p->dy_double       = NULL;
            p->current_profile = NULL;
            p->current_ypos    = NULL;
            p->current_interpolated  = NULL;
        }
    else
        {
            /* Virtual resampling */
            p->spatial_bins    = uves_extract_profile_get_nbins(slit_length, sampling_factor);
            p->slit_length     = slit_length;
            p->sampling_factor = sampling_factor;
            p->spatial_bins    = uves_extract_profile_get_nbins(slit_length, sampling_factor);
            p->is_zero_degree  = cpl_calloc(p->spatial_bins, sizeof(bool));
            p->dy_poly         = cpl_calloc(p->spatial_bins, sizeof(polynomial *));
            p->dy_double       = cpl_calloc(p->spatial_bins, sizeof(double));
            p->current_profile = cpl_calloc(p->spatial_bins, sizeof(double));
            p->current_ypos    = cpl_calloc(p->spatial_bins, sizeof(double));
            p->current_interpolated = cpl_calloc(slit_length + 3, sizeof(double));

            /* Not used */
            p->dfda  = NULL;
            p->M     = 0;
            p->y0    = NULL;
            p->sigma = NULL;
        }
    
    return p;
}


/*----------------------------------------------------------------------------*/
/**
   @brief    Delete an optimal profile
   @param    p            double pointer to profile
   @return   void

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

void
uves_extract_profile_delete(uves_extract_profile **p)
{
    if (*p == NULL) return;
    
    if ((*p)->constant) 
        {
            /* nothing to clean */
        }
    else if((*p)->f != NULL) 
        {
#if ORDER_PER_ORDER
/* Then leak some memory */
#else
            uves_polynomial_delete(&((*p)->y0));
            uves_polynomial_delete(&((*p)->sigma));
            uves_polynomial_delete(&((*p)->red_chisq));
#endif
        }
    else
        {
            /* Virtual resampling */
            int i;
            for (i = 0; i < (*p)->spatial_bins; i++)
                {
                    uves_polynomial_delete(& ((*p)->dy_poly[i]) );
                }
            cpl_free((*p)->is_zero_degree);
            cpl_free((*p)->dy_poly);
            cpl_free((*p)->dy_double);
            cpl_free((*p)->current_profile);
            cpl_free((*p)->current_ypos);
            cpl_free((*p)->current_interpolated);
        }
    
    cpl_free(*p);
    *p = NULL;
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Get number of bins (virtual resampling)
  @param    slit_length         extraction slit length
  @param    sampling_factor     oversampling factor
  @return   number of bins across slit
**/
/*----------------------------------------------------------------------------*/
int
uves_extract_profile_get_nbins(double slit_length, int sampling_factor)
{
    return uves_round_double(slit_length + 3) * sampling_factor;
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Convert to y position (virtual resampling)
  @param    pos                   geometry
  @param    bin                   bin number to convert
  @param    sampling_factor       oversampling factor
  @return   y-position

  This function is the inverse of uves_extract_profile_get_bin()
  
**/
/*----------------------------------------------------------------------------*/
double
uves_extract_profile_get_y(uves_iterate_position *pos,
                           double bin,
                           int sampling_factor)
{
    return bin*1.0/sampling_factor + (pos->ycenter - pos->sg.length/2 - 1);
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Convert to bin (virtual resampling)
  @param    pos                 (spatial) position to convert
  @param    sampling_factor     oversampling factor
  @return   bin number
  
  The function converts from (floating point)
  pixel coordinate to (always positive, floating point) bin coordinate.

**/
/*----------------------------------------------------------------------------*/
double
uves_extract_profile_get_bin(const uves_iterate_position *pos,
                             int sampling_factor)
{
    return sampling_factor*(pos->y - (pos->ycenter - pos->sg.length/2 - 1));
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Initialize evalutation of spatial profile
  @param    p          spatial profile
  @param    pos        current bin
  @param    warnings  (messaging only) warnings printed for this order,
                      or NULL
  
  This function must be called before calling @c uves_extract_profile_evaluate().
  (as a 'performance hack').

  The profile is manually normalized to 1, which
  cannot be guaranteed by the zero-resampling algorithm (e.g. sigma == 0.1
  => sum ~= 0)
  or by the virtual resampling algorithms (the polynomials might give *any*
  normalization because of statistical noise).

**/
/*----------------------------------------------------------------------------*/
void
uves_extract_profile_set(const uves_extract_profile *p, 
                         uves_iterate_position *pos,
                         int *warnings)
{
    if (p->constant) {
        ((uves_extract_profile *)p)->current_area = pos->yhigh - pos->ylow + 1;
    }
    else if (p->f != NULL)
        /* Zero */
        {
            double min_sigma = 0.1;

            /* const cast: The profile itself doesn't change */
#if ORDER_PER_ORDER
            check( ((uves_extract_profile *)p)->current_y0 = 
                   pos->ycenter + uves_polynomial_evaluate_1d(p->y0[pos->order-pos->minorder],
                                  pos->x),
                   "Error evaluating polynomial");
#else
            check( ((uves_extract_profile *)p)->current_y0 = 
                   pos->ycenter + uves_polynomial_evaluate_2d(p->y0, pos->x, pos->order),
                   "Error evaluating polynomial");
#endif
            
#if ORDER_PER_ORDER
            check( ((uves_extract_profile *)p)->current_sigma =
                   uves_polynomial_evaluate_1d(p->sigma[pos->order-pos->minorder], pos->x),
                   "Error evaluating polynomial");
#else
            check( ((uves_extract_profile *)p)->current_sigma =
                   uves_polynomial_evaluate_2d(p->sigma, pos->x, pos->order),
                   "Error evaluating polynomial");
#endif
            
            /* Make sure that the inferred 
             * sigma is always 0.1 pixel or more. 
             * Smaller values are unrealistic (undersampled profile) and cause
             * numerical problems (~zero profile area), anyway.
             */
            
            if (p->current_sigma < min_sigma)
                {
                    /* Print only 1 warning per order */
                    if (warnings != NULL && *warnings == 0)
                        {
                            (*warnings)++;
                            uves_msg_warning("Inferred spatial profile width (one sigma) is only "
                                             "%e pixels at (order, x) = (%d, %d). "
                                             "Setting sigma = %.2f pixels",
                                             p->current_sigma, pos->order, pos->x, min_sigma);
                        }
                    
                    ((uves_extract_profile *)p)->current_sigma = min_sigma;
                }

            /* If the profile is well sampled, the 'area' calculated
               below would be 1, but for undersampled profiles (sigma
           much less than 1 pixel) the
               result might differ substantially. Therefore, compute
               the actual sum, and use the correction factor
               later in uves_extract_profile_evaluate().

               The empirical area depends critically upon
               the fractional part of y, so we must do it for every bin.
            */
            {
                double area = 0;
                
                ((uves_extract_profile *)p)->current_area = 1;
                
                area = 0;
                for (pos->y = pos->ylow; pos->y <= pos->yhigh; pos->y++)
                    {
            /* For analytical profiles the results of
               uves_extract_profile_evaluate()
               may range from 1e-300 to ~1

               Such a large range (300 orders of magnitude) is a
               source of problems in the weighted extraction of flat-fields,
               where the resulting flux may end up being only ~1e-300,
               which is "unphysical" and causes infinities after division.

               To always stay on the middle of the road, one might
               decide to approximate small values of the profile to zero,
               for example all values less than 1e-10

               And this would be the place to do it:
            */
                        area += uves_extract_profile_evaluate(p, pos);
                    }

                /* This will not work:    if (area > 0)  
                   If area is very close to zero, we can still get inf.
                   when computing 1/current_area.

                   Therefore set the limit to something  much larger than machine
                   precision, and much less than 1.
                */
                if (area > 1e-10)
                    {
                        ((uves_extract_profile *)p)->current_area = area;
                    }
                else
                    /* Well... the profile must be zero everywhere.
                       To avoid dividing by zero, set the area to something else */
                    {
                        ((uves_extract_profile *)p)->current_area = 1;
                    }
            }
        }
    else
        /* Virtual */
        {
            int i;
            double sum = 0;

            for (i = 0; i < p->spatial_bins; i++)
                {
                    double prof;
                    if (p->is_zero_degree[i])
                        {
                            prof = uves_max_double(0, p->dy_double[i]);
                        }
                    else
                        {
                            /* This is slow */
                            prof = uves_max_double(
                                0, uves_polynomial_evaluate_2d(p->dy_poly[i], 
                                                               pos->x, 
                                                               pos->order));
                        }
                    
                    p->current_ypos[i] = uves_extract_profile_get_y(pos, i, p->sampling_factor);
                    p->current_profile[i] = prof;
                    sum += prof;
                }

            /* Interpolate profile at the positions needed, enforce normalization */
            i = 0;
            sum = 0;
            for (pos->y = pos->ylow; pos->y <= pos->yhigh; pos->y++) 
                {
                    double pint; /* interpolated value */
                    if (false) 
                        /* Nearest bin interpolation (steps, for testing purposes only): */
                        {
                            double bin = uves_extract_profile_get_bin(pos, p->sampling_factor);
                            pint = p->current_profile[uves_round_double(bin)];
                        }
                    else if (true)
                        /* Linear interpolation */
                        /* Interpolate linearly, flux-conserving between two nearest bins 
                         *
                         *   |-----|--|
                         *   bl    b  bu
                         *
                         *  (bl = bin_lower (integer),
                         *   bu = bin_upper (integer), 
                         *   b  = bin       (floating))
                         *
                         *  interpolated = (bu-b)*prof(bl) + (b-bl)*prof(bu)
                         */
                        {
                            double bin = uves_extract_profile_get_bin(pos, 
                                                     p->sampling_factor);
                            
                            int bin_lower = (int) bin;
                            int bin_upper = bin_lower + 1;
                            
                            double prof_lower = p->current_profile[bin_lower];
                            double prof_upper = p->current_profile[bin_upper];
                            
                            double weight = bin_upper - bin;
                            
                            pint = weight*prof_lower + (1-weight)*prof_upper;
                        }
                    else
                        {
                            pint = uves_spline_hermite(
                                pos->y,                         /* Where to interpolate */
                                p->current_ypos, p->current_profile, /* Function */
                                p->spatial_bins,
                                &i);
                        }
                    
                    p->current_interpolated[pos->y - pos->ylow] = pint;
                    sum += pint;
                }

            if ( !(sum > 0) )
                {
                    /* In the exceptional case when sum == 0, 
                       do linear extraction */
                    sum = 1;
                }

            for (pos->y = pos->ylow; pos->y <= pos->yhigh; pos->y++)
                {
                    p->current_interpolated[pos->y - pos->ylow] /= sum;
                }
        }
            
  cleanup:
    return;
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Evalute spatial profile
  @param    pos       current (x, y) position
  @param    profile   spatial profile

  @return the profile at y, guaranteed to be non-negative.
  
  The profile is evaluated at the x previously specified by
  calling @c uves_extract_profile_set().

**/
/*----------------------------------------------------------------------------*/
double
uves_extract_profile_evaluate(const uves_extract_profile *profile,
                              const uves_iterate_position *pos)
{
    double result;

    if (profile->constant) {
        result = 1.0 / profile->current_area;
    }
    else if (profile->f != NULL)
        {
            double a[5];
            
            a[0] = profile->current_y0;
            a[1] = profile->current_sigma;
            a[2] = 1/profile->current_area; /* This is to get a sum of 1
                                               when the profile is summed over
                                               all bins. */
            a[3] = 0.0;                     /* Sky offset             */
            a[4] = 0.0;                     /* Sky offset linear term */
            
            {
                /* Don't use gauss-legendre 3-point. It increases execution time,
                   and makes only extremely small (insignificant) difference on output.
                   
                   Also, the profile was measured using an unbinned analytical profile,
                   so such interpolation probably does not even make sense.
                */
                if (0)
                    {
                        double xp[3]     = {-0.387298334621, 0, 0.387298334621};
                        double weight[3] = {0.2777777777778, 0.444444444444, 0.2777777777778};
                        int i;
                        
                        result = 0;
                        for (i = 0; i < 3; i++)
                            {
                                double val;
                                double y = pos->y;
                                
                                a[0] = profile->current_y0 + xp[i];
                                profile->f(&y, a, &val);
                                result += weight[i] * val;
                            }
                    }
                else
                    {
                        double y = pos->y;
                        
                        a[0] = profile->current_y0;
                        profile->f(&y, a, &result);
                    }
            }
        }
    else
        /* Virtual */
        {
            result = profile->current_interpolated[pos->y - pos->ylow];
        }

    return result;
}

/**@}*/