airyPattern.hpp

Go to the documentation of this file.

/** \file airyPattern.hpp
 * \author Jared R. Males
 * \brief Utilities related to the Airy pattern point spread function.
 * \ingroup gen_math_files
 *
 */
 
//***********************************************************************//
// Copyright 2015, 2016, 2017 Jared R. Males (jaredmales@gmail.com)
//
// This file is part of mxlib.
//
// mxlib 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 3 of the License, or
// (at your option) any later version.
//
// mxlib 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 mxlib.  If not, see <http://www.gnu.org/licenses/>.
//***********************************************************************//
 
#ifndef math_func_airyPattern_hpp
#define math_func_airyPattern_hpp
 
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
 
#include "../constants.hpp"
#include "bessel.hpp"
#include "jinc.hpp"
 
namespace mx
{
 
namespace math
{
 
namespace func
{
 
/// The classical Airy pattern
/** Returns the intensity distribution of the Airy pattern at a given \f$ \lambda/D \f$
 *
 * References: \cite born_and_wolf, \cite mahajan_1986, https://en.wikipedia.org/wiki/Airy_disk.
 *
 * \tparam realT is the floating point type used for arithmetic.
 *
 * \ingroup gen_math_airy_pattern
 */
template <typename realT>
realT airyPattern( realT x /**< [in] is the separation in units of \f$ \lambda/D \f$. */ )
{
    return pow( 2 * jinc( pi<realT>() * x ), 2 );
}
 
/// The centrally obscured Airy pattern
/** Returns the intensity distribution of the centrally obscured Airy pattern at a given \f$ \lambda/D \f$
 *
 * References: \cite born_and_wolf, \cite mahajan_1986, https://en.wikipedia.org/wiki/Airy_disk.
 *
 * \tparam realT is the floating point type used for arithmetic.
 *
 * \ingroup gen_math_airy_pattern
 */
template <typename realT>
realT airyPattern( realT x,  ///< [in] is the separation in units of  \f$ \lambda/D \f$.
                   realT eps ///< [in] is the ratio of the circular central obscuration diameter to the diameter.
)
{
    return ( 1. / pow( 1. - eps * eps, 2 ) ) *
           pow( 2. * jinc( pi<realT>() * x ) - eps * eps * 2. * jinc( eps * pi<realT>() * x ), 2 );
}
 
/// The general centrally obscured Airy pattern, with arbitrary center and platescale.
/** Returns the intensity distribution of the centrally obscured Airy pattern at a given \f$ \lambda/D \f$.  This
 * version allows for an arbitrary center, scaling, and platescale.
 *
 * References: \cite born_and_wolf, \cite mahajan_1986, https://en.wikipedia.org/wiki/Airy_disk.
 *
 * \tparam realT is the floating point type used for arithmetic.
 *
 * \ingroup gen_math_airy_pattern
 */
template <typename realT>
realT airyPattern( realT x,  ///< [in] is the x-coordinate in units of pixels
                   realT y,  ///< [in] is the y-coordinate in units of pixels
                   realT A0, ///< [in] constant value added to the Airy pattern
                   realT A,  ///< [in] peak scale of the Airy pattern.
                   realT x0, ///< [in] is the x-center in units of pixels
                   realT y0, ///< [in] is the y-center in units of pixels
                   realT ps, ///< [in] the platescale in \f$ (\lambda/D)/pixel  \f$
                   realT eps ///< [in] is the ratio of the circular central obscuration diameter to the diameter.
)
{
    realT r = sqrt( pow( x - x0, 2 ) + pow( y - y0, 2 ) ) * ps;
 
    return A0 + A * airyPattern( r, eps );
}
 
/// Fill in an array with the 2D arbitrarily-centered classical Airy pattern
/**
 * At each pixel (x,y) of the array this computes:
 *
 * \f$ f(x,y) = A_0 + A\times Airy(\sqrt(x^2+y^2) \f$
 *
 * \tparam realT is the type to use for arithmetic
 *
 * \ingroup gen_math_airy_pattern
 */
template <typename realT>
void airyPattern2D( realT *arr,     ///< [out] is the allocated array to fill in
                    size_t nx,      ///< [in] is the size of the x dimension of the array (rows)
                    size_t ny,      ///< [in] is the size of the y dimension of the array (columns)
                    const realT A0, ///< [in] is the constant to add to the Gaussian
                    const realT A,  ///< [in] is the scaling factor (peak height = A-A0)
                    const realT x0, ///< [in] is the x-coordinate of the center
                    const realT y0, ///< [in] is the y-coordinate of the center
                    realT ps        ///< [in] the platescale in \f$ (\lambda/D)/pixel  \f$
)
{
    size_t idx;
 
    for( size_t j = 0; j < ny; ++j )
    {
        for( size_t i = 0; i < nx; ++i )
        {
            idx = i + j * nx;
 
            realT rad = sqrt( pow( i - x0, 2 ) + pow( j - y0, 2 ) ) * ps;
 
            arr[idx] = A0 + A * airyPattern( rad );
        }
    }
}
 
/// Seeing Halo Profile
/** A Moffat profile due to Roddier (1981)\cite roddier_1981, see also Racine et al (1999)\cite racine_1999, which
 * can be used to describe the seeing limited point-spread function of a telescope imaging through turbulence, or the
 * halo in partially corrected imaging.
 *
 * \returns the value of the profile at x, in units of fractional flux per unit area.
 *
 * \ingroup gen_math_airy_pattern
 */
template <typename realT>
realT seeingHalo(
    realT x,   ///< [in] the separation in the same units as fwhm.
    realT fwhm ///< [in] the fwhm in arbitrary units.  Note that this defines the area units of the density.
)
{
    return ( 0.488 / ( fwhm * fwhm ) ) * pow( 1. + ( 11. / 6. ) * pow( x / fwhm, 2 ), -11. / 6. );
}
 
/// Calculate the fraction of enclosed power at a given radius for the unobscured Airy Pattern
/** See Mahajan (1986) \cite mahajan_1986 and https://en.wikipedia.org/wiki/Airy_disk.
 *
 * \returns the fraction of enclosed power at radius x
 *
 * \ingroup gen_math_airy_pattern
 */
template <typename realT>
realT airyPatternEnclosed( realT x /**< [in] the radius */ )
{
    realT x1 = x * pi<realT>();
 
    realT b0 = bessel_j<realT>( 0, x1 );
    b0 = b0 * b0;
 
    realT b1 = bessel_j<realT>( 1, x1 );
    b1 = b1 * b1;
 
    realT encp = static_cast<realT>( 1 ) - b0 - b1;
 
    return encp;
}
 
template <typename realT>
realT apeInt( realT x, void *params )
{
    realT eps = *static_cast<realT *>( params );
 
    return bessel_j<realT>( 1, x ) * bessel_j<realT>( 1, eps * x ) / x;
}
 
/// Calculate the fraction of enclosed power at a given radius for the centrally obscured Airy Pattern
/** See Mahajan (1986) \cite mahajan_1986 and https://en.wikipedia.org/wiki/Airy_disk.
 * If eps = 0, this calls the unobscured version.
 *
 * \returns the fraction of enclosed power at radius x
 *
 * \ingroup gen_math_airy_pattern
 */
template <typename realT>
realT airyPatternEnclosed( realT x,  ///< [in] the radius
                           realT eps ///< [in] the central obscuration fraction
)
{
    if( eps == 0 )
        return airyPatternEnclosed( x );
 
    gsl_function func;
    func.function = apeInt<realT>;
    func.params = &eps;
 
    realT jint;
    realT abserr;
    size_t neval;
 
    realT x1 = x * pi<realT>();
 
    gsl_set_error_handler_off();
    gsl_integration_qng( &func, 0, x1, 1e-7, 1e-7, &jint, &abserr, &neval );
 
    realT b0 = bessel_j<realT>( 0, x1 );
    b0 = b0 * b0;
    realT b1 = bessel_j<realT>( 1, x1 );
    b1 = b1 * b1;
 
    realT b0e = bessel_j<realT>( 0, eps * x1 );
    b0e = b0e * b0e;
    realT b1e = bessel_j<realT>( 1, eps * x1 );
    b1e = b1e * b1e;
 
    realT eps2 = pow( eps, 2 );
 
    realT encp = static_cast<realT>( 1 ) - b0 - b1 + eps2 * ( static_cast<realT>( 1 ) - b0e - b1e );
 
    encp = encp - 4 * eps * jint;
    encp = encp / ( static_cast<realT>( 1 ) - eps2 );
 
    return encp;
}
 
} // namespace func
} // namespace math
} // namespace mx
 
#endif // math_func_airyPattern_hpp