/* Adaptado de Numerical Recipes in C */

#include <stdio.h>
#include <math.h>
#include "protos.h"

#define JMAX 40 /* Maximum allowed number of bisections. */

float rtbis(float x1, float x2, float xacc)
/*
  Using bisection, find the root of a function func known to lie between x1 and x2. 
  The root, returned as rtbis, will be refined until its accuracy is +/- xacc.
*/
{
  int j;
  float dx,f,fmid,xmid,rtb;
  f=func(x1);
  fmid=func(x2);
  if (f*fmid >= 0.0) printf("Raiz deve estar encurralada para a bissecao em rtbis\n");
  rtb = f < 0.0 ? (dx=x2-x1,x1) : (dx=x1-x2,x2); /* Orient the search so that f>0 lies at x+dx. */
  for (j=1;j<=JMAX;j++) { 
    fmid=func(xmid=rtb+(dx *= 0.5)); /* Bisection loop */
    if (fmid <= 0.0) rtb=xmid;
    if (fabs(dx) < xacc || fmid == 0.0) return rtb;
  }
  printf("Numero maximo de bissecoes atingido!\n");
  return 0.0; /* Never get here. */
}