Commit e31b7899 authored by Cédric Mezrag's avatar Cédric Mezrag
Browse files

refs#16

Adding simple cubic spline algorithm for NumA++
parent 461e31a0
/**
* @file CubicSpline.h
* @author Cédric Mezrag (ANL)
* @date 14 April 2017
* @version 0.1
*/
#ifndef CUBICSPLINE_H_
#define CUBICSPLINE_H_
#include <iostream>
#include <sstream>
#include <vector>
#include <cmath>
/*
* Cubic spline algorithm based on wikipedia :
* https://en.wikipedia.org/w/index.php?title=Spline_%28mathematics%29&oldid=288288033#Algorithm_for_computing_natural_cubic_splines
*
*
*/
namespace NumA{
class CubicSpline{
public:
CubicSpline(std::vector<double> x, std::vector<double> y);
virtual ~CubicSpline();
void ConstructSpline(); // construct natural cubic splines
double getSplineInsideValue(double z); //Return natural cubic splines for xmin < z < xmax and 0 if z is outside the set of point.
private:
const unsigned int m_N ; // size of the point sample
std::vector<double> m_X ; // x coordinate of the sample
std::vector<double> m_Y ; // y coordinate of the sample
std::vector<double> m_aCoeff ; // coefficient a as defined on wikipedia
std::vector<double> m_bCoeff ; // coefficient b as defined on wikipedia
std::vector<double> m_cCoeff ; // coefficient c as defined on wikipedia
std::vector<double> m_dCoeff ; // coefficient d as defined on wikipedia
bool m_SplineDefined ; // Turn to true if ConstructSpline is run at least once.
};
}
#endif
#include "../../../include/NumA/interpolation/CubicSpline.h"
#include <iostream>
#include <sstream>
#include <vector>
#include <cmath>
/*
* Cubic spline algorithm based on wikipedia :
* https://en.wikipedia.org/w/index.php?title=Spline_%28mathematics%29&oldid=288288033#Algorithm_for_computing_natural_cubic_splines
*
*
*/
namespace NumA {
CubicSpline::CubicSpline(std::vector<double> x, std::vector<double> y):
m_X(x),
m_Y(y),
m_N(x.size()-1),
m_aCoeff(y),
m_bCoeff(0),
m_cCoeff(0),
m_dCoeff(0),
m_SplineDefined(false)
{
for(unsigned int i = 0; i <= m_N; i++){
m_bCoeff.push_back(0.) ;
m_cCoeff.push_back(0.) ;
m_dCoeff.push_back(0.) ;
};
}
CubicSpline::~CubicSpline(){}
void CubicSpline::ConstructSpline(){
std::vector<double> h(0) ;
std::vector<double> alpha(0) ;
std::vector<double> l(0) ;
std::vector<double> mu(0) ;
std::vector<double> z(0) ;
for(unsigned int i = 0; i < m_N; i++){
h.push_back(m_X.at(i+1)-m_X.at(i));
};
alpha.push_back(0.) ;
for(unsigned int i=1; i < m_N; i++){
alpha.push_back( 3. / h.at(i) * ( m_aCoeff.at(i+1) - m_aCoeff.at(i) ) - 3. / h.at(i-1) * (m_aCoeff.at(i) - m_aCoeff.at(i-1)) );
};
//std::cout << " alpha is ok" << std::endl ;
l.push_back(1.);
mu.push_back(0.);
z.push_back(0.);
for(unsigned int i=1 ; i< m_N ; i++){
l.push_back( 2. * (m_X.at(i+1)-m_X.at(i-1) ) - h.at(i-1) * mu.at(i-1) );
mu.push_back( h.at(i) / l.at(i) );
z.push_back( ( alpha.at(i) - h.at(i-1) * z.at( i-1 ) ) / l.at(i) );
} ;
l.push_back(1) ;
z.push_back(0) ;
m_cCoeff.at(m_N) = 0. ;
for(int i = m_N - 1; i >= 0; i--){
m_cCoeff.at(i) = z.at(i) - mu.at(i)* m_cCoeff.at(i+1);
m_bCoeff.at(i) = ( m_aCoeff.at(i+1) - m_aCoeff.at(i) ) / h.at(i)
- 1./3. * h.at(i) * ( m_cCoeff.at(i+1) + 2 * m_cCoeff.at(i) ) ;
m_dCoeff.at(i) = ( m_cCoeff.at(i+1) - m_cCoeff.at(i) ) / 3. / h.at(i) ;
} ;
m_SplineDefined = true ;
}
double CubicSpline::getSplineInsideValue(double z){
// Test if splines have been computed before
if(m_SplineDefined){
if(z < m_X.at(0) or z > m_X.at(m_N)){
std::cout << "Warning X is outside the interpolation domain " << std::endl ;
std::cout << " Interpolation wil return zero " << std::endl ;
return 0. ;
}
else{
for(unsigned int i = 0 ; i < m_N ; i++ ){
if(z <= m_X.at(i+1) ){
return m_aCoeff.at(i) + m_bCoeff.at(i)*(z-m_X.at(i)) + m_cCoeff.at(i)* pow( (z-m_X.at(i)) , 2.)
+ m_dCoeff.at(i) * pow( (z-m_X.at(i)) , 3.) ;
};
};
};
}
else{
std::cout << "Warning splines have not been defined " << std::endl;
std::cout << " Return is set to zero " << std::endl;
return 0. ;
};
std::cout << " If this point is reached, there is a problem in computing the spline" << std::endl ;
}
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment