#ifndef
interpolation_header
#define
interpolation_header
#include
#include
#include
#include
using
namespace std;
#define
delta 1e-6
#define
very_close(x,y) fabs((x)-(y))
template
class
point2d
t x,y;
};
bool
myfunc(const point2d
class
point2d_x_less_delta
};
template
ostream & operator
<<(ostream &os,const point2d& p)
template
istream & operator >>(ostream &is,point2d& p)
class
interpolation**ooth
template
interpolation**ooth(inputiterator iter1,inputiterator iter2):
pts(iter1,iter2)
void insert(const point2d
template
void insert(inputiterator iter1,inputiterator iter2)
double getvalue(double x)
setdouble>,point2d_x_less_delta >::iterator iter=lower_bound(pts.begin(),pts.end(),point2d
if (iter==pts.begin())
int k=distance(pts.begin(),iter);
k--;
double xkp1=(*iter).x;
iter--;
double s0=(*iter).y;
double s1=g(k);
double s2=(u(k)*3-g(k)*2-g(k+1))/(xkp1-(*iter).x);
double s3=(g(k)+g(k+1)-u(k)*2.0)/((xkp1-(*iter).x)*(xkp1-(*iter).x));
s1*=x-(*iter).x;
s2*=(x-(*iter).x)*(x-(*iter).x);
s3*=(x-(*iter).x)*(x-(*iter).x)*(x-(*iter).x);
return s0+s1+s2+s3;
}
double getderivativevalue(double x)
setdouble>,point2d_x_less_delta >::iterator iter=lower_bound(pts.begin(),pts.end(),point2d
if (iter==pts.begin())
int k=distance(pts.begin(),iter);
k--;
double xkp1=(*iter).x;
iter--;
double s1=g(k);
double s2=(u(k)*3-g(k)*2-g(k+1))/(xkp1-(*iter).x);
double s3=(g(k)+g(k+1)-u(k)*2.0)/((xkp1-(*iter).x)*(xkp1-(*iter).x));
s2*=(x-(*iter).x)*2.0;
s3*=(x-(*iter).x)*(x-(*iter).x)*3.0;
return s1+s2+s3;
}
setdouble>,point2d_x_less_delta >::iterator begin()
setdouble>,point2d_x_less_delta >::iterator end()
size_t size()
private:
double u(int _pos)
if (_pos<0)
setdouble>,point2d_x_less_delta >::iterator iter1=pts.begin();
setdouble>,point2d_x_less_delta >::iterator iter2=pts.begin();
advance(iter1,_pos);
advance(iter2,_pos+1);
return ((*iter2).y-(*iter1).y)/((*iter2).x-(*iter1).x);
}
double g(int _pos)
return (fabs(ukp1-uk)*ukm1+fabs(ukm1-ukm2)*
uk
)/(fabs(ukp1-uk)+fabs(ukm1-ukm2));
}
setdouble>,point2d_x_less_delta > pts;
};
#endif