src/nr/spline.cpp

Go to the documentation of this file.

//
// C++ Implementation: spline
//
// Description: 
//
//
// Author: Roman Geus <geus@maxwell>, (C) 2004
//
// Copyright: See COPYING file that comes with this distribution
//
//

#include <cassert>
#include "spline.h"

namespace NR {

using namespace std;

void spline(const vector<double>& x,
            const vector<double>& y,
            const double yp1,
            const double ypn,
            vector<double>& y2)
{
    const double y_thresh = 0.99e30;

    unsigned int n = x.size();
    assert(y.size() == n);
    y2.resize(n);

    double p, qn, sig, un;
    vector<double> u(n-1);
    if (yp1 > y_thresh)
        y2[0] = u[0] = 0.0;
    else {
        y2[0] = -0.5;
        u[0] = (3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1);
    }
    for (unsigned int i = 1; i < n-1; ++ i) {
        sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]);
        p = sig*y2[i-1] + 2.0;
        y2[i] = (sig - 1.0)/p;
        u[i] = (y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
        u[i] = (6.0*u[i]/(x[i+1]-x[i-1]) - sig*u[i-1])/p;
    }
    if (ypn > y_thresh)
        qn = un = 0.0;
    else {
        qn = -0.5;
        un = (3.0/(x[n-1]-x[n-2]))*(ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
    }
    y2[n-1] = (un-qn*u[n-2])/(qn*y2[n-2]+1.0);
    for (int k = n-2; k >= 0; -- k)
        y2[k] = y2[k]*y2[k+1] + u[k];
}

void splint(const vector<double>& xa,
            const vector<double>& ya,
            const vector<double>& y2a,
            const double x,
            double& y)
{
    int k;
    double h, b, a;

    unsigned int n = xa.size();
    assert(ya.size() == n);
    assert(y2a.size() == n);

    // Find table index using bisection
    int klo = 0;
    int khi = n-1;
    while (khi - klo > 1) {
        k = (khi + klo) >> 1;
        if (xa[k] > x)
            khi = k;
        else
            klo = k;
    }

    h = xa[khi]-xa[klo];
    assert(h > 0.0);
    a = (xa[khi] - x)/h;
    b = (x - xa[klo])/h;
    y = a*ya[klo] + b*ya[khi] + ((a*a*a - a)*y2a[klo] + (b*b*b - b)*y2a[khi])*(h*h)/6.0;
}

} // namespace NR

---