46 #ifndef O2SCL_INTERP_H 47 #define O2SCL_INTERP_H 56 #include <boost/numeric/ublas/vector.hpp> 57 #include <boost/numeric/ublas/vector_proxy.hpp> 59 #include <o2scl/search_vec.h> 60 #include <o2scl/tridiag.h> 61 #include <o2scl/vector.h> 63 #ifndef DOXYGEN_NO_O2NS 109 #ifdef O2SCL_NEVER_DEFINED 113 #ifndef DOXYGEN_INTERNAL 136 double integ_eval(
double ai,
double bi,
double ci,
double di,
double xi,
137 double a,
double b)
const {
142 double bterm=0.5*bi*r12;
143 double cterm=ci*(r1*r1+r2*r2+r1*r2)/3.0;
144 double dterm=0.25*di*r12*(r1*r1+r2*r2);
146 return (b-a)*(ai+bterm+cterm+dterm);
168 virtual void set(
size_t size,
const vec_t &x,
const vec2_t &y)=0;
171 virtual double eval(
double x0)
const=0;
179 virtual double deriv(
double x0)
const=0;
184 virtual double deriv2(
double x0)
const=0;
187 virtual double integ(
double a,
double b)
const=0;
190 virtual const char *
type()
const=0;
192 #ifndef DOXYGEN_INTERNAL 213 #ifdef O2SCL_NEVER_DEFINED 226 virtual void set(
size_t size,
const vec_t &x,
const vec2_t &y) {
229 " than "+
szttos(this->min_size)+
" in interp_linear::"+
240 virtual double eval(
double x0)
const {
243 size_t index=this->
svx.find_const(x0,cache);
245 double x_lo=(*this->
px)[index];
246 double x_hi=(*this->
px)[index+1];
247 double y_lo=(*this->
py)[index];
248 double y_hi=(*this->
py)[index+1];
251 return y_lo+(x0-x_lo)/dx*(y_hi-y_lo);
255 virtual double deriv(
double x0)
const {
258 size_t index=this->
svx.find_const(x0,cache);
260 double x_lo=(*this->
px)[index];
261 double x_hi=(*this->
px)[index+1];
262 double y_lo=(*this->
py)[index];
263 double y_hi=(*this->
py)[index+1];
278 virtual double integ(
double a,
double b)
const {
281 size_t i, index_a, index_b;
284 if (((*this->
px)[0]<(*this->
px)[this->
sz-1] && a>b) ||
285 ((*this->
px)[0]>(*this->
px)[this->
sz-1] && a<b)) {
292 index_a=this->
svx.find_const(a,cache);
293 index_b=this->
svx.find_const(b,cache);
296 for(i=index_a; i<=index_b; i++) {
298 double x_lo=(*this->
px)[i];
299 double x_hi=(*this->
px)[i+1];
300 double y_lo=(*this->
py)[i];
301 double y_hi=(*this->
py)[i+1];
306 if (i == index_a || i == index_b) {
307 double x1=(i == index_a) ? a : x_lo;
308 double x2=(i == index_b) ? b : x_hi;
309 double D=(y_hi-y_lo)/dx;
310 result += (x2-
x1)*(y_lo+0.5*D*((x2-x_lo)+(x1-x_lo)));
312 result += 0.5*dx*(y_lo+y_hi);
316 std::string str=((std::string)
"Interval of length zero ")+
319 " in interp_linear::integ().";
324 if (flip) result=-result;
329 virtual const char *
type()
const {
return "interp_linear"; }
331 #ifndef DOXYGEN_INTERNAL 350 #ifdef O2SCL_NEVER_DEFINED 363 virtual void set(
size_t size,
const vec_t &x,
const vec2_t &y) {
366 " than "+
szttos(this->min_size)+
" in interp_nearest_neigh::"+
377 virtual double eval(
double x0)
const {
380 return (*this->
py)[index];
384 virtual double deriv(
double x0)
const {
396 virtual double integ(
double a,
double b)
const {
401 virtual const char *
type()
const {
return "interp_nearest_neigh"; }
403 #ifndef DOXYGEN_INTERNAL 410 (
const interp_nearest_neigh<vec_t,vec2_t>&);
428 #ifdef O2SCL_NEVER_DEFINED 435 typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
436 typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
437 typedef boost::numeric::ublas::slice slice;
438 typedef boost::numeric::ublas::range range;
440 #ifndef DOXYGEN_INTERNAL 456 void coeff_calc(
const ubvector &c_array,
double dy,
double dx,
457 size_t index,
double &b,
double &c2,
double &d)
const {
459 double c_i=c_array[index];
460 double c_ip1=c_array[index+1];
461 b=(dy/dx)-dx*(c_ip1+2.0*c_i)/3.0;
463 d=(c_ip1-c_i)/(3.0*dx);
484 virtual void set(
size_t size,
const vec_t &xa,
const vec2_t &ya) {
488 " than "+
szttos(this->min_size)+
" in interp_cspline::"+
492 if (size!=this->
sz) {
496 offdiag.resize(size);
509 size_t num_points=size;
510 size_t max_index=num_points-1;
511 size_t sys_size=max_index-1;
516 for (i=0; i < sys_size; i++) {
517 double h_i=xa[i+1]-xa[i];
518 double h_ip1=xa[i+2]-xa[i+1];
519 double ydiff_i=ya[i+1]-ya[i];
520 double ydiff_ip1=ya[i+2]-ya[i+1];
521 double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
522 double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
524 diag[i]=2.0*(h_ip1+h_i);
525 g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
535 ubvector_range cp1(c,range(1,c.size()));
538 (diag,offdiag,g,cp1,sys_size,p4m);
544 virtual double eval(
double x0)
const {
547 size_t index=this->
svx.find_const(x0,cache);
549 double x_lo=(*this->
px)[index];
550 double x_hi=(*this->
px)[index+1];
553 double y_lo=(*this->
py)[index];
554 double y_hi=(*this->
py)[index+1];
557 double b_i, c_i, d_i;
559 coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
561 return y_lo+delx*(b_i+delx*(c_i+delx*d_i));
565 virtual double deriv(
double x0)
const {
568 size_t index=this->
svx.find_const(x0,cache);
570 double x_lo=(*this->
px)[index];
571 double x_hi=(*this->
px)[index+1];
574 double y_lo=(*this->
py)[index];
575 double y_hi=(*this->
py)[index+1];
578 double b_i, c_i, d_i;
580 coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
582 return b_i+delx*(2.0*c_i+3.0*d_i*delx);
591 size_t index=this->
svx.find_const(x0,cache);
593 double x_lo=(*this->
px)[index];
594 double x_hi=(*this->
px)[index+1];
597 double y_lo=(*this->
py)[index];
598 double y_hi=(*this->
py)[index+1];
601 double b_i, c_i, d_i;
603 coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
605 return 2.0*c_i+6.0*d_i*delx;
609 virtual double integ(
double a,
double b)
const {
611 size_t i, index_a, index_b;
614 if (((*this->
px)[0]<(*this->
px)[this->
sz-1] && a>b) ||
615 ((*this->
px)[0]>(*this->
px)[this->
sz-1] && a<b)) {
623 index_a=this->
svx.find_const(a,cache);
624 index_b=this->
svx.find_const(b,cache);
628 for(i=index_a; i<=index_b; i++) {
630 double x_lo=(*this->
px)[i];
631 double x_hi=(*this->
px)[i+1];
632 double y_lo=(*this->
py)[i];
633 double y_hi=(*this->
py)[i+1];
638 double b_i, c_i, d_i;
639 coeff_calc(c,dy,dx,i,b_i,c_i,d_i);
640 if (i == index_a || i == index_b) {
641 double x1=(i == index_a) ? a : x_lo;
642 double x2=(i == index_b) ? b : x_hi;
643 result += this->
integ_eval(y_lo,b_i,c_i,d_i,x_lo,x1,x2);
645 result += dx*(y_lo+dx*(0.5*b_i+
646 dx*(c_i/3.0+0.25*d_i*dx)));
649 std::string str=((std::string)
"Interval of length zero ")+
652 " in interp_cspline::integ().";
658 if (flip) result*=-1.0;
664 virtual const char *
type()
const {
return "interp_cspline"; }
666 #ifndef DOXYGEN_INTERNAL 672 (
const interp_cspline<vec_t,vec2_t>&);
686 #ifdef O2SCL_NEVER_DEFINED 698 typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
699 typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
700 typedef boost::numeric::ublas::slice slice;
701 typedef boost::numeric::ublas::range range;
711 virtual const char *
type()
const {
return "interp_cspline_peri"; }
715 virtual void set(
size_t size,
const vec_t &xa,
const vec2_t &ya) {
719 " than "+
szttos(this->min_size)+
" in interp_cspline"+
723 if (size!=this->
sz) {
724 this->c.resize(size);
725 this->g.resize(size);
726 this->diag.resize(size);
727 this->offdiag.resize(size);
740 size_t num_points=size;
742 size_t max_index=num_points-1;
744 size_t sys_size=max_index;
750 double h0=xa[1]-xa[0];
751 double h1=xa[2]-xa[1];
752 double h2=xa[3]-xa[2];
753 double A=2.0*(h0+h1);
758 gx[0]=3.0*((ya[2]-ya[1])/h1-(ya[1]-ya[0])/h0);
759 gx[1]=3.0*((ya[1]-ya[2])/h2-(ya[2]-ya[1])/h1);
761 det=3.0*(h0+h1)*(h0+h1);
762 this->c[1]=( A*gx[0]-B*gx[1])/det;
763 this->c[2]=(-B*gx[0]+A*gx[1])/det;
764 this->c[0]=this->c[2];
770 for (i=0; i < sys_size-1; i++) {
771 double h_i=xa[i+1]-xa[i];
772 double h_ip1=xa[i+2]-xa[i+1];
773 double ydiff_i=ya[i+1]-ya[i];
774 double ydiff_ip1=ya[i+2]-ya[i+1];
775 double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
776 double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
777 this->offdiag[i]=h_ip1;
778 this->diag[i]=2.0*(h_ip1+h_i);
779 this->g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
784 double h_i=xa[i+1]-xa[i];
785 double h_ip1=xa[1]-xa[0];
786 double ydiff_i=ya[i+1]-ya[i];
787 double ydiff_ip1=ya[1]-ya[0];
788 double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
789 double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
790 this->offdiag[i]=h_ip1;
791 this->diag[i]=2.0*(h_ip1+h_i);
792 this->g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
795 ubvector_range cp1(this->c,range(1,this->c.size()));
798 (this->diag,this->offdiag,this->g,cp1,sys_size,p5m);
799 this->c[0]=this->c[max_index];
806 #ifndef DOXYGEN_INTERNAL 813 (
const interp_cspline_peri<vec_t,vec2_t>&);
835 typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
836 typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
837 typedef boost::numeric::ublas::slice slice;
838 typedef boost::numeric::ublas::range range;
840 #ifndef DOXYGEN_INTERNAL 856 for(
size_t i=0;i<this->
sz-1;i++) {
858 double NE=fabs(umx[3+i]-umx[2+i])+fabs(umx[1+i]-umx[i]);
865 double h_i=(*this->
px)[i+1]-(*this->
px)[i];
866 double NE_next=fabs(umx[4+i]-umx[3+i])+
867 fabs(umx[2+i]-umx[1+i]);
868 double alpha_i=fabs(umx[1+i]-umx[i])/NE;
871 if (NE_next == 0.0) {
874 alpha_ip1=fabs(umx[2+i]-umx[1+i])/NE_next;
875 tL_ip1=(1.0-alpha_ip1)*umx[2+i]+alpha_ip1*umx[3+i];
877 b[i]=(1.0-alpha_i)*umx[1+i]+alpha_i*umx[2+i];
878 c[i]=(3.0*umx[2+i]-2.0*b[i]-tL_ip1)/h_i;
879 d[i]=(b[i]+tL_ip1-2.0*umx[2+i])/(h_i*h_i);
900 virtual void set(
size_t size,
const vec_t &xa,
const vec2_t &ya) {
904 " than "+
szttos(this->min_size)+
" in interp_akima::"+
908 if (size!=this->
sz) {
923 ubvector_range m(um,range(2,um.size()));
925 for (i=0;i<=size-2;i++) {
926 m[i]=(ya[i+1]-ya[i])/(xa[i+1]-xa[i]);
929 um[0]=3.0*m[0]-2.0*m[1];
931 m[this->
sz-1]=2.0*m[size-2]-m[size-3];
932 m[size]=3.0*m[size-2]-2.0*m[size-3];
934 akima_calc(xa,size,um);
940 virtual double eval(
double x0)
const {
943 size_t index=this->
svx.find_const(x0,cache);
945 double x_lo=(*this->
px)[index];
951 return (*this->
py)[index]+delx*(bb+delx*(cc+dd*delx));
955 virtual double deriv(
double x0)
const {
958 size_t index=this->
svx.find_const(x0,cache);
960 double x_lo=(*this->
px)[index];
966 return bb+delx*(2.0*cc+3.0*dd*delx);
975 size_t index=this->
svx.find_const(x0,cache);
977 double x_lo=(*this->
px)[index];
982 return 2.0*cc+6.0*dd*delx;
986 virtual double integ(
double aa,
double bb)
const {
988 size_t i, index_a, index_b;
991 if (((*this->
px)[0]<(*this->
px)[this->
sz-1] && aa>bb) ||
992 ((*this->
px)[0]>(*this->
px)[this->
sz-1] && aa<bb)) {
1000 index_a=this->
svx.find_const(aa,cache);
1001 index_b=this->
svx.find_const(bb,cache);
1005 for(i=index_a; i<=index_b; i++) {
1007 double x_lo=(*this->
px)[i];
1008 double x_hi=(*this->
px)[i+1];
1009 double y_lo=(*this->
py)[i];
1010 double dx=x_hi-x_lo;
1014 if (i==index_a || i==index_b) {
1015 double x1=(i==index_a) ? aa : x_lo;
1016 double x2=(i==index_b) ? bb : x_hi;
1017 result += this->
integ_eval(y_lo,b[i],c[i],d[i],x_lo,x1,x2);
1019 result+=dx*(y_lo+dx*(0.5*b[i]+dx*(c[i]/3.0+0.25*d[i]*dx)));
1023 double y_hi=(*this->
py)[i+1];
1024 std::string str=((std::string)
"Interval of length zero ")+
1027 " in interp_akima::integ().";
1032 if (flip) result*=-1.0;
1038 virtual const char *
type()
const {
return "interp_akima"; }
1040 #ifndef DOXYGEN_INTERNAL 1062 typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
1063 typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
1064 typedef boost::numeric::ublas::slice slice;
1065 typedef boost::numeric::ublas::range range;
1076 virtual const char *
type()
const {
return "interp_akima_peri"; }
1079 virtual void set(
size_t size,
const vec_t &xa,
const vec2_t &ya) {
1083 " than "+
szttos(this->min_size)+
" in interp_akima"+
1087 if (size!=this->
sz) {
1088 this->b.resize(size);
1089 this->c.resize(size);
1090 this->d.resize(size);
1091 this->um.resize(size+4);
1102 ubvector_range m(this->um,range(2,this->um.size()));
1105 for (
size_t i=0;i<=size-2;i++) {
1106 m[i]=(ya[i+1]-ya[i])/(xa[i+1]-xa[i]);
1109 this->um[0]=m[this->
sz-3];
1110 this->um[1]=m[this->
sz-2];
1114 this->akima_calc(xa,size,this->um);
1119 #ifndef DOXYGEN_INTERNAL 1125 (
const interp_akima_peri<vec_t,vec2_t>&);
1139 #ifdef O2SCL_NEVER_DEFINED 1146 typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
1147 typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
1148 typedef boost::numeric::ublas::slice slice;
1149 typedef boost::numeric::ublas::range range;
1151 #ifndef DOXYGEN_INTERNAL 1167 if ((x < 0 && y > 0) || (x > 0 && y < 0)) {
1187 virtual void set(
size_t size,
const vec_t &xa,
const vec2_t &ya) {
1191 " than "+
szttos(this->min_size)+
" in interp_steffen::"+
1195 if (size!=this->
sz) {
1200 y_prime.resize(size);
1214 double h0=(xa[1]-xa[0]);
1215 double s0=(ya[1]-ya[0]) / h0;
1221 for (
size_t i=1; i < (size-1); i++) {
1226 double hi=(xa[i+1]-xa[i]);
1227 double him1=(xa[i]-xa[i-1]);
1230 double si=(ya[i+1]-ya[i]) / hi;
1231 double sim1=(ya[i]-ya[i-1]) / him1;
1234 pi=(sim1*hi + si*him1) / (him1 + hi);
1237 y_prime[i]=(copysign(1.0,sim1)+copysign(1.0,si))*
1238 std::min(fabs(sim1),std::min(fabs(si),0.5*fabs(pi)));
1249 y_prime[size-1]=(ya[size-1]-ya[size-2])/
1250 (xa[size-1]-xa[size-2]);
1253 for (
size_t i=0; i < (size-1); i++) {
1254 double hi=(xa[i+1]-xa[i]);
1255 double si=(ya[i+1]-ya[i]) / hi;
1258 a[i]=(y_prime[i] + y_prime[i+1]-2*si) / hi / hi;
1259 b[i]=(3*si-2*y_prime[i]-y_prime[i+1]) / hi;
1268 virtual double eval(
double x0)
const {
1271 size_t index=this->
svx.find_const(x0,cache);
1272 double x_lo=(*this->
px)[index];
1273 double delx=x0-x_lo;
1276 double y = d[index]+delx*(c[index]+delx*(b[index]+delx*a[index]));
1285 size_t index=this->
svx.find_const(x0,cache);
1286 double x_lo=(*this->
px)[index];
1287 double delx=x0-x_lo;
1289 return c[index]+delx*(2.0*b[index]+delx*3.0*a[index]);
1298 size_t index=this->
svx.find_const(x0,cache);
1299 double x_lo=(*this->
px)[index];
1300 double delx=x0-x_lo;
1302 return 2.0*b[index]+delx*6.0*a[index];
1306 virtual double integ(
double al,
double bl)
const {
1308 size_t i, index_a, index_b;
1311 if (((*this->
px)[0]<(*this->
px)[this->
sz-1] && al>bl) ||
1312 ((*this->
px)[0]>(*this->
px)[this->
sz-1] && al<bl)) {
1320 index_a=this->
svx.find_const(al,cache);
1321 index_b=this->
svx.find_const(bl,cache);
1325 for(i=index_a; i<=index_b; i++) {
1327 double x_lo=(*this->
px)[i];
1328 double x_hi=(*this->
px)[i+1];
1329 double y_lo=(*this->
py)[i];
1330 double y_hi=(*this->
py)[i+1];
1331 double dx=x_hi-x_lo;
1335 double x1=(i == index_a) ? al-x_lo : 0.0;
1336 double x2=(i == index_b) ? bl-x_lo : x_hi-x_lo;
1337 result += (1.0/4.0)*a[i]*(x2*x2*x2*x2-x1*x1*x1*
x1)+
1338 (1.0/3.0)*b[i]*(x2*x2*x2-x1*x1*
x1)+
1339 (1.0/2.0)*c[i]*(x2*x2-x1*
x1)+d[i]*(x2-x1);
1342 std::string str=((std::string)
"Interval of length zero ")+
1345 " in interp_steffen::integ().";
1351 if (flip) result*=-1.0;
1357 virtual const char *
type()
const {
return "interp_steffen"; }
1359 #ifndef DOXYGEN_INTERNAL 1365 (
const interp_steffen<vec_t,vec2_t>&);
1398 #ifdef O2SCL_NEVER_DEFINED 1427 virtual void set(
size_t size,
const vec_t &x,
const vec2_t &y) {
1432 " than "+
szttos(this->min_size)+
" in interp_monotonic::"+
1440 if (this->
sz!=size) {
1442 Delta.resize(size-1);
1443 alpha.resize(size-1);
1444 beta.resize(size-1);
1455 for(
size_t i=0;i<size-1;i++) {
1456 Delta[i]=(y[i+1]-y[i])/(x[i+1]-x[i]);
1458 m[i]=(Delta[i]+Delta[i-1])/2.0;
1462 m[size-1]=Delta[size-2];
1465 for(
size_t i=0;i<size-1;i++) {
1473 for(
size_t i=0;i<size-1;i++) {
1474 alpha[i]=m[i]/Delta[i];
1475 beta[i]=m[i+1]/Delta[i];
1479 for(
size_t i=0;i<size-1;i++) {
1480 double norm2=alpha[i]*alpha[i]+beta[i]*beta[i];
1482 double tau=3.0/sqrt(norm2);
1483 m[i]=tau*alpha[i]*Delta[i];
1484 m[i+1]=tau*beta[i]*Delta[i];
1492 virtual double eval(
double x0)
const {
1495 size_t index=this->
svx.find_const(x0,cache);
1497 double x_lo=(*this->
px)[index];
1498 double x_hi=(*this->
px)[index+1];
1499 double y_lo=(*this->
py)[index];
1500 double y_hi=(*this->
py)[index+1];
1502 double t=(x0-x_lo)/h;
1503 double t2=t*t, t3=t2*t;
1505 double h00=2.0*t3-3.0*t2+1.0;
1506 double h10=t3-2.0*t2+t;
1507 double h01=-2.0*t3+3.0*t2;
1509 double interp=y_lo*h00+h*m[index]*h10+y_hi*h01+h*m[index+1]*h11;
1518 size_t index=this->
svx.find_const(x0,cache);
1520 double x_lo=(*this->
px)[index];
1521 double x_hi=(*this->
px)[index+1];
1522 double y_lo=(*this->
py)[index];
1523 double y_hi=(*this->
py)[index+1];
1525 double t=(x0-x_lo)/h;
1528 double dh00=6.0*t2-6.0*t;
1529 double dh10=3.0*t2-4.0*t+1.0;
1530 double dh01=-6.0*t2+6.0*t;
1531 double dh11=3.0*t2-2.0*t;
1532 double deriv=(y_lo*dh00+h*m[index]*dh10+y_hi*dh01+
1533 h*m[index+1]*dh11)/h;
1544 size_t index=this->
svx.find_const(x0,cache);
1546 double x_lo=(*this->
px)[index];
1547 double x_hi=(*this->
px)[index+1];
1548 double y_lo=(*this->
py)[index];
1549 double y_hi=(*this->
py)[index+1];
1551 double t=(x0-x_lo)/h;
1553 double ddh00=12.0*t-6.0;
1554 double ddh10=6.0*t-4.0;
1555 double ddh01=-12.0*t+6.0;
1556 double ddh11=6.0*t-2.0;
1557 double deriv2=(y_lo*ddh00+h*m[index]*ddh10+y_hi*ddh01+
1558 h*m[index+1]*ddh11)/h/h;
1564 virtual double integ(
double a,
double b)
const {
1566 size_t i, index_a, index_b;
1569 if (((*this->
px)[0]<(*this->
px)[this->
sz-1] && a>b) ||
1570 ((*this->
px)[0]>(*this->
px)[this->
sz-1] && a<b)) {
1578 index_a=this->
svx.find_const(a,cache);
1579 index_b=this->
svx.find_const(b,cache);
1583 for(i=index_a; i<=index_b; i++) {
1585 double x_hi=(*this->
px)[i+1];
1586 double x_lo=(*this->
px)[i];
1587 double y_lo=(*this->
py)[i];
1588 double y_hi=(*this->
py)[i+1];
1600 double t=(x_hi-x_lo)/h;
1601 double t2=t*t, t3=t2*t, t4=t3*t;
1603 double ih00=t4/2.0-t3+t;
1604 double ih10=t4/4.0-2.0*t3/3.0+t2/2.0;
1605 double ih01=-t4/2.0+t3;
1606 double ih11=t4/4.0-t3/3.0;
1607 double intres=h*(y_lo*ih00+h*m[i]*ih10+y_hi*ih01+
1612 std::string str=((std::string)
"Interval of length zero ")+
1615 " in interp_monotonic::integ().";
1621 if (flip) result*=-1.0;
1627 virtual const char *
type()
const {
return "interp_monotonic"; }
1629 #ifndef DOXYGEN_INTERNAL 1635 (
const interp_monotonic<vec_t,vec2_t>&);
1653 template<
class vec_t=boost::numeric::ublas::vector<
double>,
1656 #ifndef DOXYGEN_INTERNAL 1686 O2SCL_ERR((((std::string)
"Invalid interpolation type, ")+
1699 virtual double eval(
const double x0,
size_t n,
const vec_t &x,
1702 return itp->
eval(x0);
1710 virtual double deriv(
const double x0,
size_t n,
const vec_t &x,
1713 return itp->
deriv(x0);
1720 virtual double deriv2(
const double x0,
size_t n,
const vec_t &x,
1730 virtual double integ(
const double x1,
const double x2,
size_t n,
1731 const vec_t &x,
const vec2_t &y) {
1733 return itp->
integ(x1,x2);
1756 O2SCL_ERR((((std::string)
"Invalid interpolation type, ")+
1763 #ifndef DOXYGEN_INTERNAL 1794 template<
class vec_t=boost::numeric::ublas::vector<
double>,
1798 #ifndef DOXYGEN_INTERNAL 1823 const vec2_t &y,
size_t interp_type=
itp_cspline) {
1826 O2SCL_ERR((((std::string)
"Vector endpoints equal (value=")+
1848 O2SCL_ERR((((std::string)
"Invalid interpolation type, ")+
1850 "interp_vec::interp_vec().").c_str(),
exc_einval);
1866 void set(
size_t n,
const vec_t &x,
const vec2_t &y) {
1873 void set(
size_t n,
const vec_t &x,
1874 const vec2_t &y,
size_t interp_type) {
1877 O2SCL_ERR((((std::string)
"Vector endpoints equal (value=")+
1900 O2SCL_ERR((((std::string)
"Invalid interpolation type, ")+
1920 virtual double eval(
const double x0)
const {
1922 O2SCL_ERR(
"No vector set in interp_vec::eval().",
1925 return itp->
eval(x0);
1931 O2SCL_ERR(
"No vector set in interp_vec::operator().",
1934 return itp->
eval(x0);
1938 virtual double deriv(
const double x0)
const {
1940 O2SCL_ERR(
"No vector set in interp_vec::deriv().",
1943 return itp->
deriv(x0);
1949 virtual double deriv2(
const double x0)
const {
1951 O2SCL_ERR(
"No vector set in interp_vec::deriv2().",
1958 virtual double integ(
const double x1,
const double x2)
const {
1960 O2SCL_ERR(
"No vector set in interp_vec::integ().",
1963 return itp->
integ(x1,x2);
1968 return "interp_vec";
1971 #ifndef DOXYGEN_INTERNAL 1977 (
const interp_vec<vec_t,vec2_t> &it);
1988 public interp<double[n]> {
1994 :
interp<double[n]>(interp_type) {}
2014 size_t interp_type) :
2040 (
double level,
size_t n, vec_t &x, vec2_t &y) {
2043 O2SCL_ERR2(
"Need at least two data points in ",
2050 if (y[0]==level) count++;
2053 for(
size_t i=0;i<n-1;i++) {
2055 if ((y[i]<level && y[i+1]>=level) ||
2056 (y[i]>level && y[i+1]<=level)) {
2085 (
double level,
size_t n, vec_t &x, vec2_t &y, std::vector<double> &locs) {
2088 O2SCL_ERR2(
"Need at least two data points in ",
2097 locs.push_back(x[0]);
2101 for(
size_t i=0;i<n-1;i++) {
2103 if ((y[i]<level && y[i+1]>level) ||
2104 (y[i]>level && y[i+1]<level)) {
2107 double x0=x[i]+(x[i+1]-x[i])*(level-y[i])/(y[i+1]-y[i]);
2109 }
else if (y[i+1]==level) {
2110 locs.push_back(x[i+1]);
2120 template<
class ovec_t,
class vec2_t>
2124 for(
size_t i=0;i<n;i++) grid[i]=((
double)i);
2126 for(
size_t i=0;i<n;i++) dv[i]=oi.
deriv(((
double)i));
2133 template<
class ovec_t,
class vec2_t>
2137 for(
size_t i=0;i<n;i++) grid[i]=((
double)i);
2139 for(
size_t i=0;i<n;i++) dv[i]=oi.
deriv2(((
double)i));
2146 template<
class vec_t,
class vec2_t,
class vec3_t>
2150 for(
size_t i=0;i<n;i++) dv[i]=oi.
deriv(vx[i]);
2157 template<
class vec_t,
class vec2_t,
class vec3_t>
2161 for(
size_t i=0;i<n;i++) dv[i]=oi.
deriv(vx[i]);
2167 template<
class ovec_t>
2170 for(
size_t i=0;i<n;i++) grid[i]=((
double)i);
2172 return oi.
integ(0.0,((
double)(n-1)));
2188 template<
class vec_t,
class vec2_t>
2196 double total=si.
integ(x[0],x[n-1],n,x,y);
2204 template<
class vec_t,
class vec2_t,
class vec3_t>
2213 for(
size_t i=1;i<n;i++) {
2214 iy[i]=si.
integ(x[0],x[i],n,x,y);
2223 template<
class ovec_t>
2225 ovec_t &v,
size_t interp_type) {
2227 for(
size_t i=0;i<n;i++) grid[i]=((
double)i);
2229 return oi.
integ(0.0,x2);
2235 template<
class vec_t,
class vec2_t>
2237 const vec2_t &y,
double x2,
2244 double total=si.
integ(x[0],x2,n,x,y);
2287 (
double sum,
size_t n, vec_t &x, vec2_t &y,
double &lev,
2288 int boundaries=0,
int verbose=0,
bool err_on_fail=
true) {
2291 O2SCL_ERR2(
"Need at least two data points in ",
2298 std::vector<double>
x2, y2;
2300 if (boundaries==1) {
2302 std::cout <<
"Fix left boundary to zero." << std::endl;
2306 x2[0]=x[0]-(x[1]-x[0])/1.0e6;
2308 for(
size_t i=0;i<n;i++) {
2313 }
else if (boundaries==2) {
2315 std::cout <<
"Fix right boundary to zero." << std::endl;
2319 for(
size_t i=0;i<n;i++) {
2323 x2[n]=x[n-1]+(x[n-1]-x[n-2])/1.0e6;
2326 }
else if (boundaries==3) {
2328 std::cout <<
"Fix both boundaries to zero." << std::endl;
2332 x2[0]=x[0]-(x[1]-x[0])/1.0e6;
2334 for(
size_t i=0;i<n;i++) {
2338 x2[n+1]=x[n-1]+(x[n-1]-x[n-2])/1.0e6;
2343 std::cout <<
"No boundary extrapolation." << std::endl;
2347 for(
size_t i=0;i<n;i++) {
2357 vector_sort<ubvector,double>(ysort.size(),ysort);
2363 std::vector<double> ylist;
2364 for(
size_t i=0;i<ysort.size()-1;i++) {
2365 if (ysort[i]!=ysort[i+1]) {
2366 ylist.push_back((ysort[i+1]+3.0*ysort[i])/4.0);
2367 ylist.push_back((ysort[i+1]*3.0+ysort[i])/4.0);
2377 std::vector<double> xi, yi;
2381 for(
size_t k=0;k<ylist.size();k++) {
2382 double lev_tmp=ylist[k];
2383 std::vector<double> locs;
2385 if ((locs.size()%2)!=0) {
2388 std::cout << k <<
" " << lev_tmp <<
" " << 0.0 <<
" " 2389 << locs.size() <<
" (fail)" << std::endl;
2392 double sum_temp=0.0;
2393 for(
size_t i=0;i<locs.size()/2;i++) {
2394 double x0=locs[2*i];
2395 double x1=locs[2*i+1];
2396 sum_temp+=itp2.
integ(x0,x1,n2,x2,y2);
2398 xi.push_back(sum_temp);
2399 yi.push_back(lev_tmp);
2401 std::cout << k <<
" " << lev_tmp <<
" " << sum_temp <<
" " 2402 << locs.size() <<
" ";
2403 for(
size_t i=0;i<locs.size();i++) {
2404 std::cout << locs[i] <<
" ";
2406 std::cout << std::endl;
2418 lev=itp2.
eval(sum,xi.size(),xi,yi);
2426 (
size_t n, vec_t &x, vec2_t &y,
double intl, std::vector<double> &locs,
2427 int boundaries=0,
int verbose=0,
bool err_on_fail=
true) {
2432 std::cout <<
"Total integral: " << total << std::endl;
2436 std::cout <<
"Target integral: " << intl << std::endl;
2441 boundaries,verbose,err_on_fail);
2444 O2SCL_ERR2(
"Failed to find a level which enclosed the ",
2445 "specified integral in vector_region_int().",
2451 std::cout <<
"Level from vector_invert: " << lev << std::endl;
2457 std::cout <<
"Locations from vector_find_level: ";
2458 for(
size_t i=0;i<locs.size();i++) {
2459 std::cout << locs[i];
2460 if (i!=locs.size()-1) {
2464 std::cout << std::endl;
2472 (
size_t n, vec_t &x, vec2_t &y,
double frac, std::vector<double> &locs,
2473 int boundaries=0,
int verbose=0,
bool err_on_fail=
true) {
2476 verbose,err_on_fail);
2486 (
size_t n, vec_t &x, vec2_t &y,
double frac,
double &low,
double &high,
2487 int boundaries=0,
int verbose=0,
bool err_on_fail=
true) {
2489 std::vector<double> locs;
2491 verbose,err_on_fail);
2492 if (locs.size()==0 || ret!=0) {
2494 O2SCL_ERR2(
"Zero level crossings or vector_region_fracint() ",
2513 (
size_t n, vec_t &x, vec2_t &y,
double frac,
double &low,
double &high,
2514 int boundaries=0,
int verbose=0,
bool err_on_fail=
true) {
2516 std::vector<double> locs;
2518 verbose,err_on_fail);
2519 if (locs.size()==0 || ret!=0) {
2521 O2SCL_ERR2(
"Zero level crossings or vector_region_int() ",
2536 template<
class vec_t,
class vec2_t,
class vec3_t,
class vec4_t>
2538 vec3_t &new_x, vec4_t &new_y,
size_t n_pts,
2539 size_t interp_type) {
2541 if (x.size()!=y.size()) {
2542 O2SCL_ERR2(
"The x and y vectors must have the same size ",
2546 O2SCL_ERR2(
"Number of points must be at least 2 ",
2558 new_y.resize(n_pts);
2560 for(
size_t i=0;i<n_pts;i++) {
2561 new_y[i]=itp.
eval(new_x[i]);
2571 template<
class vec_t,
class vec2_t,
class vec3_t,
class vec4_t>
2573 vec3_t &new_x, vec4_t &new_y,
size_t n_pts,
2574 size_t interp_type1,
size_t interp_type2,
2575 double acc=1.0e-4) {
2577 if (x.size()!=y.size()) {
2578 O2SCL_ERR2(
"The x and y vectors must have the same size ",
2582 O2SCL_ERR2(
"Number of points must be at least 2 ",
2594 std::vector<double> new_y2(n_pts);
2595 new_y.resize(n_pts);
2598 for(
size_t i=0;i<n_pts;i++) {
2599 new_y[i]=itp1.
eval(new_x[i]);
2600 new_y2[i]=itp2.
eval(new_x[i]);
2603 double max_y, min_y;
2605 for(
size_t i=0;i<n_pts;i++) {
2606 if (fabs(new_y2[i]-new_y[i])/(max_y-min_y)>acc) {
2618 template<
class vec_t,
class vec2_t>
2621 static const size_t n2=11;
2624 std::vector<double> xn, yn;
2628 double min_y, max_y;
2630 for(
size_t i=0;i<n2;i++) {
2631 yn[i]=(yn[i]-min_y)/(max_y-min_y);
2637 for(
size_t i=0;i<n2;i++) {
2638 double ty=((double)i)/10.0;
2639 chi2+=pow(yn[i]-ty,2.0);
2654 template<
class vec_t,
class vec2_t>
2656 if (x.size()!=y.size()) {
2657 O2SCL_ERR2(
"The x and y vectors must have the same size ",
2664 O2SCL_ERR2(
"Vector size must be at least 2 ",
2673 bool x_positive=
true;
2674 bool y_positive=
true;
2675 for(
size_t i=0;i<n;i++) {
2676 if (x[i]<=0.0) x_positive=
false;
2677 if (y[i]<=0.0) y_positive=
false;
2680 if (x_positive==
false && y_positive==
false)
return;
2685 std::vector<double> lx(n), ly(n);
2687 for(
size_t i=0;i<n;i++) {
2692 for(
size_t i=0;i<n;i++) {
2703 if (chi2_xy<chi2_x && chi2_xy<chi2_y && chi2_xy<chi2) {
2706 }
else if (chi2_x<chi2_xy && chi2_x<chi2_y && chi2_x<chi2) {
2708 }
else if (chi2_y<chi2_xy && chi2_y<chi2_x && chi2_y<chi2) {
2713 if (chi2_x<chi2) log_x=
true;
2717 if (chi2_y<chi2) log_y=
true;
2732 template<
class vec_t>
2734 std::vector<double> x(y.size());
2735 for(
size_t i=0;i<y.size();i++) {
2743 #ifndef DOXYGEN_NO_O2NS double linear_or_log_chi2(const vec_t &x, const vec2_t &y)
Rebin, rescale, sort, and match to .
virtual const char * type() const
Return the type, "interp_steffen".
virtual const char * type() const
Return the type, "interp_cspline_peri".
Interpolation class for pre-specified vector.
o2scl_linalg::ubvector_5_mem p5m
Memory for the tridiagonalization.
void vector_deriv2_xy_interp(size_t n, vec_t &vx, vec2_t &vy, vec3_t &dv, size_t interp_type=itp_linear)
Compute second derivative at all points from an interpolation object.
virtual double eval(double x0) const
Give the value of the function .
virtual void set(size_t size, const vec_t &x, const vec2_t &y)=0
Initialize interpolation routine.
virtual double eval(const double x0, size_t n, const vec_t &x, const vec2_t &y)
Give the value of the function, , as specified as the first n elements of vectors x and y...
interp_vec(size_t n, const vec_t &x, const vec2_t &y, size_t interp_type=itp_cspline)
Create an interpolation object with interpolation type itp_cspline based on the first n entries of ve...
virtual const char * type() const
Return the type, "interp_linear".
double integ_eval(double ai, double bi, double ci, double di, double xi, double a, double b) const
An internal function to assist in computing the integral for both the cspline and Akima types...
virtual double operator()(double x0) const
Give the value of the function .
virtual const char * type() const =0
Return the type.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
interp(size_t interp_type=itp_cspline)
Create with base interpolation object it.
void solve_tridiag_sym(const vec_t &diag, const vec2_t &offdiag, const vec3_t &b, vec4_t &x, size_t N, mem_t &m)
Solve a symmetric tridiagonal linear system with user-specified memory.
void vector_deriv_interp(size_t n, ovec_t &v, vec2_t &dv, size_t interp_type=itp_linear)
Compute derivative at all points from an interpolation object.
void akima_calc(const vec_t &x_array, size_t size, ubvector &umx)
For initializing the interpolation.
int vector_invert_enclosed_sum(double sum, size_t n, vec_t &x, vec2_t &y, double &lev, int boundaries=0, int verbose=0, bool err_on_fail=true)
Compute the endpoints which enclose the regions whose integral is equal to sum.
size_t min_size
The minimum size of the vectors to interpolate between.
A specialization of interp for C-style double arrays.
virtual double deriv(double x0) const
Give the value of the derivative .
interp_base< vec_t, vec2_t > * itp
Pointer to base interpolation object.
virtual double integ(double a, double b) const =0
Give the value of the integral .
virtual double deriv2(const double x0) const
Give the value of the second derivative .
void coeff_calc(const ubvector &c_array, double dy, double dx, size_t index, double &b, double &c2, double &d) const
Compute coefficients for cubic spline interpolation.
Akima spline interpolation with periodic boundary conditions (GSL)
virtual double integ(double al, double bl) const
Give the value of the integral .
void clear()
Manually clear the pointer to the user-specified vector.
virtual const char * type() const
Return the type, "interp_cspline".
virtual const char * type() const
Return the type, "interp_akima".
double vector_integ_ul_xy_interp(size_t n, const vec_t &x, const vec2_t &y, double x2, size_t interp_type=itp_linear)
Compute the integral over y(x) using interpolation up to a specified upper limit. ...
int vector_bound_fracint(size_t n, vec_t &x, vec2_t &y, double frac, double &low, double &high, int boundaries=0, int verbose=0, bool err_on_fail=true)
Find the boundaries of the region enclosing a integral.
virtual const char * type() const
Return the type, "interp_monotonic".
double copysign(const double x, const double y)
Flip the sign of x if x and y have different signs.
virtual double eval(double x0) const
Give the value of the function .
invalid argument supplied by user
const vec_t * px
Independent vector.
virtual double eval(double x0) const
Give the value of the function .
virtual double integ(double a, double b) const
Give the value of the integral .
virtual double deriv(double x0) const
Give the value of the derivative .
Monotonicity-preserving interpolation.
Base low-level interpolation class [abstract base].
interp_array_vec(size_t nv, const arr_t &x, const arr_t &y, size_t interp_type)
Create with base interpolation object it.
virtual const char * type() const
Return the type, "interp_akima_peri".
double vector_integ_xy_interp(size_t n, const vec_t &x, const vec2_t &y, size_t interp_type=itp_linear)
Compute the integral over y(x) using interpolation.
virtual double deriv(double x0) const
Give the value of the derivative .
virtual double deriv(const double x0) const
Give the value of the derivative .
Akima spline for periodic boundary conditions.
void vector_min(size_t n, const vec_t &data, size_t &index, data_t &val)
Compute the minimum of the first n elements of a vector.
virtual double integ(double a, double b) const
Give the value of the integral .
virtual double deriv2(double x0) const
Give the value of the second derivative .
interp_vec()
Blank interpolator.
void vector_deriv2_interp(size_t n, ovec_t &v, vec2_t &dv, size_t interp_type=itp_linear)
Compute second derivative at all points from an interpolation object.
virtual double deriv(double x0) const =0
Give the value of the derivative .
virtual double integ(double aa, double bb) const
Give the value of the integral .
virtual double deriv2(double x0) const
Give the value of the second derivative .
Steffen's monotonic method.
virtual double deriv2(double x0) const
Give the value of the second derivative .
virtual double integ(double a, double b) const
Give the value of the integral .
Steffen's monotonicity-preserving interpolation.
interp_akima()
Create a base interpolation object with or without periodic boundary conditions.
Cubic spline interpolation with periodic boundary conditions (GSL)
virtual const char * type() const
Return the type, "interp_nearest_neigh".
virtual double operator()(double x0) const
Give the value of the function .
virtual double deriv(double x0) const
Give the value of the derivative .
virtual double integ(const double x1, const double x2, size_t n, const vec_t &x, const vec2_t &y)
Give the value of the integral , where is specified in the first n elements of vectors x and y...
virtual double integ(double a, double b) const
Give the value of the integral .
Cubic spline for natural boundary conditions.
double vector_integ_interp(size_t n, ovec_t &v, size_t interp_type)
Integral of a vector from interpolation object.
virtual double deriv(double x0) const
Give the value of the derivative .
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
ubvector beta
Staggered ratio.
Allocation object for 5 arrays of equal size.
virtual const char * type() const
Return the type, "interp_vec".
interp_array()
Create an interpolator using the default base interpolation objects.
void linear_or_log(vec_t &x, vec2_t &y, bool &log_x, bool &log_y)
Attempt to determine if data represented by (x,y) would be better plotted on a semi-log or log-log pl...
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
int vector_region_fracint(size_t n, vec_t &x, vec2_t &y, double frac, std::vector< double > &locs, int boundaries=0, int verbose=0, bool err_on_fail=true)
Find the region enclosing a partial integral.
size_t itype
Interpolation type.
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
virtual double integ(const double x1, const double x2) const
Give the value of the integral .
void vector_minmax_value(size_t n, vec_t &data, data_t &min, data_t &max)
Compute the minimum and maximum of the first n elements of a vector.
virtual double eval(double x0) const
Give the value of the function .
const vec2_t * py
Dependent vector.
virtual double eval(double x0) const
Give the value of the function .
o2scl_linalg::ubvector_4_mem p4m
Memory for the tridiagonalization.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
void vector_deriv_xy_interp(size_t n, vec_t &vx, vec2_t &vy, vec3_t &dv, size_t interp_type=itp_linear)
Compute derivative at all points from an interpolation object.
Allocation object for 4 arrays of equal size.
void set_type(size_t interp_type)
Set base interpolation type.
size_t vector_level_count(double level, size_t n, vec_t &x, vec2_t &y)
Count level crossings.
Akima spline for natural boundary conditions.
int vector_bound_int(size_t n, vec_t &x, vec2_t &y, double frac, double &low, double &high, int boundaries=0, int verbose=0, bool err_on_fail=true)
Find the boundaries of the region enclosing a integral.
Akima spline interpolation (GSL)
size_t ordered_lookup(const double x0) const
Find the index of x0 in the ordered array x.
interp_steffen()
Create a base interpolation object.
virtual double eval(double x0) const =0
Give the value of the function .
virtual double deriv2(double x0) const
Give the value of the second derivative .
void vector_max(size_t n, const vec_t &data, size_t &index, data_t &val)
Compute the maximum of the first n elements of a vector.
int vector_region_int(size_t n, vec_t &x, vec2_t &y, double intl, std::vector< double > &locs, int boundaries=0, int verbose=0, bool err_on_fail=true)
Find the region enclosing an integral.
interp_cspline()
Create a base interpolation object with natural or periodic boundary conditions.
void rebin_xy(const vec_t &x, const vec2_t &y, vec3_t &new_x, vec4_t &new_y, size_t n_pts, size_t interp_type)
From an (x,y) pair, create a new (x,y) pair using interpolation where the new x vector is uniformly s...
virtual double deriv(double x0) const
Give the value of the derivative .
Cubic spline for periodic boundary conditions.
Monotonicity-preserving interpolation (unfinished)
virtual double deriv2(double x0) const
Give the value of the second derivative (always zero)
virtual double deriv2(const double x0, size_t n, const vec_t &x, const vec2_t &y)
Give the value of the second derivative, , where is specified in the first n elements of vectors x ...
ubvector Delta
Finite differences.
Cubic spline interpolation (GSL)
interp_array(size_t interp_type)
Create with base interpolation objects it and rit.
Nearest-neighbor interpolation.
virtual double deriv2(double x0) const =0
Give the value of the second derivative .
void solve_cyc_tridiag_sym(const vec_t &diag, const vec2_t &offdiag, const vec3_t &b, vec4_t &x, size_t N, mem_t &m)
Solve a symmetric cyclic tridiagonal linear system with user specified memory.
void set_vec(size_t nn, const vec_t &x)
Set the vector to be searched.
double vector_integ_ul_interp(size_t n, double x2, ovec_t &v, size_t interp_type)
Compute the integral of a vector using interpolation up to a specified upper limit.
search_vec< const vec_t > svx
To perform binary searches.
static const double x2[5]
virtual double eval(const double x0) const
Give the value of the function .
static const double x1[5]
Linear interpolation (GSL)
std::string szttos(size_t x)
Convert a size_t to a string.
virtual double deriv2(double x0) const
Give the value of the second derivative (always zero)
void vector_find_level(double level, size_t n, vec_t &x, vec2_t &y, std::vector< double > &locs)
Perform inverse linear interpolation.
virtual double eval(double x0) const
Give the value of the function .
A specialization of o2scl::interp_vec for C-style arrays.
interp_base< vec_t, vec2_t > * itp
Base interpolation object.
void vector_sort_double(size_t n, vec_t &data)
Sort a vector of doubles (in increasing order)
Interpolation class for general vectors.
virtual double deriv(const double x0, size_t n, const vec_t &x, const vec2_t &y)
Give the value of the derivative, , where is specified in the first n elements of vectors x and y...