spline.h

/*************************************************************************************
 * MechSys - A C++ library to simulate (Continuum) Mechanical Systems                *
 * Copyright (C) 2005 Dorival de Moraes Pedroso <dorival.pedroso at gmail.com>       *
 * Copyright (C) 2005 Raul Dario Durand Farfan  <raul.durand at gmail.com>           *
 *                                                                                   *
 * This file is part of MechSys.                                                     *
 *                                                                                   *
 * MechSys 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.                                                                          *
 *                                                                                   *
 * MechSys 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      *
 * MechSys; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, *
 * Fifth Floor, Boston, MA 02110-1301, USA                                           *
 *************************************************************************************/

#ifndef MECHSYS_SPLINE_H
#define MECHSYS_SPLINE_H

#ifdef HAVE_CONFIG_H
  #include "config.h"
#else
  #ifndef REAL
    #define REAL double
  #endif
#endif

#define H_ZERO 1.0e-4

#include "util/exception.h"
#include "util/numstreams.h"

using Util::_4;
using Util::_6_3;

using std::cout;
using std::endl;

class Spline
{
    friend std::ostream & operator<< (std::ostream & os, Spline const & S);
public:
    // Constructor & Destructor
     Spline (REAL const X[], REAL const Y[], int Size);
    ~Spline ();
    // Methods
    REAL Eval       (REAL x) const;
    bool CheckRange (REAL x) const;
    REAL YRange     ()       const;
private:
    int    _n;
    REAL * _x; // size == _n
    REAL * _a; // size == _n
    REAL * _b; // size == _n - 1
    REAL * _c; // size == _n
    REAL * _d; // size == _n - 1
}; // class Spline




inline Spline::Spline(REAL const X[], REAL const Y[], int Size) // {{{
    : _n (Size)
{
    // Check
    if (_n<2) throw new Fatal(_("Spline::Spline: The size or arrays (%d) must be 2 or greater"), _n);

    // Allocate internal arrays
    _x = new REAL [_n];
    _a = new REAL [_n];
    _b = new REAL [_n-1];
    _c = new REAL [_n];
    _d = new REAL [_n-1];

    // Initialize arrays (x must be increasing)
    if (X[_n-1]<X[0]) // X decreasing
    {
        for (int i=0; i<_n; ++i)
        {
            _x[_n-1-i] = X[i];
            _a[_n-1-i] = Y[i];
        }
    }
    else
    {
        for (int i=0; i<_n; ++i)
        {
            _x[i] = X[i];
            _a[i] = Y[i];
        }
    }

    // Allocate temporary arrays
    REAL * h = new REAL [_n-1];
    REAL * A = new REAL [_n-1];
    REAL * l = new REAL [_n  ];
    REAL * m = new REAL [_n-1];
    REAL * z = new REAL [_n  ];

    // Calculate alpha (A)
    A[0] = 0.0; // not used
    for (int i=0; i<_n-1; ++i)
    {
        h[i] = _x[i+1] - _x[i];
        if (fabs(h[i])<H_ZERO) throw new Fatal(_("Spline::Spline: h == 0"));
        if (i>0)
            A[i] = 3.0*(_a[i+1]-_a[i])/h[i] - 3.0*(_a[i]-_a[i-1])/h[i-1];
    }

    // Solve tridiagonal system by Crout Factorization method
    l[0] = 1.0;
    m[0] = 0.0;
    z[0] = 0.0;
    for (int i=1; i<_n-1; ++i)
    {
        l[i] = 2.0*(_x[i+1]-_x[i-1]) - h[i-1]*m[i-1];
        m[i] = h[i]/l[i];
        z[i] = (A[i]-h[i-1]*z[i-1])/l[i];
    }
     l[_n-1] = 1.0;
     z[_n-1] = 0.0;
    _c[_n-1] = 0.0;
    for (int i=_n-2; i>=0; --i)
    {
        _c[i] = z[i] - m[i]*_c[i+1];
        _b[i] = (_a[i+1]-_a[i])/h[i] - h[i]*(_c[i+1]+2.0*_c[i])/3.0;
        _d[i] = (_c[i+1]-_c[i])/(3.0*h[i]);
    }

    // Clean up
    delete [] h;
    delete [] A;
    delete [] l;
    delete [] m;
    delete [] z;

} // }}}

inline Spline::~Spline() // {{{
{
    delete [] _x;
    delete [] _a;
    delete [] _b;
    delete [] _c;
    delete [] _d;
} // }}}

inline REAL Spline::Eval(REAL x) const // {{{
{
    // Find bounds for x
    int k;
    int lo = 0;
    int hi = _n-1;
    while (hi-lo>1)
    {
        k = (hi+lo) >> 1;
        if (_x[k]>x) hi = k;
        else         lo = k;
    }

    // Check
    //std::cout << "x_lo=" << _x[lo] << "  x=" << x << "  x_hi=" << _x[hi] << std::endl;
    if (x<_x[lo]) throw new Fatal(_("Spline::Eval: Value (x=%.6f < %.6f) is not valid"), x, _x[lo]);
    if (x>_x[hi]) throw new Fatal(_("Spline::Eval: Value (x=%.6f > %.6f) is not valid"), x, _x[hi]);

    // Interpolate
    REAL h = x-_x[lo];
    return _a[lo] + _b[lo]*h + _c[lo]*h*h + _d[lo]*h*h*h;

} // }}}

inline bool Spline::CheckRange(REAL x) const // {{{
{
    if (x<_x[0]   ) return false;
    if (x>_x[_n-1]) return false;
    return true;
} // }}}

inline REAL Spline::YRange() const // {{{
{
    return _a[_n-1]-_a[0];
} // }}}

std::ostream & operator<< (std::ostream & os, Spline const & S) // {{{
{
    os << _4()<< "i" << _6_3()<< "x" << _6_3()<< "a" << _6_3()<< "b" << _6_3()<< "c" << _6_3()<< "d" << std::endl;
    for (int i=0; i<S._n; ++i)
    {
        os << _4()<< i << _6_3()<< S._x[i] << _6_3()<< S._a[i];
        if (i<S._n-1) os << _6_3()<< S._b[i] << _6_3()<< S._c[i] << _6_3()<< S._d[i] << std::endl;
        else os << std::endl;
    }
    return os;
} // }}}

#endif // MECHSYS_SPLINE_H

// vim:fdm=marker

---