/************************************************************************
Leah Jager
February 2003

This program is the generic recipe for the Van Wijngaarden-Dekker-Brent
method from "Numerical Recipes in C", p361. 

Press, W. H., Teukolsky, S. A., Vetterling, W. T., Flannery, B. P. 
     "Numerical Recipes in C." Cambridge University Press, 2002. 

The function zbrent finds the root of a function func known to lie between 
x1 and x2 using Brent's method.  The root is returned as zbrent, and is 
refined until accuracy is better than tol.                              
***************************************************************************/

#include <stdio.h>
#include <math.h>
/*  #include "nrutil.h"  replace with inline declaration for SIGN  */
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
#define ITMAX 100
#define EPS 3.0e-8


void nrerror(char error_text[])
     /*  Numerical Recipes standard error handler  */
{
  fprintf(stderr, "Numerical Recipes run-time error...\n");
  fprintf(stderr, "%s\n", error_text);
  fprintf(stderr, "...now exiting to system...\n");
  exit(1);
}

/***********************************************************************/
double zbrent(double (func)(double, double, double, int, double *), double nuisance, double nuisance2, int n, double *data,double x1, double x2, double tol)
{
  int iter;
  
  double a=x1,b=x2,c=x2,d,e,min1,min2;
  /*
  printf("\n a and b are %lf and %lf", a,b);
  */
  double fa=(func)(a,nuisance,nuisance2,n,data);

  double fb=(func)(b,nuisance,nuisance2,n,data),fc,p,q,r,s,tol1,xm;
  /*
  printf("\n fa is %lf, fb is %lf", fa, fb);
  */
  
  if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)){
    /*   
    printf("\t here we are returning -1");   
    */
    return -1.0;
  }
  /*
    nrerror("Root must be bracketed in zbrent");
  */ 
  fc=fb;
  for (iter=1;iter<=ITMAX;iter++) {
   
    if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
      c=a;                           
      fc=fa;                       
      e=d=b-a;
    }

    if (fabs(fc) < fabs(fb)) {
      a=b;
      b=c;
      c=a;
      fa=fb;
      fb=fc;
      fc=fa;
    }

    tol1=2.0*EPS*fabs(b)+0.5*tol;     
    xm=0.5*(c-b);

    if (fabs(xm) <= tol1 || fb == 0.0) { 
      return b;
    }
    if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
      s=fb/fa;                 
      if (a == c) {
	p=2.0*xm*s;
	q=1.0-s;
      } else {
	q=fa/fc;
	r=fb/fc;
	p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
	q=(q-1.0)*(r-1.0)*(s-1.0);
      }

      if (p > 0.0) q = -q;      
      p=fabs(p);
      min1=3.0*xm*q-fabs(tol1*q);
      min2=fabs(e*q);

      if (2.0*p < (min1 < min2 ? min1 : min2)) {
	e=d;                      
	d=p/q;
      } else {
	d=xm;                     
	e=d;
      }
    }
    else {
      d=xm;                 
      e=d;
    }
    a=b;                   
    fa=fb;
    if (fabs(d) > tol1)      
      b += d;
    else 
      b += SIGN(tol1,xm);
    fb=(*func)(b,nuisance,nuisance2,n,data);
  }
  nrerror("Maximum number of iterations exceeded in zbrent");										    
  return 0.0;            
}