#define TNNearlyZero  1E-16

double UserFunction(double X,double (*FuncAddr)(double));

double _abs(double x);

void Secant(double         Guess1   ,
			double	       Guess2   ,
			double	       Tol      ,
			int	           MaxIter  ,
			double         *Root    ,
			double         *Value   ,
			int            *Iter    ,
			int  *Error   ,
			double (*FuncPtr)(double));

/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-             Input: Guess1, Guess2, Tol, MaxIter                          -*/
/*-            Output: Root, Value, Iter, Error                              -*/
/*-                                                                          -*/
/*-           Purpose: This unit provides a procedure for finding a single   -*/
/*-                    real root of a user specified function, given a       -*/
/*-                    specified initial guess.  The procedure implements    -*/
/*-                    the secant method for finding a single                -*/
/*-                    root of a function.                                   -*/
/*-                    The user must specify the desired tolerance           -*/
/*-                    in the answer.                                        -*/
/*-                                                                          -*/
/*-  Global Variables: Guess1  : real;    initial approx #1                  -*/
/*-                    Guess2  : real;    initial approx #2                  -*/
/*-                    Tol     : real;    tolerance in answer                -*/
/*-                    MaxIter : integer; maximum number of iterations       -*/
/*-                    Root    : real;    real part of calculated roots      -*/
/*-                    Value   : real;    value of the polynomial at root    -*/
/*-                    Iter    : real;    number of iterations it took       -*/
/*-                                       to find root                       -*/
/*-                    Error   : byte;    flags if something went wrong      -*/
/*-                                                                          -*/
/*-            Errors: 1: Iter >= MaxIter                                    -*/
/*-                    2: The slope was zero at some point                   -*/
/*-                    3: Tol <= 0                                           -*/
/*-                    4: MaxIter < 0                                        -*/
/*-                                                                          -*/
/*----------------------------------------------------------------------------*/

void Initial_Secant (double Guess1,
			 double	  Guess2,
			 double	  Tol,
			 int	  MaxIter,
			 double  *OldX,
			 double  *X,
			 double  *OldY,
			 double  *Y,
			 char  *Found,
			 int  *Iter,
			 int  *Error,
			 double (*FuncPtr)(double))

/*-------------------------------------------------------------*/
/*- Input:  Guess1, Guess2, Tol, MaxIter                      -*/
/*- Output: OldX, X, OldY, Y, Found, Iter, Error              -*/
/*-                                                           -*/
/*- This void sets the initial values of the above       -*/
/*- variables. If OldY or Y is zero, ) a root has been     -*/
/*- Found and Found = TRUE.  This void also checks       -*/
/*- the tolerance (Tol) and the maximum number of iterations  -*/
/*- (MaxIter) for errors.                                     -*/
/*-------------------------------------------------------------*/

{
  *Found = 0;
  *Iter = 0;
  *Error = 0;
  *OldX = Guess1;
  *X = Guess2;
  *OldY = UserFunction(*OldX, FuncPtr);
  *Y = UserFunction(*X, FuncPtr);
  if ( _abs(*OldY) <= TNNearlyZero )         /* OldX is the root */
	{
	  *X = *OldX;
	  *Y = *OldY;
	  *Found = 1;
	}
  else
	if ( _abs(*Y) <= TNNearlyZero )
	  *Found = 1;                         /*  X is the root  */
	else
	  if ( _abs((*OldY) - (*Y)) <= TNNearlyZero )
		*Error = 2;                  /* Slope of line is zero; no intercept */
  if ( Tol <= 0 )
	*Error = 3;
  if ( MaxIter < 0 )
	*Error = 4;
} /* function Initial_Secant  */

int TestForRoot_Secant (double X, double OldX, double y, double  Tol)

/*-----------------------------------------------------------------*/
/*- These are the stopping criteria.  Four different ones are     -*/
/*- provided.  If you wish to change the active criteria, simply  -*/
/*- comment off the current criteria (including the preceding OR) -*/
/*- and remove the comment brackets from the criteria (including  -*/
/*- the following OR) you wish to be active.                      -*/
/*-----------------------------------------------------------------*/

{
  return                              /*---------------------------*/
	(_abs(y) <= TNNearlyZero)          /*- Y = 0                   -*/
									  /*-                         -*/
		   ||                         /*-                         -*/
									  /*-                         -*/
	(_abs(X - OldX) < _abs(OldX*Tol))   /*- Relative change in X    -*/
									  /*-                         -*/
/*         ||                      */ /*-                         -*/
/*                                 */ /*-                         -*/
/*  (_abs(OldX - X) < Tol)          */ /*- Absolute change in X    -*/
/*                                 */ /*-                         -*/
/*         ||                      */ /*-                         -*/
/*                                 */ /*-                         -*/
/*  (_abs(y) <= Tol)                */ /*- Absolute change in Y    -*/
		   ;							  /*---------------------------*/

/*-----------------------------------------------------------------------*/
/*- The first criteria simply checks to see if ( the value of the         -*/
/*- function is zero.  You should probably always keep this criteria    -*/
/*- active.                                                             -*/
/*-                                                                     -*/
/*- The second criteria checks the relative error in X. This criteria   -*/
/*- evaluates the fractional change in X between interations. Note      -*/
/*- that X has been multiplied through the inequality to avoid divide   -*/
/*- by zero errors.                                                     -*/
/*-                                                                     -*/
/*- The third criteria checks the absolute difference in X between      -*/
/*- iterations.                                                         -*/
/*-                                                                     -*/
/*- The fourth criteria checks the absolute difference between          -*/
/*- the value of the function and zero.                                 -*/
/*-----------------------------------------------------------------------*/

} /* function TestForRoot_Secant  */


void Secant(double         Guess1   ,
			double	       Guess2   ,
			double	       Tol      ,
			int	           MaxIter  ,
			double         *Root    ,
			double         *Value   ,
			int            *Iter    ,
			int  *Error   ,
			double (*FuncPtr)(double))
{

  char Found;            /* Flags that a root has been Found */
  double  OldX, OldY, X, Y, NewX, NewY;
					   /* Iteration variables */


  /* function Secant */
  Initial_Secant (Guess1, Guess2, Tol, MaxIter, &OldX, &X, &OldY, &Y,
		  &Found, Iter, Error, FuncPtr);

  while ( !Found && !(*Error) && (*Iter<MaxIter) )
  {
	(*Iter)++;
	NewX = X - Y * (X - OldX) / (Y - OldY);    /* The secant formula */
	NewY = UserFunction(NewX, FuncPtr);
	Found = TestForRoot_Secant (NewX, OldX, NewY, Tol);
	OldX = X;
	OldY = Y;
	X = NewX;
	Y = NewY;
  }
  *Root = X;
  *Value = Y;
  if ( !(Found) && !(*Error) && ((*Iter) >= MaxIter) )
	*Error = 1;
} /* function Secant */

double UserFunction(double X, double (*ProcAddr)(double))
{
 return ProcAddr(X);
}

double _abs(double x)
{
 return ( x < 0 ? -x : x);
}
