GridCubicSpline.h

//
//  Copyright (C) 2017, Alexander Stukowski and Paul Erhart
//
//  This file is part of atomicrex.
//
//  Atomicrex 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.
//
//  Atomicrex 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/>.
//

#pragma once

#include "../../Atomicrex.h"

namespace atomicrex {

class GridCubicSpline
{
public:

    GridCubicSpline() : _N(0), _delta(0) {}

    void init(int N, double delta, double derivStart = 0, double derivEnd = 0) {
        BOOST_ASSERT(N >= 2);
        BOOST_ASSERT(delta > 0);
        _N = N;
        _derivStart = derivStart;
        _derivEnd = derivEnd;
        _delta = delta;
        _xcutoff = delta * (N-1);
        Y.resize(N);
        Y2.resize(N);
    }

    void setKnot(int i, double y) {
        BOOST_ASSERT(i >=0 && i < _N);
        Y[i] = y;
    }

    int knotCount() const { return _N; }

    double derivStart() const { return _derivStart; }

    double derivEnd() const { return _derivEnd; }

    double knotX(int i) const { BOOST_ASSERT(i >=0 && i < _N); return _delta * i; }

    double knotY(int i) const { BOOST_ASSERT(i >=0 && i < _N); return Y[i]; }

    double knotY2(int i) const { BOOST_ASSERT(i >=0 && i < _N); return Y2[i]; }

    bool parse(std::istream& stream);

    void prepareSpline();

    double eval(double x) const;

    double eval(double x, double& deriv) const;

    double cutoff() const { return _xcutoff; }

    void writeGnuplot(const FPString& filename, const FPString& title) const;

private:

    int _N;                     // Number of spline knots
    double _delta;              // Distance between spline knots.
    double _xcutoff;            // The X coordinate of the last knot.
    std::vector<double> Y;      // Function values at spline knots
    std::vector<double> Y2;     // Second derivatives at spline knots
    double _derivStart;         // First derivative at knot 0
    double _derivEnd;           // First derivative at knot (N-1)
};

} // End of namespace