1 #ifdef HAVE_STDLIB_H
2 #include <stdlib.h>
3 #endif
4 #include <math.h>
5 #include "yagi.h"
6
7 #define ITMAX 100
8 #define EPS 3.0e-8
9
zbrent(double (* error_3dB)(double x),double x1,double x2,double tol)10 double zbrent(double (*error_3dB) (double x),double x1,double x2, double tol)
11 {
12 int iter;
13 double a=x1,b=x2,c,d,e,min1,min2;
14 double fa=(*error_3dB)(a),fb=(*error_3dB)(b),fc,p,q,r,s,tol1,xm;
15 void nrerror();
16 /* printf("x1=%f fa=%f x2=%f fb=%f\n", x1,fa,x2,fb); */
17 if (fb*fa > 0.0) nrerror("Root must be bracketed in ZBRENT - try the -e or -h options to output");
18 fc=fb;
19 for (iter=1;iter<=ITMAX;iter++) {
20 if (fb*fc > 0.0) {
21 c=a;
22 fc=fa;
23 e=d=b-a;
24 }
25 if (fabs(fc) < fabs(fb)) {
26 a=b;
27 b=c;
28 c=a;
29 fa=fb;
30 fb=fc;
31 fc=fa;
32 }
33 tol1=2.0*EPS*fabs(b)+0.5*tol;
34 xm=0.5*(c-b);
35 if (fabs(xm) <= tol1 || fb == 0.0) return b;
36 if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
37 s=fb/fa;
38 if (a == c) {
39 p=2.0*xm*s;
40 q=1.0-s;
41 } else {
42 q=fa/fc;
43 r=fb/fc;
44 p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
45 q=(q-1.0)*(r-1.0)*(s-1.0);
46 }
47 if (p > 0.0) q = -q;
48 p=fabs(p);
49 min1=3.0*xm*q-fabs(tol1*q);
50 min2=fabs(e*q);
51 if (2.0*p < (min1 < min2 ? min1 : min2)) {
52 e=d;
53 d=p/q;
54 } else {
55 d=xm;
56 e=d;
57 }
58 } else {
59 d=xm;
60 e=d;
61 }
62 a=b;
63 fa=fb;
64 if (fabs(d) > tol1)
65 b += d;
66 else
67 b += (xm > 0.0 ? fabs(tol1) : -fabs(tol1));
68 fb=(*error_3dB)(b);
69 }
70 nrerror("Maximum number of iterations exceeded in ZBRENT");
71 }
72
73 #undef ITMAX
74 #undef EPS
75