diff options
Diffstat (limited to 'mass_profile/spline.hpp')
-rw-r--r-- | mass_profile/spline.hpp | 110 |
1 files changed, 0 insertions, 110 deletions
diff --git a/mass_profile/spline.hpp b/mass_profile/spline.hpp deleted file mode 100644 index 1027b61..0000000 --- a/mass_profile/spline.hpp +++ /dev/null @@ -1,110 +0,0 @@ -#ifndef SPLINE_HPP -#define SPLINE_HPP - -#include <vector> -#include <cstdlib> -#include <cassert> -#include <cmath> -#include <limits> - -template <typename T> -class spline -{ -public: - std::vector<T> x_list; - std::vector<T> y_list; - std::vector<T> y2_list; - -public: - void push_point(T x,T y) - { - if(!x_list.empty()) - { - assert(x>*(x_list.end()-1)); - } - x_list.push_back(x); - y_list.push_back(y); - } - - T get_value(T x) - { - if(x<=x_list[0]) - { - return y_list[0]; - } - if(x>=x_list.back()) - { - return y_list.back(); - } - assert(x_list.size()==y2_list.size()); - assert(x>x_list[0]); - assert(x<x_list.back()); - int n1,n2; - n1=0; - n2=x_list.size()-1; - while((n2-n1)!=1) - { - //cerr<<n1<<"\t"<<n2<<endl; - if(x_list[n1+1]<=x) - { - n1++; - } - if(x_list[n2-1]>x) - { - n2--; - } - } - T h=x_list[n2]-x_list[n1]; - double a=(x_list[n2]-x)/h; - double b=(x-x_list[n1])/h; - return a*y_list[n1]+b*y_list[n2]+((a*a*a-a)*y2_list[n1]+ - (b*b*b-b)*y2_list[n2])*(h*h)/6.; - - } - - void gen_spline(T y2_0,T y2_N) - { - int n=x_list.size(); - y2_list.resize(0); - y2_list.resize(x_list.size()); - std::vector<T> u(x_list.size()); - if(std::abs(y2_0)<std::numeric_limits<T>::epsilon()) - { - y2_list[0]=0; - u[0]=0; - } - else - { - y2_list[0]=-.5; - u[0]=(3./(x_list[1]-x_list[0]))*((y_list[1]-y_list[0])/(x_list[1]-x_list[0])-y2_0); - } - for(int i=1;i<n-1;++i) - { - double sig=(x_list[i]-x_list[i-1])/(x_list[i+1]-x_list[i-1]); - double p=sig*y2_list[i-1]+2.; - y2_list[i]=(sig-1.)/p; - u[i]=(y_list[i+1]-y_list[i])/(x_list[i+1]-x_list[i]) - -(y_list[i]-y_list[i-1])/(x_list[i]-x_list[i-1]); - u[i]=(6.*u[i]/(x_list[i+1]-x_list[i-1])-sig*u[i-1])/p; - } - double qn,un; - if(std::abs(y2_N)<std::numeric_limits<T>::epsilon()) - { - qn=un=0; - } - else - { - qn=.5; - un=(3./(x_list[n-1]-x_list[n-2]))*(y2_N-(y_list[n-1]-y_list[n-2])/(x_list[n-1]-x_list[n-2])); - - } - y2_list[n-1]=(un-qn*u[n-2])/(qn*y2_list[n-2]+1.); - for(int i=n-2;i>=0;--i) - { - y2_list[i]=y2_list[i]*y2_list[i+1]+u[i]; - } - } - -}; - -#endif /* SPLINE_HPP */ |