/*
 * spline.h
 *
 * simple cubic spline interpolation library without external
 * dependencies
 *
 * ---------------------------------------------------------------------
 * Copyright (C) 2011, 2014 Tino Kluge (ttk448 at gmail.com)
 *
 *  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, see <http://www.gnu.org/licenses/>.
 * ---------------------------------------------------------------------
 *
 */


#ifndef TK_SPLINE_H
#define TK_SPLINE_H

#include <cstdio>
//#include <cassert>
#include <vector>
#include <algorithm>


namespace tk
{

  // band matrix solver
  class band_matrix
  {
  private:
      std::vector< std::vector<double> > m_upper;  // upper band
      std::vector< std::vector<double> > m_lower;  // lower band
  public:
      band_matrix() {};                             // constructor
      band_matrix(int dim, int n_u, int n_l);       // constructor
      ~band_matrix() {};                            // destructor
      bool resize(int dim, int n_u, int n_l);      // init with dim,n_u,n_l
      int dim() const;                             // matrix dimension
      int num_upper() const
      {
          return m_upper.size()-1;
      }
      int num_lower() const
      {
          return m_lower.size()-1;
      }
      // access operator
      double & operator () (int i, int j);            // write
      double   operator () (int i, int j) const;      // read
      // we can store an additional diogonal (in m_lower)
      double& saved_diag(int i);
      double  saved_diag(int i) const;
      void lu_decompose();
      std::vector<double> r_solve(const std::vector<double>& b);
      std::vector<double> l_solve(const std::vector<double>& b);
      std::vector<double> lu_solve(const std::vector<double>& b,
                                   bool is_lu_decomposed=false);
      bool success;

  };


  // spline interpolation
  class spline
  {
  public:
      enum bd_type {
          first_deriv = 1,
          second_deriv = 2
      };

  private:
      std::vector<double> m_x,m_y;            // x,y coordinates of points
      // interpolation parameters
      // f(x) = a*(x-x_i)^3 + b*(x-x_i)^2 + c*(x-x_i) + y_i
      std::vector<double> m_a,m_b,m_c;        // spline coefficients
      double  m_b0, m_c0;                     // for left extrapol
      bd_type m_left, m_right;
      double  m_left_value, m_right_value;
      bool    m_force_linear_extrapolation;

  public:
      // set default boundary condition to be zero curvature at both ends
      spline(): m_left(second_deriv), m_right(second_deriv),
          m_left_value(0.0), m_right_value(0.0),
          m_force_linear_extrapolation(false)
      {
          ;
      }

      // optional, but if called it has to come be before set_points()
      void set_boundary(bd_type left, double left_value,
                        bd_type right, double right_value,
                        bool force_linear_extrapolation=false);
      bool set_points(const std::vector<double>& x,
                      const std::vector<double>& y, bool cubic_spline=true);
      double operator() (double x) const;
  };

} // namespace tk

#endif /* TK_SPLINE_H */
