gsl_interpolator.hpp

Go to the documentation of this file.

/*
 * This file is part of BubbleProfiler.
 *
 * BubbleProfiler 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.
 *
 * BubbleProfiler 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 BubbleProfiler.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef BUBBLEPROFILER_GSL_INTERPOLATOR_HPP_INCLUDED
#define BUBBLEPROFILER_GSL_INTERPOLATOR_HPP_INCLUDED

#include <Eigen/Core>

#include <gsl/gsl_spline.h>
#include <gsl/gsl_version.h>

#include <memory>

namespace BubbleProfiler {

Eigen::VectorXd interpolate_f_at(const Eigen::VectorXd& x_old,
                                 const Eigen::VectorXd& y_old,
                                 const Eigen::VectorXd& x_new);

Eigen::VectorXd interpolate_f_at(const Eigen::VectorXd& x_old,
                                 const Eigen::VectorXd& y_old,
                                 int n_old, const Eigen::VectorXd& x_new,
                                 int n_new);

class GSL_interpolator {
public:
   enum class Interpolation_type {
      GSL_linear,
      GSL_polynomial,
      GSL_cspline,
      GSL_cspline_periodic,
      GSL_akima,
      GSL_akima_periodic
#ifdef GSL_MAJOR_VERSION
#if (GSL_MAJOR_VERSION >= 2)
      , GSL_steffen
#endif
#endif
   };

   GSL_interpolator();
   ~GSL_interpolator() = default;
   GSL_interpolator(const GSL_interpolator&);
   GSL_interpolator& operator=(const GSL_interpolator&);
   GSL_interpolator(GSL_interpolator&&) = default;
   GSL_interpolator& operator=(GSL_interpolator&&) = default;

   void set_interpolation_type(Interpolation_type t) { interp_type = t; }

   void construct_interpolant(const Eigen::VectorXd&, const Eigen::VectorXd&);

   void construct_interpolant(const Eigen::VectorXd&, const Eigen::VectorXd&,
                              int);

   double evaluate_at(double) const;

   double evaluate_deriv_at(double) const;

   double evaluate_second_deriv_at(double) const;

private:
   struct GSL_accel_deleter {
      void operator()(gsl_interp_accel* a) const {
         if (a) {
            gsl_interp_accel_free(a);
         }
      }
   };

   struct GSL_spline_deleter {
      void operator()(gsl_spline* s) const {
         if (s) {
            gsl_spline_free(s);
         }
      }
   };

   Interpolation_type interp_type{Interpolation_type::GSL_cspline};
   std::unique_ptr<gsl_interp_accel, GSL_accel_deleter> accelerator{nullptr};
   std::unique_ptr<gsl_spline, GSL_spline_deleter> interpolant{nullptr};

   const gsl_interp_type* to_gsl_interp_type(Interpolation_type) const;
   void reset(Interpolation_type, int);
   void check_interpolant_computed() const;
   void check_within_bounds(double) const;
};

} // namespace BubbleProfiler

#endif