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> 62 #ifndef DOXYGEN_NO_O2NS 106 #ifdef O2SCL_NEVER_DEFINED 110 #ifndef DOXYGEN_INTERNAL 133 double integ_eval(
double ai,
double bi,
double ci,
double di,
double xi,
134 double a,
double b)
const {
139 double bterm=0.5*bi*r12;
140 double cterm=ci*(r1*r1+r2*r2+r1*r2)/3.0;
141 double dterm=0.25*di*r12*(r1*r1+r2*r2);
143 return (b-a)*(ai+bterm+cterm+dterm);
165 virtual void set(
size_t size,
const vec_t &x,
const vec2_t &y)=0;
168 virtual double eval(
double x0)
const=0;
176 virtual double deriv(
double x0)
const=0;
181 virtual double deriv2(
double x0)
const=0;
184 virtual double integ(
double a,
double b)
const=0;
187 virtual const char *
type()
const=0;
189 #ifndef DOXYGEN_INTERNAL 210 #ifdef O2SCL_NEVER_DEFINED 223 virtual void set(
size_t size,
const vec_t &x,
const vec2_t &y) {
226 " than "+
szttos(this->min_size)+
" in interp_linear::"+
237 virtual double eval(
double x0)
const {
240 size_t index=this->
svx.find_const(x0,cache);
242 double x_lo=(*this->
px)[index];
243 double x_hi=(*this->
px)[index+1];
244 double y_lo=(*this->
py)[index];
245 double y_hi=(*this->
py)[index+1];
248 return y_lo+(x0-x_lo)/dx*(y_hi-y_lo);
252 virtual double deriv(
double x0)
const {
255 size_t index=this->
svx.find_const(x0,cache);
257 double x_lo=(*this->
px)[index];
258 double x_hi=(*this->
px)[index+1];
259 double y_lo=(*this->
py)[index];
260 double y_hi=(*this->
py)[index+1];
275 virtual double integ(
double a,
double b)
const {
278 size_t i, index_a, index_b;
281 if (((*this->
px)[0]<(*this->
px)[this->
sz-1] && a>b) ||
282 ((*this->
px)[0]>(*this->
px)[this->
sz-1] && a<b)) {
289 index_a=this->
svx.find_const(a,cache);
290 index_b=this->
svx.find_const(b,cache);
293 for(i=index_a; i<=index_b; i++) {
295 double x_lo=(*this->
px)[i];
296 double x_hi=(*this->
px)[i+1];
297 double y_lo=(*this->
py)[i];
298 double y_hi=(*this->
py)[i+1];
303 if (i == index_a || i == index_b) {
304 double x1=(i == index_a) ? a : x_lo;
305 double x2=(i == index_b) ? b : x_hi;
306 double D=(y_hi-y_lo)/dx;
307 result += (x2-
x1)*(y_lo+0.5*D*((x2-x_lo)+(x1-x_lo)));
309 result += 0.5*dx*(y_lo+y_hi);
313 std::string str=((std::string)
"Interval of length zero ")+
316 " in interp_linear::integ().";
321 if (flip) result=-result;
326 virtual const char *
type()
const {
return "interp_linear"; }
328 #ifndef DOXYGEN_INTERNAL 351 #ifdef O2SCL_NEVER_DEFINED 358 typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
359 typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
360 typedef boost::numeric::ublas::slice slice;
361 typedef boost::numeric::ublas::range range;
363 #ifndef DOXYGEN_INTERNAL 379 void coeff_calc(
const ubvector &c_array,
double dy,
double dx,
380 size_t index,
double &b,
double &c2,
double &d)
const {
382 double c_i=c_array[index];
383 double c_ip1=c_array[index+1];
384 b=(dy/dx)-dx*(c_ip1+2.0*c_i)/3.0;
386 d=(c_ip1-c_i)/(3.0*dx);
407 virtual void set(
size_t size,
const vec_t &xa,
const vec2_t &ya) {
411 " than "+
szttos(this->min_size)+
" in interp_cspline::"+
415 if (size!=this->
sz) {
419 offdiag.resize(size);
432 size_t num_points=size;
433 size_t max_index=num_points-1;
434 size_t sys_size=max_index-1;
439 for (i=0; i < sys_size; i++) {
440 double h_i=xa[i+1]-xa[i];
441 double h_ip1=xa[i+2]-xa[i+1];
442 double ydiff_i=ya[i+1]-ya[i];
443 double ydiff_ip1=ya[i+2]-ya[i+1];
444 double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
445 double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
447 diag[i]=2.0*(h_ip1+h_i);
448 g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
458 ubvector_range cp1(c,range(1,c.size()));
461 (diag,offdiag,g,cp1,sys_size,p4m);
467 virtual double eval(
double x0)
const {
470 size_t index=this->
svx.find_const(x0,cache);
472 double x_lo=(*this->
px)[index];
473 double x_hi=(*this->
px)[index+1];
476 double y_lo=(*this->
py)[index];
477 double y_hi=(*this->
py)[index+1];
480 double b_i, c_i, d_i;
482 coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
484 return y_lo+delx*(b_i+delx*(c_i+delx*d_i));
488 virtual double deriv(
double x0)
const {
491 size_t index=this->
svx.find_const(x0,cache);
493 double x_lo=(*this->
px)[index];
494 double x_hi=(*this->
px)[index+1];
497 double y_lo=(*this->
py)[index];
498 double y_hi=(*this->
py)[index+1];
501 double b_i, c_i, d_i;
503 coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
505 return b_i+delx*(2.0*c_i+3.0*d_i*delx);
514 size_t index=this->
svx.find_const(x0,cache);
516 double x_lo=(*this->
px)[index];
517 double x_hi=(*this->
px)[index+1];
520 double y_lo=(*this->
py)[index];
521 double y_hi=(*this->
py)[index+1];
524 double b_i, c_i, d_i;
526 coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
528 return 2.0*c_i+6.0*d_i*delx;
532 virtual double integ(
double a,
double b)
const {
534 size_t i, index_a, index_b;
537 if (((*this->
px)[0]<(*this->
px)[this->
sz-1] && a>b) ||
538 ((*this->
px)[0]>(*this->
px)[this->
sz-1] && a<b)) {
546 index_a=this->
svx.find_const(a,cache);
547 index_b=this->
svx.find_const(b,cache);
551 for(i=index_a; i<=index_b; i++) {
553 double x_lo=(*this->
px)[i];
554 double x_hi=(*this->
px)[i+1];
555 double y_lo=(*this->
py)[i];
556 double y_hi=(*this->
py)[i+1];
561 double b_i, c_i, d_i;
562 coeff_calc(c,dy,dx,i,b_i,c_i,d_i);
563 if (i == index_a || i == index_b) {
564 double x1=(i == index_a) ? a : x_lo;
565 double x2=(i == index_b) ? b : x_hi;
566 result += this->
integ_eval(y_lo,b_i,c_i,d_i,x_lo,x1,x2);
568 result += dx*(y_lo+dx*(0.5*b_i+
569 dx*(c_i/3.0+0.25*d_i*dx)));
572 std::string str=((std::string)
"Interval of length zero ")+
575 " in interp_cspline::integ().";
581 if (flip) result*=-1.0;
587 virtual const char *
type()
const {
return "interp_cspline"; }
589 #ifndef DOXYGEN_INTERNAL 595 (
const interp_cspline<vec_t,vec2_t>&);
609 #ifdef O2SCL_NEVER_DEFINED 621 typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
622 typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
623 typedef boost::numeric::ublas::slice slice;
624 typedef boost::numeric::ublas::range range;
634 virtual const char *
type()
const {
return "interp_cspline_peri"; }
638 virtual void set(
size_t size,
const vec_t &xa,
const vec2_t &ya) {
642 " than "+
szttos(this->min_size)+
" in interp_cspline"+
646 if (size!=this->
sz) {
647 this->c.resize(size);
648 this->g.resize(size);
649 this->diag.resize(size);
650 this->offdiag.resize(size);
663 size_t num_points=size;
665 size_t max_index=num_points-1;
667 size_t sys_size=max_index;
673 double h0=xa[1]-xa[0];
674 double h1=xa[2]-xa[1];
675 double h2=xa[3]-xa[2];
676 double A=2.0*(h0+h1);
681 gx[0]=3.0*((ya[2]-ya[1])/h1-(ya[1]-ya[0])/h0);
682 gx[1]=3.0*((ya[1]-ya[2])/h2-(ya[2]-ya[1])/h1);
684 det=3.0*(h0+h1)*(h0+h1);
685 this->c[1]=( A*gx[0]-B*gx[1])/det;
686 this->c[2]=(-B*gx[0]+A*gx[1])/det;
687 this->c[0]=this->c[2];
693 for (i=0; i < sys_size-1; i++) {
694 double h_i=xa[i+1]-xa[i];
695 double h_ip1=xa[i+2]-xa[i+1];
696 double ydiff_i=ya[i+1]-ya[i];
697 double ydiff_ip1=ya[i+2]-ya[i+1];
698 double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
699 double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
700 this->offdiag[i]=h_ip1;
701 this->diag[i]=2.0*(h_ip1+h_i);
702 this->g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
707 double h_i=xa[i+1]-xa[i];
708 double h_ip1=xa[1]-xa[0];
709 double ydiff_i=ya[i+1]-ya[i];
710 double ydiff_ip1=ya[1]-ya[0];
711 double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
712 double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
713 this->offdiag[i]=h_ip1;
714 this->diag[i]=2.0*(h_ip1+h_i);
715 this->g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
718 ubvector_range cp1(this->c,range(1,this->c.size()));
721 (this->diag,this->offdiag,this->g,cp1,sys_size,p5m);
722 this->c[0]=this->c[max_index];
729 #ifndef DOXYGEN_INTERNAL 736 (
const interp_cspline_peri<vec_t,vec2_t>&);
758 typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
759 typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
760 typedef boost::numeric::ublas::slice slice;
761 typedef boost::numeric::ublas::range range;
763 #ifndef DOXYGEN_INTERNAL 779 for(
size_t i=0;i<this->
sz-1;i++) {
781 double NE=fabs(umx[3+i]-umx[2+i])+fabs(umx[1+i]-umx[i]);
788 double h_i=(*this->
px)[i+1]-(*this->
px)[i];
789 double NE_next=fabs(umx[4+i]-umx[3+i])+
790 fabs(umx[2+i]-umx[1+i]);
791 double alpha_i=fabs(umx[1+i]-umx[i])/NE;
794 if (NE_next == 0.0) {
797 alpha_ip1=fabs(umx[2+i]-umx[1+i])/NE_next;
798 tL_ip1=(1.0-alpha_ip1)*umx[2+i]+alpha_ip1*umx[3+i];
800 b[i]=(1.0-alpha_i)*umx[1+i]+alpha_i*umx[2+i];
801 c[i]=(3.0*umx[2+i]-2.0*b[i]-tL_ip1)/h_i;
802 d[i]=(b[i]+tL_ip1-2.0*umx[2+i])/(h_i*h_i);
823 virtual void set(
size_t size,
const vec_t &xa,
const vec2_t &ya) {
827 " than "+
szttos(this->min_size)+
" in interp_akima::"+
831 if (size!=this->
sz) {
846 ubvector_range m(um,range(2,um.size()));
848 for (i=0;i<=size-2;i++) {
849 m[i]=(ya[i+1]-ya[i])/(xa[i+1]-xa[i]);
852 um[0]=3.0*m[0]-2.0*m[1];
854 m[this->
sz-1]=2.0*m[size-2]-m[size-3];
855 m[size]=3.0*m[size-2]-2.0*m[size-3];
857 akima_calc(xa,size,um);
863 virtual double eval(
double x0)
const {
866 size_t index=this->
svx.find_const(x0,cache);
868 double x_lo=(*this->
px)[index];
874 return (*this->
py)[index]+delx*(bb+delx*(cc+dd*delx));
878 virtual double deriv(
double x0)
const {
881 size_t index=this->
svx.find_const(x0,cache);
883 double x_lo=(*this->
px)[index];
889 return bb+delx*(2.0*cc+3.0*dd*delx);
898 size_t index=this->
svx.find_const(x0,cache);
900 double x_lo=(*this->
px)[index];
905 return 2.0*cc+6.0*dd*delx;
909 virtual double integ(
double aa,
double bb)
const {
911 size_t i, index_a, index_b;
914 if (((*this->
px)[0]<(*this->
px)[this->
sz-1] && aa>bb) ||
915 ((*this->
px)[0]>(*this->
px)[this->
sz-1] && aa<bb)) {
923 index_a=this->
svx.find_const(aa,cache);
924 index_b=this->
svx.find_const(bb,cache);
928 for(i=index_a; i<=index_b; i++) {
930 double x_lo=(*this->
px)[i];
931 double x_hi=(*this->
px)[i+1];
932 double y_lo=(*this->
py)[i];
937 if (i==index_a || i==index_b) {
938 double x1=(i==index_a) ? aa : x_lo;
939 double x2=(i==index_b) ? bb : x_hi;
940 result += this->
integ_eval(y_lo,b[i],c[i],d[i],x_lo,x1,x2);
942 result+=dx*(y_lo+dx*(0.5*b[i]+dx*(c[i]/3.0+0.25*d[i]*dx)));
946 double y_hi=(*this->
py)[i+1];
947 std::string str=((std::string)
"Interval of length zero ")+
950 " in interp_akima::integ().";
955 if (flip) result*=-1.0;
961 virtual const char *
type()
const {
return "interp_akima"; }
963 #ifndef DOXYGEN_INTERNAL 985 typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
986 typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
987 typedef boost::numeric::ublas::slice slice;
988 typedef boost::numeric::ublas::range range;
999 virtual const char *
type()
const {
return "interp_akima_peri"; }
1002 virtual void set(
size_t size,
const vec_t &xa,
const vec2_t &ya) {
1006 " than "+
szttos(this->min_size)+
" in interp_akima"+
1010 if (size!=this->
sz) {
1011 this->b.resize(size);
1012 this->c.resize(size);
1013 this->d.resize(size);
1014 this->um.resize(size+4);
1025 ubvector_range m(this->um,range(2,this->um.size()));
1028 for (
size_t i=0;i<=size-2;i++) {
1029 m[i]=(ya[i+1]-ya[i])/(xa[i+1]-xa[i]);
1032 this->um[0]=m[this->
sz-3];
1033 this->um[1]=m[this->
sz-2];
1037 this->akima_calc(xa,size,this->um);
1042 #ifndef DOXYGEN_INTERNAL 1048 (
const interp_akima_peri<vec_t,vec2_t>&);
1062 #ifdef O2SCL_NEVER_DEFINED 1069 typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
1070 typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
1071 typedef boost::numeric::ublas::slice slice;
1072 typedef boost::numeric::ublas::range range;
1074 #ifndef DOXYGEN_INTERNAL 1090 if ((x < 0 && y > 0) || (x > 0 && y < 0)) {
1110 virtual void set(
size_t size,
const vec_t &xa,
const vec2_t &ya) {
1114 " than "+
szttos(this->min_size)+
" in interp_steffen::"+
1118 if (size!=this->
sz) {
1123 y_prime.resize(size);
1137 double h0=(xa[1]-xa[0]);
1138 double s0=(ya[1]-ya[0]) / h0;
1144 for (
size_t i=1; i < (size-1); i++) {
1149 double hi=(xa[i+1]-xa[i]);
1150 double him1=(xa[i]-xa[i-1]);
1153 double si=(ya[i+1]-ya[i]) / hi;
1154 double sim1=(ya[i]-ya[i-1]) / him1;
1157 pi=(sim1*hi + si*him1) / (him1 + hi);
1160 y_prime[i]=(copysign(1.0,sim1)+copysign(1.0,si))*
1161 std::min(fabs(sim1),std::min(fabs(si),0.5*fabs(pi)));
1172 y_prime[size-1]=(ya[size-1]-ya[size-2])/
1173 (xa[size-1]-xa[size-2]);
1176 for (
size_t i=0; i < (size-1); i++) {
1177 double hi=(xa[i+1]-xa[i]);
1178 double si=(ya[i+1]-ya[i]) / hi;
1181 a[i]=(y_prime[i] + y_prime[i+1]-2*si) / hi / hi;
1182 b[i]=(3*si-2*y_prime[i]-y_prime[i+1]) / hi;
1191 virtual double eval(
double x0)
const {
1194 size_t index=this->
svx.find_const(x0,cache);
1195 double x_lo=(*this->
px)[index];
1196 double delx=x0-x_lo;
1199 double y = d[index]+delx*(c[index]+delx*(b[index]+delx*a[index]));
1208 size_t index=this->
svx.find_const(x0,cache);
1209 double x_lo=(*this->
px)[index];
1210 double delx=x0-x_lo;
1212 return c[index]+delx*(2.0*b[index]+delx*3.0*a[index]);
1221 size_t index=this->
svx.find_const(x0,cache);
1222 double x_lo=(*this->
px)[index];
1223 double delx=x0-x_lo;
1225 return 2.0*b[index]+delx*6.0*a[index];
1229 virtual double integ(
double al,
double bl)
const {
1231 size_t i, index_a, index_b;
1234 if (((*this->
px)[0]<(*this->
px)[this->
sz-1] && al>bl) ||
1235 ((*this->
px)[0]>(*this->
px)[this->
sz-1] && al<bl)) {
1243 index_a=this->
svx.find_const(al,cache);
1244 index_b=this->
svx.find_const(bl,cache);
1248 for(i=index_a; i<=index_b; i++) {
1250 double x_lo=(*this->
px)[i];
1251 double x_hi=(*this->
px)[i+1];
1252 double y_lo=(*this->
py)[i];
1253 double y_hi=(*this->
py)[i+1];
1254 double dx=x_hi-x_lo;
1258 double x1=(i == index_a) ? al-x_lo : 0.0;
1259 double x2=(i == index_b) ? bl-x_lo : x_hi-x_lo;
1260 result += (1.0/4.0)*a[i]*(x2*x2*x2*x2-x1*x1*x1*
x1)+
1261 (1.0/3.0)*b[i]*(x2*x2*x2-x1*x1*
x1)+
1262 (1.0/2.0)*c[i]*(x2*x2-x1*
x1)+d[i]*(x2-x1);
1265 std::string str=((std::string)
"Interval of length zero ")+
1268 " in interp_steffen::integ().";
1274 if (flip) result*=-1.0;
1280 virtual const char *
type()
const {
return "interp_steffen"; }
1282 #ifndef DOXYGEN_INTERNAL 1288 (
const interp_steffen<vec_t,vec2_t>&);
1321 #ifdef O2SCL_NEVER_DEFINED 1350 virtual void set(
size_t size,
const vec_t &x,
const vec2_t &y) {
1355 " than "+
szttos(this->min_size)+
" in interp_monotonic::"+
1363 if (this->
sz!=size) {
1365 Delta.resize(size-1);
1366 alpha.resize(size-1);
1367 beta.resize(size-1);
1378 for(
size_t i=0;i<size-1;i++) {
1379 Delta[i]=(y[i+1]-y[i])/(x[i+1]-x[i]);
1381 m[i]=(Delta[i]+Delta[i-1])/2.0;
1385 m[size-1]=Delta[size-2];
1388 for(
size_t i=0;i<size-1;i++) {
1396 for(
size_t i=0;i<size-1;i++) {
1397 alpha[i]=m[i]/Delta[i];
1398 beta[i]=m[i+1]/Delta[i];
1402 for(
size_t i=0;i<size-1;i++) {
1403 double norm2=alpha[i]*alpha[i]+beta[i]*beta[i];
1405 double tau=3.0/sqrt(norm2);
1406 m[i]=tau*alpha[i]*Delta[i];
1407 m[i+1]=tau*beta[i]*Delta[i];
1415 virtual double eval(
double x0)
const {
1418 size_t index=this->
svx.find_const(x0,cache);
1420 double x_lo=(*this->
px)[index];
1421 double x_hi=(*this->
px)[index+1];
1422 double y_lo=(*this->
py)[index];
1423 double y_hi=(*this->
py)[index+1];
1425 double t=(x0-x_lo)/h;
1426 double t2=t*t, t3=t2*t;
1428 double h00=2.0*t3-3.0*t2+1.0;
1429 double h10=t3-2.0*t2+t;
1430 double h01=-2.0*t3+3.0*t2;
1432 double interp=y_lo*h00+h*m[index]*h10+y_hi*h01+h*m[index+1]*h11;
1441 size_t index=this->
svx.find_const(x0,cache);
1443 double x_lo=(*this->
px)[index];
1444 double x_hi=(*this->
px)[index+1];
1445 double y_lo=(*this->
py)[index];
1446 double y_hi=(*this->
py)[index+1];
1448 double t=(x0-x_lo)/h;
1451 double dh00=6.0*t2-6.0*t;
1452 double dh10=3.0*t2-4.0*t+1.0;
1453 double dh01=-6.0*t2+6.0*t;
1454 double dh11=3.0*t2-2.0*t;
1455 double deriv=(y_lo*dh00+h*m[index]*dh10+y_hi*dh01+
1456 h*m[index+1]*dh11)/h;
1467 size_t index=this->
svx.find_const(x0,cache);
1469 double x_lo=(*this->
px)[index];
1470 double x_hi=(*this->
px)[index+1];
1471 double y_lo=(*this->
py)[index];
1472 double y_hi=(*this->
py)[index+1];
1474 double t=(x0-x_lo)/h;
1476 double ddh00=12.0*t-6.0;
1477 double ddh10=6.0*t-4.0;
1478 double ddh01=-12.0*t+6.0;
1479 double ddh11=6.0*t-2.0;
1480 double deriv2=(y_lo*ddh00+h*m[index]*ddh10+y_hi*ddh01+
1481 h*m[index+1]*ddh11)/h/h;
1487 virtual double integ(
double a,
double b)
const {
1489 size_t i, index_a, index_b;
1492 if (((*this->
px)[0]<(*this->
px)[this->
sz-1] && a>b) ||
1493 ((*this->
px)[0]>(*this->
px)[this->
sz-1] && a<b)) {
1501 index_a=this->
svx.find_const(a,cache);
1502 index_b=this->
svx.find_const(b,cache);
1506 for(i=index_a; i<=index_b; i++) {
1508 double x_hi=(*this->
px)[i+1];
1509 double x_lo=(*this->
px)[i];
1510 double y_lo=(*this->
py)[i];
1511 double y_hi=(*this->
py)[i+1];
1523 double t=(x_hi-x_lo)/h;
1524 double t2=t*t, t3=t2*t, t4=t3*t;
1526 double ih00=t4/2.0-t3+t;
1527 double ih10=t4/4.0-2.0*t3/3.0+t2/2.0;
1528 double ih01=-t4/2.0+t3;
1529 double ih11=t4/4.0-t3/3.0;
1530 double intres=h*(y_lo*ih00+h*m[i]*ih10+y_hi*ih01+
1535 std::string str=((std::string)
"Interval of length zero ")+
1538 " in interp_monotonic::integ().";
1544 if (flip) result*=-1.0;
1550 virtual const char *
type()
const {
return "interp_monotonic"; }
1552 #ifndef DOXYGEN_INTERNAL 1558 (
const interp_monotonic<vec_t,vec2_t>&);
1576 template<
class vec_t=boost::numeric::ublas::vector<
double>,
1579 #ifndef DOXYGEN_INTERNAL 1607 O2SCL_ERR((((std::string)
"Invalid interpolation type, ")+
1618 virtual double eval(
const double x0,
size_t n,
const vec_t &x,
1621 return itp->
eval(x0);
1625 virtual double deriv(
const double x0,
size_t n,
const vec_t &x,
1628 return itp->
deriv(x0);
1634 virtual double deriv2(
const double x0,
size_t n,
const vec_t &x,
1641 virtual double integ(
const double x1,
const double x2,
size_t n,
1642 const vec_t &x,
const vec2_t &y) {
1644 return itp->
integ(x1,x2);
1665 O2SCL_ERR((((std::string)
"Invalid interpolation type, ")+
1672 #ifndef DOXYGEN_INTERNAL 1703 template<
class vec_t=boost::numeric::ublas::vector<
double>,
1707 #ifndef DOXYGEN_INTERNAL 1729 const vec2_t &y,
size_t interp_type=
itp_cspline) {
1732 O2SCL_ERR((((std::string)
"Vector endpoints equal (value=")+
1752 O2SCL_ERR((((std::string)
"Invalid interpolation type, ")+
1754 "interp_vec::interp_vec().").c_str(),
exc_einval);
1769 void set(
size_t n,
const vec_t &x,
const vec2_t &y) {
1776 void set(
size_t n,
const vec_t &x,
1777 const vec2_t &y,
size_t interp_type) {
1780 O2SCL_ERR((((std::string)
"Vector endpoints equal (value=")+
1801 O2SCL_ERR((((std::string)
"Invalid interpolation type, ")+
1821 virtual double eval(
const double x0)
const {
1823 O2SCL_ERR(
"No vector set in interp_vec::eval().",
1826 return itp->
eval(x0);
1832 O2SCL_ERR(
"No vector set in interp_vec::operator().",
1835 return itp->
eval(x0);
1839 virtual double deriv(
const double x0)
const {
1841 O2SCL_ERR(
"No vector set in interp_vec::deriv().",
1844 return itp->
deriv(x0);
1850 virtual double deriv2(
const double x0)
const {
1852 O2SCL_ERR(
"No vector set in interp_vec::deriv2().",
1859 virtual double integ(
const double x1,
const double x2)
const {
1861 O2SCL_ERR(
"No vector set in interp_vec::integ().",
1864 return itp->
integ(x1,x2);
1869 return "interp_vec";
1872 #ifndef DOXYGEN_INTERNAL 1878 (
const interp_vec<vec_t,vec2_t> &it);
1889 public interp<double[n]> {
1895 :
interp<double[n]>(interp_type) {}
1915 size_t interp_type) :
1941 (
double level,
size_t n, vec_t &x, vec2_t &y) {
1944 O2SCL_ERR2(
"Need at least two data points in ",
1951 if (y[0]==level) count++;
1954 for(
size_t i=0;i<n-1;i++) {
1956 if ((y[i]<level && y[i+1]>=level) ||
1957 (y[i]>level && y[i+1]<=level)) {
1986 (
double level,
size_t n, vec_t &x, vec2_t &y, std::vector<double> &locs) {
1989 O2SCL_ERR2(
"Need at least two data points in ",
1998 locs.push_back(x[0]);
2002 for(
size_t i=0;i<n-1;i++) {
2004 if ((y[i]<level && y[i+1]>level) ||
2005 (y[i]>level && y[i+1]<level)) {
2008 double x0=x[i]+(x[i+1]-x[i])*(level-y[i])/(y[i+1]-y[i]);
2010 }
else if (y[i+1]==level) {
2011 locs.push_back(x[i+1]);
2021 template<
class ovec_t,
class vec2_t>
2025 for(
size_t i=0;i<n;i++) grid[i]=((
double)i);
2027 for(
size_t i=0;i<n;i++) dv[i]=oi.
deriv(((
double)i));
2034 template<
class ovec_t,
class vec2_t>
2038 for(
size_t i=0;i<n;i++) grid[i]=((
double)i);
2040 for(
size_t i=0;i<n;i++) dv[i]=oi.
deriv2(((
double)i));
2047 template<
class vec_t,
class vec2_t,
class vec3_t>
2051 for(
size_t i=0;i<n;i++) dv[i]=oi.
deriv(vx[i]);
2058 template<
class vec_t,
class vec2_t,
class vec3_t>
2062 for(
size_t i=0;i<n;i++) dv[i]=oi.
deriv(vx[i]);
2068 template<
class ovec_t>
2071 for(
size_t i=0;i<n;i++) grid[i]=((
double)i);
2073 return oi.
integ(0.0,((
double)(n-1)));
2089 template<
class vec_t,
class vec2_t>
2097 double total=si.
integ(x[0],x[n-1],n,x,y);
2105 template<
class vec_t,
class vec2_t,
class vec3_t>
2114 for(
size_t i=1;i<n;i++) {
2115 iy[i]=si.
integ(x[0],x[i],n,x,y);
2124 template<
class ovec_t>
2126 ovec_t &v,
size_t interp_type) {
2128 for(
size_t i=0;i<n;i++) grid[i]=((
double)i);
2130 return oi.
integ(0.0,x2);
2136 template<
class vec_t,
class vec2_t>
2144 double total=si.
integ(x[0],x2,n,x,y);
2185 (
double sum,
size_t n, vec_t &x, vec2_t &y,
double &lev,
2186 int boundaries=0,
int verbose=0,
bool err_on_fail=
true) {
2189 O2SCL_ERR2(
"Need at least two data points in ",
2196 std::vector<double>
x2, y2;
2198 if (boundaries==1) {
2201 x2[0]=x[0]-(x[1]-x[0])/1.0e6;
2203 for(
size_t i=0;i<n;i++) {
2208 }
else if (boundaries==2) {
2211 for(
size_t i=0;i<n;i++) {
2215 x2[n]=x[n-1]+(x[n-1]-x[n-2])/1.0e6;
2218 }
else if (boundaries==3) {
2221 x2[0]=x[0]-(x[1]-x[0])/1.0e6;
2223 for(
size_t i=0;i<n;i++) {
2227 x2[n+1]=x[n-1]+(x[n-1]-x[n-2])/1.0e6;
2233 for(
size_t i=0;i<n;i++) {
2243 vector_sort<ubvector,double>(ysort.size(),ysort);
2249 std::vector<double> ylist;
2250 for(
size_t i=0;i<ysort.size()-1;i++) {
2251 if (ysort[i]!=ysort[i+1]) {
2252 ylist.push_back((ysort[i+1]+3.0*ysort[i])/4.0);
2253 ylist.push_back((ysort[i+1]*3.0+ysort[i])/4.0);
2263 std::vector<double> xi, yi;
2267 for(
size_t k=0;k<ylist.size();k++) {
2268 double lev_tmp=ylist[k];
2269 std::vector<double> locs;
2271 if ((locs.size()%2)!=0) {
2274 std::cout << lev_tmp <<
" " << 0.0 <<
" " 2275 << locs.size() <<
" (fail)" << std::endl;
2278 double sum_temp=0.0;
2279 for(
size_t i=0;i<locs.size()/2;i++) {
2280 double x0=locs[2*i];
2281 double x1=locs[2*i+1];
2282 sum_temp+=itp.
integ(x0,x1,n2,x2,y2);
2284 xi.push_back(sum_temp);
2285 yi.push_back(lev_tmp);
2287 std::cout << lev_tmp <<
" " << sum_temp <<
" " 2288 << locs.size() <<
" ";
2289 for(
size_t i=0;i<locs.size();i++) {
2290 std::cout << locs[i] <<
" ";
2292 std::cout << std::endl;
2304 lev=itp2.
eval(sum,xi.size(),xi,yi);
2312 (
size_t n, vec_t &x, vec2_t &y,
double intl, std::vector<double> &locs,
2313 int boundaries=0,
int verbose=0,
bool err_on_fail=
true) {
2318 std::cout <<
"Total integral: " << total << std::endl;
2322 std::cout <<
"Target integral: " << intl << std::endl;
2327 boundaries,verbose,err_on_fail);
2330 O2SCL_ERR2(
"Failed to find a level which enclosed the ",
2331 "specified integral in vector_region_int().",
2337 std::cout <<
"Level from vector_invert: " << lev << std::endl;
2343 std::cout <<
"Locations from vector_find_level: ";
2344 for(
size_t i=0;i<locs.size();i++) {
2345 std::cout << locs[i];
2346 if (i!=locs.size()-1) {
2350 std::cout << std::endl;
2358 (
size_t n, vec_t &x, vec2_t &y,
double frac, std::vector<double> &locs,
2359 int boundaries=0,
int verbose=0,
bool err_on_fail=
true) {
2362 verbose,err_on_fail);
2372 (
size_t n, vec_t &x, vec2_t &y,
double frac,
double &low,
double &high,
2373 int boundaries=0,
int verbose=0,
bool err_on_fail=
true) {
2375 std::vector<double> locs;
2377 verbose,err_on_fail);
2378 if (locs.size()==0 || ret!=0) {
2380 O2SCL_ERR2(
"Zero level crossings or vector_region_fracint() ",
2399 (
size_t n, vec_t &x, vec2_t &y,
double frac,
double &low,
double &high,
2400 int boundaries=0,
int verbose=0,
bool err_on_fail=
true) {
2402 std::vector<double> locs;
2404 verbose,err_on_fail);
2405 if (locs.size()==0 || ret!=0) {
2407 O2SCL_ERR2(
"Zero level crossings or vector_region_int() ",
2419 #ifndef DOXYGEN_NO_O2NS 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 .
interp_vec(size_t n, const vec_t &x, const vec2_t &y, size_t interp_type=itp_cspline)
Create with base interpolation object it.
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 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 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 .
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.
#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 .
virtual double eval(double x0) const
Give the value of the function .
const vec2_t * py
Dependent vector.
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)
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.
Cubic spline for periodic boundary conditions.
Monotonicity-preserving interpolation (unfinished)
virtual double deriv2(const double x0, size_t n, const vec_t &x, const vec2_t &y)
Give the value of the second derivative .
ubvector Delta
Finite differences.
Cubic spline interpolation (GSL)
interp_array(size_t interp_type)
Create with base interpolation objects it and rit.
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.
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 .