/* newrap.c: Newton-Raphson iteration */
/************************************************************/
#include <math.h>
#include "mtb.h"
/************************************************************/

int NewRap (double f(double),    double fp(double),
            double *xp,  int monitor, double eps, int maxiter)
{
  double fx, fpx;
  int i, converged, fpnearzero;
  const double nearzero = 1E-20;

  for (i=0; i <  maxiter; i++) {
    fx  = f(*xp);
    fpx = fp(*xp) ;
    if (monitor)
      printf ("i=%3d, x=%20.16f, f=%20.16f, fp=%20.16f\n",
               i, *xp, fx, fpx);
    if ( converged = (fabs(fx) < eps ) ) break;
    if ( fpnearzero = (fabs(fpx) < nearzero) ) break;
    *xp = *xp - fx / fpx;
  }

  if (converged)  return NR_CONV;
  else if (i == maxiter) return NR_NOCONV;
  else return NR_FP0;
}

/************************************************************/
void NewRapMsg (int rc)
{
  const char*  NR_MSG [] =
    {"NewRap converged",
     "NewRap did not converge",
     "Derivative near zero in NewRap",
     "Invalid return code" };

  const char myname[] = "NewRapMsg:";
  const int n = sizeof(NR_MSG) / sizeof(NR_MSG[0]) - 1;

  if (rc < 0 || rc >= n)
       printf ("%s %s, rc=%d\n", myname, NR_MSG[n], rc);
  else printf ("%s %s\n",        myname, NR_MSG[rc]   );
}
/************************************************************/