Line data Source code
1 : /** \file gslInterpolation.hpp
2 : * \brief Wrappers for using the GNU Scientific Library 1-D interpolation functions
3 : *
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : * \ingroup gen_math_files
7 : *
8 : */
9 :
10 : //***********************************************************************//
11 : // Copyright 2015-2022 Jared R. Males (jaredmales@gmail.com)
12 : //
13 : // This file is part of mxlib.
14 : //
15 : // mxlib is free software: you can redistribute it and/or modify
16 : // it under the terms of the GNU General Public License as published by
17 : // the Free Software Foundation, either version 3 of the License, or
18 : // (at your option) any later version.
19 : //
20 : // mxlib is distributed in the hope that it will be useful,
21 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
22 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 : // GNU General Public License for more details.
24 : //
25 : // You should have received a copy of the GNU General Public License
26 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
27 : //***********************************************************************//
28 :
29 : #ifndef gslInterpolation_hpp
30 : #define gslInterpolation_hpp
31 :
32 : #include <gsl/gsl_interp.h>
33 : #include <gsl/gsl_errno.h>
34 :
35 : #include <iostream>
36 : #include <vector>
37 :
38 : namespace mx
39 : {
40 : namespace math
41 : {
42 :
43 : /// Interpolate a 1-D data X vs Y discrete function onto a new X axis
44 : /**
45 : * \param interpT one of the <a
46 : * href="https://www.gnu.org/software/gsl/manual/html_node/Interpolation-Types.html#Interpolation-Types">gsl
47 : * interpolation types</a>. \param [in] xin the input x-axis \param [in] yin the input y-values \param [in] Nin the
48 : * size of the input x-y data series \param [in] xout the desired x-axis \param [out] yout the output interpolated
49 : * y-values, pre-allocated \param [in] Nout the size of the output x-y axis
50 : *
51 : * \returns 0 on success
52 : * \returns -1 on a gsl error.
53 : *
54 : * \todo report errors iaw mxlib standard in gsl_interpolate
55 : *
56 : * \ingroup interpolation
57 : */
58 : template <typename realT>
59 0 : int gsl_interpolate(
60 : const gsl_interp_type *interpT, realT *xin, realT *yin, size_t Nin, realT *xout, realT *yout, size_t Nout )
61 : {
62 : static_assert( std::is_same<double, typename std::remove_cv<realT>::type>::value,
63 : "GSL Interpolation only works with double" );
64 :
65 0 : gsl_set_error_handler_off();
66 :
67 0 : gsl_interp *interp = gsl_interp_alloc( interpT, Nin );
68 0 : if( interp == nullptr )
69 : {
70 0 : std::cerr << "gsl_interpolate: gsl_interp_alloc failed\n";
71 0 : return -1;
72 : }
73 :
74 0 : gsl_interp_accel *acc = gsl_interp_accel_alloc();
75 0 : if( acc == nullptr )
76 : {
77 0 : std::cerr << "gsl_interpolate: gsl_interp_accel_alloc failed\n";
78 0 : return -1;
79 : }
80 0 : int gsl_errno = gsl_interp_init( interp, xin, yin, Nin );
81 :
82 0 : if( gsl_errno != 0 )
83 : {
84 0 : std::cerr << "gsl_interpolate: error from gsl_interp_init [" << gsl_strerror( gsl_errno ) << "]\n";
85 0 : gsl_interp_free( interp );
86 0 : gsl_interp_accel_free( acc );
87 0 : return -1;
88 : }
89 :
90 0 : gsl_errno = 0;
91 0 : for( size_t i = 0; i < Nout; ++i )
92 : {
93 : // Don't error check, let it set NAN
94 0 : gsl_errno += gsl_interp_eval_e( interp, xin, yin, xout[i], acc, &yout[i] );
95 : }
96 :
97 0 : gsl_interp_free( interp );
98 0 : gsl_interp_accel_free( acc );
99 :
100 0 : if( gsl_errno )
101 : {
102 0 : std::cerr << "gsl_interpolate: error(s) reported by gsl_interp_eval_e\n";
103 0 : return -1;
104 : }
105 :
106 0 : return 0;
107 : }
108 :
109 : /// Interpolate a 1-D data X vs Y discrete function onto a new X axis (vector version)
110 : /**
111 : * \param interpT of the <a href="https://www.gnu.org/software/gsl/doc/html/interp.html#d-interpolation-types">gsl
112 : * interpolation types</a>. \param [in] xin the input x-axis \param [in] yin the input y-values \param [in] xout the
113 : * desired x-axis \param [out] yout the output interpolated y-values, does not need to be allocated
114 : *
115 : * \returns 0 on success or -1 if GSL reports an interpolation error.
116 : *
117 : * \ingroup interpolation
118 : */
119 : template <typename realT>
120 0 : int gsl_interpolate( const gsl_interp_type *interpT,
121 : std::vector<realT> &xin,
122 : std::vector<realT> &yin,
123 : std::vector<realT> &xout,
124 : std::vector<realT> &yout )
125 : {
126 0 : yout.resize( xout.size() );
127 0 : return gsl_interpolate( interpT, xin.data(), yin.data(), xin.size(), xout.data(), yout.data(), xout.size() );
128 : }
129 :
130 : } // namespace math
131 : } // namespace mx
132 :
133 : #endif // gslInterpolation_hpp
|
