/* $Id: irplib_strehl.c,v 1.43 2009-11-18 21:37:48 llundin Exp $
 *
 * This file is part of the irplib package
 * Copyright (C) 2002,2003 European Southern Observatory
 *
 * This program 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, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
 */

/*
 * $Author: llundin $
 * $Date: 2009-11-18 21:37:48 $
 * $Revision: 1.43 $
 * $Name: not supported by cvs2svn $
 */

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

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

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

#include <cpl.h>

#include "irplib_utils.h"
#include "irplib_strehl.h"

/*----------------------------------------------------------------------------*/
/**
 * @defgroup irplib_strehl   Functions to compute the Strehl
 */
/*----------------------------------------------------------------------------*/

/*-----------------------------------------------------------------------------
                                   Define
 -----------------------------------------------------------------------------*/

#ifndef IRPLIB_STREHL_RAD_CENTRAL
#define IRPLIB_STREHL_RAD_CENTRAL 5
#endif

#define IRPLIB_DISK_BG_MIN_PIX_NB    30
#define IRPLIB_DISK_BG_REJ_LOW       0.1
#define IRPLIB_DISK_BG_REJ_HIGH      0.1

/*-----------------------------------------------------------------------------
                                   Functions prototypes
 -----------------------------------------------------------------------------*/

static cpl_image * irplib_strehl_generate_otf(double, double, double, double,
                                              int, double);
static double PSF_H1(double, double, double);
static double PSF_H2(double, double);
static double PSF_G(double, double);
static double PSF_sinc(double);
static double PSF_TelOTF(double, double);
static cpl_error_code update_bad_pixel_map(cpl_image* im);


/*-----------------------------------------------------------------------------
                                   Functions code
 -----------------------------------------------------------------------------*/
/**@{*/

/**
  @brief    Update the bad pixel map for the image - mark all pixels with NaN value as a bad
  @param    im             Image with pixel-type float or double
  @return   0 iff ok
  */
cpl_error_code update_bad_pixel_map(cpl_image* im)
{
	int szx = cpl_image_get_size_x(im);
	int szy = cpl_image_get_size_y(im);
	int x = 0;
	cpl_mask* bpm = cpl_image_get_bpm(im);

	for (x = 1; x <=szx; x++)
	{
		int y = 0;
		for(y = 1; y <= szy; y++)
		{
			int isnull = 0;
			double value = cpl_image_get(im, x, y, &isnull);
			if (isnan(value))
			{
				cpl_mask_set(bpm, x, y, CPL_BINARY_1);
			}
		}
	}
	return cpl_error_get_code();
}
/**
  @brief    Compute the strehl ratio in an image.
@param    im             Image with pixel-type float
  @param    m1             Diameter of the M1 mirror [m]
  @param    m2             Diameter of the M2 mirror [m]
  @param    lam            Central wavelength [micron]
  @param    dlam           Filter bandwidth [micron]
  @param    pscale         Positive pixel scale
  @param    size           Size of image to be used for internal PSF
  @param    xpos           The x position of the ring center
  @param    ypos           The y position of the ring center
  @param    r1             The Star Radius, r1 > 0
  @param    r2             The Internal Radius, r2 > 0
  @param    r3             The External Radius, r3 > r2
  @param    noise_box_sz   Pass -1 for default values
  @param    noise_nsamples Pass -1 for default values
  @param    strehl         Pointer to the Strehl Error (positive on success)
  @param    strehl_err     Pointer to the Strehl Error (non-negative on success)
  @param    star_bg        Pointer to the Star Background
  @param    star_peak      Pointer to the Star Peak (positive on success)
  @param    star_flux      Pointer to the Star Flux (positive on success)
  @param    psf_peak       Pointer to the PSF Peak (positive on success)
  @param    psf_flux       Pointer to the PSF Flux (1 on success)
  @param    bg_noise       Pointer to the Background Noise
  @return   0 iff ok
  @note     The output is undefined on error. On success the Strehl Ratio may
            exceed 1. Before computation, all NaN values in the image would be marked in the
            bad pixel map. That means, parameter im would changed (bad pixel map would be updated / created).
  */

cpl_error_code irplib_strehl_mark_bad_and_compute(cpl_image *   im,
                                     double              m1,
                                     double              m2,
                                     double              lam,
                                     double              dlam,
                                     double              pscale,
                                     int                 size,
                                     double              xpos,
                                     double              ypos,
                                     double              r1,
                                     double              r2,
                                     double              r3,
                                     int                 noise_box_sz,
                                     int                 noise_nsamples,
                                     double          *   strehl,
                                     double          *   strehl_err,
                                     double          *   star_bg,
                                     double          *   star_peak,
                                     double          *   star_flux,
                                     double          *   psf_peak,
                                     double          *   psf_flux,
                                     double          *   bg_noise)
{
	cpl_ensure_code(!update_bad_pixel_map(im), cpl_error_get_code());
	return irplib_strehl_compute(im, m1, m2, lam, dlam, pscale, size, xpos, ypos,
            r1,
            r2,
            r3,
            noise_box_sz,
            noise_nsamples,
            strehl,
            strehl_err,
            star_bg,
            star_peak,
            star_flux,
            psf_peak,
            psf_flux,
            bg_noise);
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Compute the strehl ratio in an image
  @param    im             Image with pixel-type float
  @param    m1             Diameter of the M1 mirror [m]
  @param    m2             Diameter of the M2 mirror [m]
  @param    lam            Central wavelength [micron]
  @param    dlam           Filter bandwidth [micron]
  @param    pscale         Positive pixel scale [Arcsecond/pixel]
  @param    size           Size of image to be used for internal PSF [pixel]
  @param    xpos           The x position of the ring center [pixel]
  @param    ypos           The y position of the ring center [pixel]
  @param    r1             The Star Radius, r1 > 0 [Arcsecond]
  @param    r2             The Internal Radius, r2 > 0 [Arcsecond]
  @param    r3             The External Radius, r3 > r2 [Arcsecond]
  @param    noise_box_sz   Pass -1 for default values
  @param    noise_nsamples Pass -1 for default values
  @param    strehl         Pointer to the Strehl Error (positive on success)
  @param    strehl_err     Pointer to the Strehl Error (non-negative on success)
  @param    star_bg        Pointer to the Star Background
  @param    star_peak      Pointer to the Star Peak (positive on success)
  @param    star_flux      Pointer to the Star Flux (positive on success)
  @param    psf_peak       Pointer to the PSF Peak (positive on success)
  @param    psf_flux       Pointer to the PSF Flux (1 on success)
  @param    bg_noise       Pointer to the Background Noise
  @return   0 iff ok
  @note     The output is undefined on error. On success the Strehl Ratio may
            exceed 1. Strehl computation could fail if the image contains NaN values
            and those values are not marked in the bad pixel map. In that case the function
            rplib_strehl_mark_bad_and_compute could be used instead.
 */
/*----------------------------------------------------------------------------*/
cpl_error_code irplib_strehl_compute(const cpl_image *   im,
                                     double              m1,
                                     double              m2,
                                     double              lam,
                                     double              dlam,
                                     double              pscale,
                                     int                 size,
                                     double              xpos,
                                     double              ypos,
                                     double              r1,
                                     double              r2,
                                     double              r3,
                                     int                 noise_box_sz,
                                     int                 noise_nsamples,
                                     double          *   strehl,
                                     double          *   strehl_err,
                                     double          *   star_bg,
                                     double          *   star_peak,
                                     double          *   star_flux,
                                     double          *   psf_peak,
                                     double          *   psf_flux,
                                     double          *   bg_noise)
{
    cpl_image  * psf;
    double       star_radius, max_radius;

    /* FIXME: Arbitrary choice of image border */
    const double window_size = (double)(IRPLIB_STREHL_RAD_CENTRAL);

    /* Determined empirically by C. Lidman for Strehl error computation */
    const double strehl_error_coefficient = CPL_MATH_PI * 0.007 / 0.0271;
    double       ring[4];
    /* cpl_flux_get_noise_ring() must succeed with this many tries */
    int          ring_tries = 3;
    cpl_errorstate prestate;

    /* Check compile-time constant */
    cpl_ensure_code(window_size > 0.0,  CPL_ERROR_ILLEGAL_INPUT);

    /* Test inputs */
    cpl_ensure_code(im != NULL,         CPL_ERROR_NULL_INPUT);
    cpl_ensure_code(strehl != NULL,     CPL_ERROR_NULL_INPUT);
    cpl_ensure_code(strehl_err != NULL, CPL_ERROR_NULL_INPUT);
    cpl_ensure_code(star_bg != NULL,    CPL_ERROR_NULL_INPUT);
    cpl_ensure_code(star_peak != NULL,  CPL_ERROR_NULL_INPUT);
    cpl_ensure_code(star_flux != NULL,  CPL_ERROR_NULL_INPUT);
    cpl_ensure_code(psf_peak != NULL,   CPL_ERROR_NULL_INPUT);
    cpl_ensure_code(psf_flux != NULL,   CPL_ERROR_NULL_INPUT);

    cpl_ensure_code(pscale > 0.0,      CPL_ERROR_ILLEGAL_INPUT);

    cpl_ensure_code(r1 > 0.0,      CPL_ERROR_ILLEGAL_INPUT);
    cpl_ensure_code(r2 > 0.0,      CPL_ERROR_ILLEGAL_INPUT);
    cpl_ensure_code(r3 > r2,       CPL_ERROR_ILLEGAL_INPUT);

    /* Computing a Strehl ratio is a story between an ideal PSF */
    /* and a candidate image supposed to approximate this ideal PSF. */

    /* Generate first appropriate PSF to find max peak */
    psf = irplib_strehl_generate_psf(m1, m2, lam, dlam, pscale, size);
    cpl_ensure_code(psf != NULL,      CPL_ERROR_ILLEGAL_OUTPUT);

    /* Compute flux in PSF and find max peak */
    *psf_peak = cpl_image_get_max(psf);
    cpl_image_delete(psf);

    assert( *psf_peak > 0.0); /* The ideal PSF has a positive maximum */
    *psf_flux = 1.0; /* The psf flux, cpl_image_get_flux(psf), is always 1 */

    /* Measure the background in the candidate image */
    *star_bg = irplib_strehl_ring_background(im, xpos, ypos, r2/pscale, r3/pscale,
                                             IRPLIB_BG_METHOD_AVER_REJ);

    /* Compute star_radius in pixels */
    star_radius = r1/pscale;

    /* Measure the flux on the candidate image */
    *star_flux = irplib_strehl_disk_flux(im, xpos, ypos, star_radius, *star_bg);

    cpl_ensure_code(*star_flux > 0.0,      CPL_ERROR_ILLEGAL_OUTPUT);

    /* Find the peak value on the central part of the candidate image */
    max_radius = window_size < star_radius ? window_size : star_radius;
    cpl_ensure_code(!irplib_strehl_disk_max(im, xpos, ypos, max_radius,
                                            star_peak), cpl_error_get_code());
    *star_peak -= *star_bg;

    cpl_ensure_code(*star_peak > 0.0,      CPL_ERROR_ILLEGAL_OUTPUT);

    /* Compute Strehl */
    /* (StarPeak / StarFlux) / (PsfPeak / PsfFlux) */
    *strehl = (*star_peak * *psf_flux ) / ( *star_flux * *psf_peak);

    if (*strehl > 1)
        cpl_msg_warning(cpl_func, "Extreme Strehl-ratio=%g, star_peak=%g, "
                        "star_flux=%g, psf_peak=%g, psf_flux=%g", *strehl,
                        *star_peak, *star_flux, *psf_peak, *psf_flux);

    /* Compute Strehl error */
    /* computation could fail if the image contains pixels with NaN value*/
    ring[0] = xpos;
    ring[1] = ypos;
    ring[2] = r2/pscale;
    ring[3] = r3/pscale;

    /* FIXME: With CPL 5.1 the recoverable error
       will be CPL_ERROR_DATA_NOT_FOUND */
    prestate = cpl_errorstate_get();
    while (cpl_flux_get_noise_ring(im, ring, noise_box_sz, noise_nsamples,
                                   bg_noise, NULL) && --ring_tries > 0);
    if (ring_tries > 0) {
        cpl_errorstate_set(prestate); /* Recover, if an error happened */
    } else {
        return cpl_error_set_where(cpl_func);
    }

    *strehl_err = strehl_error_coefficient * (*bg_noise) * pscale *
        star_radius * star_radius / *star_flux;

    /* This check should not be able to fail, but just to be sure */
    cpl_ensure_code(*strehl_err >= 0.0,       CPL_ERROR_ILLEGAL_OUTPUT);

    return CPL_ERROR_NONE;
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Compute the flux from a disk
  @param    im          Image with pixel-type float
  @param    xpos        the x position of the disk center
  @param    ypos        the y position of the disk center
  @param    rad         the radius
  @param    bg          the star background
  @return   The computed flux or 0.0 on error
  @note     (xpos, ypos) and may be outside the image, if so then a sufficiently
            small rad will cause no pixels to be encircled, in which case the
            returned flux is zero.

  The nan values check has been added for SINFONI. The input image in
  input should not contain nan values, but it does in  some cases...
 */
/*----------------------------------------------------------------------------*/
double irplib_strehl_disk_flux(const cpl_image * im,
                               double            xpos,
                               double            ypos,
                               double            rad,
                               double            bg)
{
    const double    sqr = rad * rad;
    double          sqrest;
    const float *   pim;
    double          flux = 0.0;
    double          yj, xi;
    int             nx, ny;
    int             lx, ly, ux, uy;
    int             i, j;


    /* Check entries */
    cpl_ensure(im != NULL, CPL_ERROR_NULL_INPUT, 0.0);
    cpl_ensure(cpl_image_get_type(im) == CPL_TYPE_FLOAT,
               CPL_ERROR_UNSUPPORTED_MODE, 0.0);
    cpl_ensure(rad > 0.0, CPL_ERROR_ILLEGAL_INPUT, 0.0);

    nx = cpl_image_get_size_x(im);
    ny = cpl_image_get_size_y(im);

    /* Round down */
    lx = (int)(xpos - rad);
    ly = (int)(ypos - rad);
    if (lx < 0) lx = 0;
    if (ly < 0) ly = 0;

    /* Round up */
    ux = (int)(xpos + rad) + 1;
    uy = (int)(ypos + rad) + 1;
    if (ux > (nx-1)) ux = nx-1;
    if (uy > (ny-1)) uy = ny-1;

    pim = cpl_image_get_data_float_const(im);
    for (j=ly ; j<uy ; j++) {
        yj = (double)j - ypos;
        sqrest = sqr - yj * yj;
        for (i=lx; i<ux ; i++) {
            xi = (double)i - xpos;
            if (sqrest >= xi * xi && irplib_isnan(pim[i+j*nx]) == 0) {
                flux += (double)pim[i+j*nx] - bg;
            }
        }
    }
    return flux;
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Compute the background in the image from a specified disk
  @param    im          Image with pixel-type float
  @param    xpos        the x position of the ring center
  @param    ypos        the y position of the ring center
  @param    rad_int     the internal radius
  @param    rad_ext     the external radius
  @param    mode        IRPLIB_BG_METHOD_AVER_REJ or IRPLIB_BG_METHOD_MEDIAN
  @return   The computed background or 0.0 on error.

  The nan values check has been added for SINFONI. The input image in
  input should not contain nan values, but it does in  some cases...
 */
/*----------------------------------------------------------------------------*/
double irplib_strehl_ring_background(const cpl_image *       im,
                                     double                  xpos,
                                     double                  ypos,
                                     double                  rad_int,
                                     double                  rad_ext,
                                     irplib_strehl_bg_method mode)
{
    int             npix;
    const double    sqr_int = rad_int * rad_int;
    const double    sqr_ext = rad_ext * rad_ext;
    double          dist;
    cpl_vector  *   pix_arr;
    const float *   pim;
    double          flux = 0.0;
    double          yj, xi;
    int             lx, ly, ux, uy;
    int             nx, ny;
    int             i, j;

    /* Check entries */
    cpl_ensure(im != NULL, CPL_ERROR_NULL_INPUT, 0.0);
    cpl_ensure(cpl_image_get_type(im) == CPL_TYPE_FLOAT,
               CPL_ERROR_UNSUPPORTED_MODE, 0.0);
    cpl_ensure(rad_int > 0.0, CPL_ERROR_ILLEGAL_INPUT, 0.0);
    cpl_ensure(rad_ext > rad_int, CPL_ERROR_ILLEGAL_INPUT, 0.0);

    cpl_ensure(mode == IRPLIB_BG_METHOD_AVER_REJ ||
               mode == IRPLIB_BG_METHOD_MEDIAN,
               CPL_ERROR_UNSUPPORTED_MODE, 0.0);

    nx = cpl_image_get_size_x(im);
    ny = cpl_image_get_size_y(im);

    /* Round down */
    lx = (int)(xpos - rad_ext);
    ly = (int)(ypos - rad_ext);
    if (lx < 0) lx = 0;
    if (ly < 0) ly = 0;

    /* Round up */
    ux = (int)(xpos + rad_ext) + 1;
    uy = (int)(ypos + rad_ext) + 1;
    if (ux > (nx-1)) ux = nx-1;
    if (uy > (ny-1)) uy = ny-1;

    npix = (ux - lx + 1) * (uy - ly + 1);
    cpl_ensure(npix >= IRPLIB_DISK_BG_MIN_PIX_NB, CPL_ERROR_DATA_NOT_FOUND, 0.0);

    /* Allocate pixel array to hold values in the ring */
    pix_arr = cpl_vector_new(npix);

    /* Count number of pixels in the ring */
    /* Retrieve all pixels which belong to the ring */
    pim = cpl_image_get_data_float_const(im);
    npix = 0;
    for (j=ly ; j<uy ; j++) {
        yj = (double)j - ypos;
        for (i=lx ; i<ux; i++) {
            xi = (double)i - xpos;
            dist = yj * yj + xi * xi;
            if (sqr_int <= dist && dist <= sqr_ext &&
                irplib_isnan(pim[i+j*nx]) == 0) {
                cpl_vector_set(pix_arr, npix, (double)pim[i+j*nx]);
                npix++;
            }
        }
    }

    if (npix < IRPLIB_DISK_BG_MIN_PIX_NB) {
        cpl_vector_delete(pix_arr);
        cpl_ensure(0, CPL_ERROR_DATA_NOT_FOUND, 0.0);
    }

    /* Should not be able to fail now */

    /* Resize pixel array to actual number of values within the ring */
    cpl_vector_set_size(pix_arr, npix);

    if (mode == IRPLIB_BG_METHOD_AVER_REJ) {
        const int low_ind  = (int)((double)npix * IRPLIB_DISK_BG_REJ_LOW);
        const int high_ind = (int)((double)npix
                                   * (1.0 - IRPLIB_DISK_BG_REJ_HIGH));

        /* Sort the array */
        cpl_vector_sort(pix_arr, 1);

        for (i=low_ind ; i<high_ind ; i++) {
            flux += cpl_vector_get(pix_arr, i);
        }
        if (high_ind - low_ind > 1) flux /= (double)(high_ind - low_ind);
    } else /* if (mode == IRPLIB_BG_METHOD_MEDIAN) */ {
        flux = cpl_vector_get_median(pix_arr);
    }

    cpl_vector_delete(pix_arr);

    return flux;
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Compute the ideal PSF for a given telescope+instrument.
  @param    m1          Diameter of the M1 mirror [meter]
  @param    m2          Diameter of the M2 mirror [meter]
  @param    lam         Central wavelength [micron]
  @param    dlam        Filter bandwidth [micron]
  @param    size        Generated image size (image will be square) [pixel]
  @param    pscale      Pixel scale on the sky [Arcsecond/pixel]
  @return   1 newly generated image.

  This function computes the ideal PSF for a given telescope and instrument.
  The PSF is computed by first generated the ideal OTF for the provided
  conditions, and applying a Fourier transform to it to bring it back to real
  space. The returned PSF is normalized to unity flux, to help Strehl ratio
  computations.

  The image halves of the returned PSF are swapped in both directions.

 */
/*----------------------------------------------------------------------------*/
cpl_image * irplib_strehl_generate_psf(
        double   m1,
        double   m2,
        double   lam,
        double   dlam,
        double   pscale,
        int      size)
{
    cpl_image * otf_image = irplib_strehl_generate_otf(m1, m2, lam, dlam,
            size, pscale);

    if (otf_image == NULL) return NULL;

    /* Transform back to real space
       - Normalization is unnecessary, due to the subsequent normalisation.
       - An OTF is point symmetric about its center, i.e. it is even,
         i.e. the real space image is real.
       - Because of this a forward FFT works as well.
       - If the PSF ever needs to have its images halves swapped add
         CPL_FFT_SWAP_HALVES to the FFT call.
     */

    if (cpl_image_fft(otf_image, NULL, CPL_FFT_UNNORMALIZED) ||

        /* Compute absolute values of PSF */
        cpl_image_abs(otf_image) ||

        /* Normalize PSF to get flux=1  */
        cpl_image_normalise(otf_image, CPL_NORM_FLUX)) {

        cpl_image_delete(otf_image);
        return NULL;
    }

    return otf_image;
}

/**@}*/

/*----------------------------------------------------------------------------*/
/**
  @brief    Generate an image of an ideal Optical Transfer Function.
  @param    m1          Diameter of the M1 mirror [m]
  @param    m2          Diameter of the M2 mirror [m], m2 < m1
  @param    lam         Central wavelength [micron]
  @param    dlam        Filter bandwidth [micron]
  @param    size        Generated image size (image will be square).
  @param    pscale      Pixel scale on the sky [Arcsecond/pixel]
  @return   1 newly generated cpl_image

  Based on the paper "Amplitude estimation from speckle interferometry" by
  Christian Perrier in "Diffraction-limited imaging with very large telescopes",
  NATO ASI Series C, Vol. 274, edited by D. Alloin and J.-M. Mariotti, 1989
  (p. 99).
 */
/*----------------------------------------------------------------------------*/
static cpl_image * irplib_strehl_generate_otf(
        double  m1,
        double  m2,
        double  lam,
        double  dlam,
        int     size,
        double  pscale)
{
    cpl_image   *   otf_image;
    double      *   otf_data;
    double          obs_ratio ;  /* m1 / m2    */
    double          f_max ;      /* cut-off frequency        */
    int             pix0 ;      /* Pixel corresponding to the zero frequency */
    double          a, x, y;
    double          f, rsq, fc, invfc, lambda;
    double          sincy;
    double          invsize;
    register int    pos;
    int             i, j, k;


    cpl_ensure(m2   > 0.0,      CPL_ERROR_ILLEGAL_INPUT, NULL);
    cpl_ensure(m1   > m2,       CPL_ERROR_ILLEGAL_INPUT, NULL);
    cpl_ensure(lam  > 0.0,      CPL_ERROR_ILLEGAL_INPUT, NULL);
    cpl_ensure(dlam > 0.0,      CPL_ERROR_ILLEGAL_INPUT, NULL);
    cpl_ensure(size > 0,        CPL_ERROR_ILLEGAL_INPUT, NULL);
    cpl_ensure(pscale > 0.0,    CPL_ERROR_ILLEGAL_INPUT, NULL);

    /* Convert pixel scale from sec to radians, microns in meters    */
    pscale /= (double)206265;
    lam /= (double)1.0e6;
    dlam /= (double)1.0e6;

    /* Obscuration ratio    */
    obs_ratio = m2 / m1;

    /* Pixel corresponding to the zero frequency    */
    pix0 = size/2;
    invsize = (double)1.0 / (double)size;

    /* Cut-off frequency in pixels  */
    f_max = m1 * pscale * (double)size / lam;

    /* Allocate for output image    */
    otf_image = cpl_image_new(size, size, CPL_TYPE_DOUBLE);
    if (otf_image==NULL) return NULL;
    otf_data = cpl_image_get_data_double(otf_image);

    /* Now compute the OTF  */
    /* OPTIMIZED CODE !!! LIMITED READABILITY !!!   */

    for (k=1 ; k<=9 ; k++) {    /* iteration on the wavelength  */
        /* Compute intermediate cut-off frequency   */
        lambda = (double)(lam - dlam*(double)(k-5)/8.0);
        fc = (double)f_max * (double)lam / lambda;
        invfc = 1.0 / fc;

        /* Convolution with the detector pixels */
        pos = 0;
        for (j=0 ; j<size ; j++) {
            y = (double)(j-pix0);
            sincy = PSF_sinc(CPL_MATH_PI * y * invsize);
            for (i=0 ; i<size ; i++) {
                x = (double)(i-pix0);
                rsq = x*x + y*y;
                if (rsq < fc*fc) {
                    if (rsq < 0.01)
                        a = 1.0;
                    else {
                        f = sqrt(rsq) * invfc;
                        a = PSF_TelOTF(f,obs_ratio) *
                            PSF_sinc(CPL_MATH_PI * x * invsize) * sincy;
                    }
                } else {
                    a = 0.0;
                }
                otf_data[pos++] += a / 9.0;
            }
        }
    }
    return otf_image;
}

/*----------------------------------------------------------------------------*
 * H1 function
 *----------------------------------------------------------------------------*/
static double PSF_H1(
        double  f,
        double  u,
        double  v)
{
    const double e = fabs(1.0-v) > 0.0 ? -1.0 : 1.0; /* e = 1.0 iff v = 1.0 */

    return((v*v/CPL_MATH_PI)*acos((f/v)*(1.0+e*(1.0-u*u)/(4.0*f*f))));
}

/*----------------------------------------------------------------------------*
 * H2 function
 *----------------------------------------------------------------------------*/
static double PSF_H2(double  f,
                     double  u)
{
    const double tmp1 = (2.0 * f) / (1.0 + u);
    const double tmp2 = (1.0 - u) / (2.0 * f);

    return -1.0 * (f/CPL_MATH_PI) * (1.0+u)
        * sqrt((1.0-tmp1*tmp1)*(1.0-tmp2*tmp2));
}

/*----------------------------------------------------------------------------*
 * G function
 *----------------------------------------------------------------------------*/
static double PSF_G(double  f,
                    double  u)
{
    if (f <= (1.0-u)/2.0) return(u*u);
    if (f >= (1.0+u)/2.0) return(0.0);
    else return(PSF_H1(f,u,1.0) + PSF_H1(f,u,u) + PSF_H2(f,u));
}

/*----------------------------------------------------------------------------*
 * sinc function
 *----------------------------------------------------------------------------*/
static double PSF_sinc(double x)
{
  return fabs(x) > fabs(sin(x)) ? sin(x)/x : 1.0;
}

/*----------------------------------------------------------------------------*
 * Telescope OTF function
 *----------------------------------------------------------------------------*/
static double PSF_TelOTF(double  f,
                         double  u)
{
    return((PSF_G(f,1.0)+u*u*PSF_G(f/u,1.0)-2.0*PSF_G(f,u))/(1.0-u*u));
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Find the peak flux inside a disk
  @param    self        Image with pixel-type float
  @param    xpos        The x position of the disk center
  @param    ypos        The y position of the disk center
  @param    radius      The positive radius
  @param    ppeak       Pointer to the Star Peak
  @return   0 iff successful
  @note *ppeak is undefined on error

 */
/*----------------------------------------------------------------------------*/
cpl_error_code irplib_strehl_disk_max(const cpl_image * self,
                                             double            xpos,
                                             double            ypos,
                                             double            radius,
                                             double          * ppeak)
{

    const double    sqr = radius * radius;
    double          sqrest;
    const float *   pself;
    float           peak = FLT_MAX;  /* Avoid (false) uninit warning */
    double          yj, xi;
    int             nx, ny;
    int             lx, ly, ux, uy;
    int             i, j;
    cpl_boolean     first = CPL_TRUE;


    /* Check entries */
    cpl_ensure_code(ppeak != NULL, CPL_ERROR_NULL_INPUT);
    cpl_ensure_code(self  != NULL, CPL_ERROR_NULL_INPUT);
    cpl_ensure_code(cpl_image_get_type(self) == CPL_TYPE_FLOAT,
                    CPL_ERROR_UNSUPPORTED_MODE);
    cpl_ensure_code(radius > 0.0, CPL_ERROR_ILLEGAL_INPUT);

    nx = cpl_image_get_size_x(self);
    ny = cpl_image_get_size_y(self);

    /* Round down */
    lx = (int)(xpos - radius);
    ly = (int)(ypos - radius);
    if (lx < 0) lx = 0;
    if (ly < 0) ly = 0;

    /* Round up */
    ux = (int)(xpos + radius) + 1;
    uy = (int)(ypos + radius) + 1;
    if (ux > (nx-1)) ux = nx-1;
    if (uy > (ny-1)) uy = ny-1;

    pself = cpl_image_get_data_float_const(self);
    for (j=ly ; j<uy ; j++) {
        yj = (double)j - ypos;
        sqrest = sqr - yj * yj;
        for (i=lx; i<ux ; i++) {
            xi = (double)i - xpos;
            if (sqrest >= xi * xi && irplib_isnan(pself[i+j*nx]) == 0) {
                if (first || pself[i+j*nx] > peak) {
                    first = CPL_FALSE;
                    peak = pself[i+j*nx];
                }
            }
        }
    }

    cpl_ensure_code(!first, CPL_ERROR_DATA_NOT_FOUND);

    *ppeak = (double)peak;

    return CPL_ERROR_NONE;
}