71 typedef std::vector<std::vector<double> >
dMatrix;
73 typedef std::vector<std::vector<int> >
iMatrix;
75 typedef std::vector<std::vector<bool> >
bMatrix;
77 double Interpola(
double vizq,
double vder,
double axid,
double xif);
115 return x * x * x * x;
128 return sqrt(sqrt(x));
141 return sqrt(
pow3(x));
157 template<
class T,
class U>
158 inline T poww(T x, U y) {
159 return abs(x) > std::numeric_limits<T>::epsilon() ? pow(x, y) : 0;
163 inline T sqrtw(T x) {
164 return abs(x) > std::numeric_limits<T>::epsilon() ? sqrt(x) : 0;
170 return pow(x, 1 / 3) > std::numeric_limits<T>::epsilon() ? pow(x, 1 / 3) : 0;
175 inline bool DoubEqZero(T x) {
176 return abs(x) > std::numeric_limits<T>::epsilon() ? false :
true;
179 double QuadraticEqP(
double A,
double B,
double C);
181 inline double QuadraticEqN(
double A,
double B,
double C);
184 inline T SQR(
const T a) {
189 inline const T &Max(
const T &a,
const T &b) {
190 return b > a ? (b) : (a);
193 float Max(
const double &a,
const float &b);
195 float Max(
const float &a,
const double &b);
198 inline const T &Min(
const T &a,
const T &b) {
199 return b < a ? (b) : (a);
202 inline float Min(
const double &a,
const float &b);
204 inline float Min(
const float &a,
const double &b);
207 inline T Sign(
const T &a,
const T &b) {
208 return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
211 inline float Sign(
const float &a,
const double &b);
213 inline float Sign(
const double &a,
const float &b);
216 inline void Swap(T &a, T &b) {
223 inline T MaxComponent(std::vector<T> &x) {
225 for(
Uint i = 1; i < x.size(); i++) {
233 inline T MinComponent(std::vector<T> &x) {
235 for(
Uint i = 1; i < x.size(); i++) {
271 void operator()(
double X,
double Y);
282 void operator()(
double Ang,
double Mod);
286 int n, mm, jsav, cor, dj;
287 const double *xx, *yy;
293 double interp(
double x);
295 int locate(
const double x);
296 int hunt(
const double x);
298 double virtual rawinterp(
int jlo,
double x) = 0;
308 double rawinterp(
int j,
double x);
320 void sety2(
const double *xv,
const double *yv);
321 double rawinterp(
int j,
double x);
332 double rawinterp(
int j,
double x);
352 inline double zbrent(T& func,
const double& x1,
const double& x2,
const double& tol) {
353 const int ITMAX = 100;
354 const double EPS = std::numeric_limits<double>::epsilon();
370 if((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) {
371 if(fabs(fa) < fabs(fb)) {
375 else if(fabs(fa) > fabs(fb)) {
384 for(
int iter = 0; iter < ITMAX; iter++) {
385 if((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
390 if(fabs(fc) < fabs(fb)) {
398 tol1 = 2.0 * EPS * fabs(b) + 0.5 * tol;
400 if(fabs(xm) <= tol1 || fb == 0.0) {
404 if(fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
412 p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0));
413 q = (q - 1.0) * (r - 1.0) * (s - 1.0);
418 double min1 = 3.0 * xm * q - fabs(tol1 * q);
419 double min2 = fabs(e * q);
420 if(2.0 * p < (min1 < min2 ? min1 : min2)) {
459 inline double FindRoot(T & func,
const double x1,
const double x2) {
460 return zbrent(func, x1, x2, 1e-10);
464 inline bool zbrac(T & func,
double & x1,
double & x2) {
465 const int NTRY = 200;
466 const double FACTOR = 0.1;
468 throw(
"Bad initial range in zbrac");
469 double f1 = func(x1);
470 double f2 = func(x2);
471 for(
int j = 0; j < NTRY; j++) {
474 if(fabs(f1) < fabs(f2))
475 f1 = func(x1 += FACTOR * (x1 - x2));
477 f2 = func(x2 += FACTOR * (x2 - x1));
482 inline bool zbrac2(T & func,
double & x1,
double & x2,
const double & min,
const double & max) {
483 const int NTRY = 200;
486 throw(
"Bad initial range in zbrac");
493 double f1 = func(x1);
494 double f2 = func(x2);
495 for(
int j = 0; j < NTRY; j++) {
498 if(fabs(f1) < fabs(f2))
499 f1 = func(x1 = (0.9 * x1 + 0.1 * x1lim));
501 f2 = func(x2 = (0.9 * x2 + 0.1 * x2lim));
503 std::cout <<
"Do not exist solution" << std::endl;
504 std::cout <<
"x1 :" << x1 <<
" f1: " << f1 <<
" x1max: " << x1lim << std::endl;
505 std::cout <<
"x2 :" << x2 <<
" f2: " << f2 <<
" x2max: " << x2lim << std::endl;
512 inline void zbrak(T & fx,
const double x1,
const double x2,
const int n,
dVector & xb1,
dVector & xb2,
int & nroot) {
517 double dx = (x2 - x1) / n;
519 double fp = fx(x += dx);
520 for(
int i = 0; i < n; i++) {
521 double fc = fx(x += dx);
526 dVector tempvec1(xb1), tempvec2(xb2);
529 for(
int j = 0; j < nb; j++) {
530 xb1[j] = tempvec1[j];
531 xb2[j] = tempvec2[j];
541 inline double rtbis(T & func,
const double x1,
const double x2,
const double xacc) {
544 double dx, xmid, rtb;
547 double fmid = func(x2);
548 if(f * fmid >= 0.0) {
549 if((fabs(f) < fabs(fmid)) && fabs(f) < xacc) {
551 }
else if((fabs(f) > fabs(fmid)) && fabs(fmid) < xacc)
553 throw(
"Root must be bracketed for bisection in rtbis");
555 rtb = f < 0.0 ? (dx = x2 - x1, x1) : (dx = x1 - x2, x2);
556 for(
int j = 0; j < JMAX; j++) {
557 fmid = func(xmid = rtb + (dx *= 0.5));
560 if(fabs(dx) < xacc || fmid == 0.0) {
564 throw(
"Too many bisections in rtbis");
568 inline double rtflsp(T & func,
const double x1,
const double x2,
const double xacc) {
569 const int MAXIT = 1000;
573 double fl = func(x1);
574 double fh = func(x2);
576 throw(
"Root must be bracketed in rtflsp");
586 for(
int j = 0; j < MAXIT; j++) {
587 double rtf = xl + dx * fl / (fl - fh);
588 double f = func(rtf);
599 if(fabs(del) < xacc || f == 0.0) {
603 throw(
"Maximum number of iterations exceeded in rtflsp");
607 inline double rtsec(T & func,
const double x1,
const double x2,
const double xacc) {
608 const int MAXIT = 100;
612 double fl = func(x1);
614 if(fabs(fl) < fabs(f)) {
622 for(
int j = 0; j < MAXIT; j++) {
623 double dx = (xl - rts) * f / (f - fl);
628 if(fabs(dx) < xacc || f == 0.0) {
632 throw(
"Maximum number of iterations exceede in rtsec");
636 inline double zriddr(T & func,
const double x1,
const double x2,
const double xacc) {
637 const int MAXIT = 100;
638 double fl = func(x1);
639 double fh = func(x2);
640 if((fl > 0.0 && fh < 0.0) || (fl < 0.0 && fh > 0.0)) {
643 double ans = -9.99e99;
644 for(
int j = 0; j < MAXIT; j++) {
645 double xm = 0.5 * (xl + xh);
646 double fm = func(xm);
647 double s = sqrt(fm * fm - fl * fh);
651 double xnew = xm + (xm - xl) * ((fl >= fh ? 1.0 : -1.0) * fm / s);
652 if(fabs((
double)(xnew - ans <= xacc))) {
656 double fnew = func(ans);
660 if(Sign(fm, fnew) != fm) {
665 }
else if(Sign(fl, fnew) != fl) {
668 }
else if(Sign(fh, fnew) != fh) {
672 throw(
"never get here.");
673 if(fabs(xh - xl) <= xacc) {
677 throw(
"zriddr exceed maximum iterations");
683 throw(
"root must be bracketed in zriddr.");
688 inline double rtnewt(T & funcd,
const double x1,
const double x2,
const double xacc) {
689 const int JMAX = 20.;
690 double rtn = 0.5 * (x1 + x2);
691 for(
int j = 0; j < JMAX; j++) {
692 double f = funcd(rtn);
693 double df = funcd.df(rtn);
696 if((x1 - rtn) * (rtn - x2) < 0.0)
697 throw(
"Jumped out of brackets in rtnewt");
701 throw(
"Maximum number of iterations exceede in rtnewt");
705 inline double rtsafe(T & funcd,
const double x1,
const double x2,
const double xacc) {
706 const int MAXIT = 100;
710 double fl = funcd(x1);
711 double fh = funcd(x2);
712 if((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0)) {
713 throw(
"Root mus be bracketed in rt safe");
726 double rts = 0.5 * (x1 + x2);
727 double dxold = fabs(x2 - x1);
729 double f = funcd(rts);
730 double df = funcd.df(rts);
731 for(
int j = 0; j < MAXIT; j++) {
732 if((((rts - xh) * df - f) * ((rts - xl) * df - f) > 0.0) || (fabs(2.0 * f) > fabs(dxold * df))) {
734 dx = 0.5 * (xh - xl);
748 double f = funcd(rts);
749 double df = funcd.df(rts);
755 throw(
"Maximum number of iterations exceede in rtsafe");