/* S4ZERO.C
xref:
input:
output:
does: Scan for a zero of a function in an interval
*/
/************************************************************/

#include "mtb.h"

/************************************************************/

int s4zero (double f(double), double lower, double upper, int n,
            double *ap,  double *bp)

        /* scan the interval [lower, upper] for a zero of f() */
        /* evaluate f() at n points in the interval */
        /* return 1 if found a zero, 0 otherwise */
        /* [*ap, *bp] is the interval */
        /* returns 0 if n<1                         */
{
  double d, aa, bb, fa, fb;
  int i, found;

  if (n < 1)  return (0);
  if (lower>upper) {d=upper; upper=lower; lower=d;} /* swap */

  d = (upper - lower) / (n-1);
  aa = lower;
  fa = f(aa);

  for (i = 1; i<=n; i++) {
    bb = aa+d;
    fb = f(bb);
    if (found=(fa*fb <= 0)) break;
    aa = bb;  fa = fb;    /*  new a <- old b */
  }

  if (found) {*ap = aa; *bp = bb;}   /* success */
  return (found);
}
/************************************************************/

int s4zeros (double f(double), double lower, double upper, int n,
            double a[],  double b[] )

        /* scan the interval [lower, upper] for zeros of f() */
        /* evaluate f() at n points in the interval */
        /* return number of roots found (m) */
        /* a[i], b[i]  bracket i-th root, i=1,...,m */
        /* returns 0 if n<1                         */
{
  double d, aa, bb, fa, fb;
  int i, m;

  if (n < 1)  return (0);
  if (lower>upper) {d=upper; upper=lower; lower=d;} /* swap */

  d = (upper - lower) / (n-1);
  aa = lower;
  fa = f(aa);
  m = 0;

  for (i = 1; i<=n; i++) {
    bb = aa+d;
    fb = f(bb);
    if (fa*fb <= 0) {
        a[++m]=aa; b[m]=bb;
    }
    aa = bb;  fa = fb;    /*  new a <- old b */
  }

  return m;
}