gslInterpolation.hpp

Go to the documentation of this file.

/** \file gslInterpolation.hpp
 * \brief Wrappers for using the GNU Scientific Library 1-D interpolation functions
 *
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 * \ingroup gen_math_files
 *
 */
 
//***********************************************************************//
// Copyright 2015-2022 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 gslInterpolation_hpp
#define gslInterpolation_hpp
 
#include <gsl/gsl_interp.h>
#include <gsl/gsl_errno.h>
 
#include <iostream>
#include <vector>
 
namespace mx
{
namespace math
{
 
/// Interpolate a 1-D data X vs Y discrete function onto a new X axis
/**
 * \param interpT one of the <a
 * href="https://www.gnu.org/software/gsl/manual/html_node/Interpolation-Types.html#Interpolation-Types">gsl
 * interpolation types</a>. \param [in]  xin the input x-axis \param [in]  yin the input y-values \param [in]  Nin the
 * size of the input x-y data series \param [in]  xout the desired x-axis \param [out] yout the output interpolated
 * y-values, pre-allocated \param [in]  Nout the size of the output x-y axis
 *
 * \returns 0 on success
 * \returns -1 on a gsl error.
 *
 * \todo report errors iaw mxlib standard in gsl_interpolate
 *
 * \ingroup interpolation
 */
template <typename realT>
int gsl_interpolate(
    const gsl_interp_type *interpT, realT *xin, realT *yin, size_t Nin, realT *xout, realT *yout, size_t Nout )
{
    static_assert( std::is_same<double, typename std::remove_cv<realT>::type>::value,
                   "GSL Interpolation only works with double" );
 
    gsl_set_error_handler_off();
 
    gsl_interp *interp = gsl_interp_alloc( interpT, Nin );
    if( interp == nullptr )
    {
        std::cerr << "gsl_interpolate: gsl_interp_alloc failed\n";
        return -1;
    }
 
    gsl_interp_accel *acc = gsl_interp_accel_alloc();
    if( acc == nullptr )
    {
        std::cerr << "gsl_interpolate: gsl_interp_accel_alloc failed\n";
        return -1;
    }
    int gsl_errno = gsl_interp_init( interp, xin, yin, Nin );
 
    if( gsl_errno != 0 )
    {
        std::cerr << "gsl_interpolate: error from gsl_interp_init [" << gsl_strerror( gsl_errno ) << "]\n";
        gsl_interp_free( interp );
        gsl_interp_accel_free( acc );
        return -1;
    }
 
    gsl_errno = 0;
    for( size_t i = 0; i < Nout; ++i )
    {
        // Don't error check, let it set NAN
        gsl_errno += gsl_interp_eval_e( interp, xin, yin, xout[i], acc, &yout[i] );
    }
 
    gsl_interp_free( interp );
    gsl_interp_accel_free( acc );
 
    if( gsl_errno )
    {
        std::cerr << "gsl_interpolate: error(s) reported by gsl_interp_eval_e\n";
        return -1;
    }
 
    return 0;
}
 
/// Interpolate a 1-D data X vs Y discrete function onto a new X axis (vector version)
/**
 * \param interpT of the <a href="https://www.gnu.org/software/gsl/doc/html/interp.html#d-interpolation-types">gsl
 * interpolation types</a>. \param [in] xin the input x-axis \param [in] yin the input y-values \param [in] xout the
 * desired x-axis \param [out] yout the output interpolated y-values, does not need to be allocated
 *
 * \returns 0 on success or -1 if GSL reports an interpolation error.
 *
 * \ingroup interpolation
 */
template <typename realT>
int gsl_interpolate( const gsl_interp_type *interpT,
                     std::vector<realT> &xin,
                     std::vector<realT> &yin,
                     std::vector<realT> &xout,
                     std::vector<realT> &yout )
{
    yout.resize( xout.size() );
    return gsl_interpolate( interpT, xin.data(), yin.data(), xin.size(), xout.data(), yout.data(), xout.size() );
}
 
} // namespace math
} // namespace mx
 
#endif // gslInterpolation_hpp