/*                                                                              *
 *   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;
}

/**@}*/