double fmin( string _s, double _u, double _v, double _tol,
  matrix _p1, matrix _p2, matrix _p3, matrix _p4, matrix _p5)
{
/*
 * Minimization of function s using interval halving.
 *
 * Return value is position of minimum s(b) in the interval u..v
 * within relative tolerance tol.
 *
 * p1..p5 are up to 5 extra parameters to be passed directly
 * to the function being minimized.
 *
 * R. Perry, June 1998
 */

  double b1, b2, b3, b4, b5;	/* positions and */
  double s1, s2, s3, s4, s5;	/*  function values */

  double bmin, min;		/* position and value of current minimum */
  double h;			/* distance between points b1..b2, etc. */
  int jmin;			/* index of current minimum */

/*
 * The interval being searched is  b1 <= b <= b5, with b3 being a point
 * in the middle of the interval.  The variable jmin indicates whether
 * the current minimum is at b3 (jmin=3) or at one of the interval endpoints
 * b1 (jmin=1) or b5 (jmin=5).
 *
 * If the current minimum is at b3, then the function is evaluated at
 * b2, a point midway between b1 & b3.  If s(b2) is lower than the current
 * minimum, the interval is reduced to b1..b3; otherwise, the function is
 * evaluated at b4, a point midway between b3 & b5.  If s(b4) is lower than
 * the current minimum, the interval is reduced to b3..b5.  If neither
 * b2 nor b4 produce a lower minimum, the interval is reduced to b2..b4.
 *
 * If the current minimum is at one of the interval endpoints b1 or b5
 * then the interval is reduced to b1..b3 or b3..b5 and the function is
 * evaluated at the midpoint of the new interval.
 */

/**** initialize ****/

	h = ( _v - _u ) / 2;
	b1 = _u;
	b3 = _u + h;
	b5 = _v;

	s1 = feval( _s, b1, _p1, _p2, _p3, _p4, _p5);
	s3 = feval( _s, b3, _p1, _p2, _p3, _p4, _p5);
	s5 = feval( _s, b5, _p1, _p2, _p3, _p4, _p5);

	min = s1;
	jmin = 1;
	bmin = b1;

	if( s3 < min) {
	  min = s3;
	  jmin = 3;
	  bmin = b3;
	}

	if( s5 < min) {
	  min = s5;
	  jmin = 5;
	  bmin = b5;
	}

	h /= 2;

/**** main loop ****/

	while( abs(h) >= _tol*(abs(bmin)+_tol) ) {

	  if( jmin == 3) { /* current minimum is at midpoint */

	    b2 = b1 + h;
	    s2 = feval( _s, b2, _p1, _p2, _p3, _p4, _p5);

	    if( s2 < min ) { /* reduce interval to b1..b3 */

	      min = s2;
	      bmin = b2;
	      b5 = b3;
	      b3 = b2;

	    } else { /* try b4 */

	      b4 = b3 + h;
	      s4 = feval( _s, b4, _p1, _p2, _p3, _p4, _p5);

	      if( s4 < min) { /* reduce interval to b3..b5 */

		min = s4;
		bmin = b4;
		b1 = b3;
		b3 = b4;

	      } else { /* reduce interval to b2..b4 */

		b1 = b2;
		b5 = b4;

	      }

	    }

	  } else { /* current minimum is at an endpoint */

	    if( jmin == 1) /* reduce interval to b1..b3 */
	      b5 = b3;
	    else /* reduce interval to b3..b5 */
	      b1 = b3;

	    b3 = b1 + h; /* new midpoint */
	    s3 = feval( _s, b3, _p1, _p2, _p3, _p4, _p5);

	    if( s3 < min) { /* new minimum is at midpoint */
	      min = s3;
	      jmin = 3;
	      bmin = b3;
	    }

	  }

	  h /= 2;

	}

	return  bmin;
}
