C+
C NAME:
C	nrZBrentd
C PURPOSE:
C	Find node (zero) of function of one variable
C CATEGORY:
C	Math: finding roots
C CALLING SEQUENCE:
	double precision function nrZBrentd(FUNC,X1,X2,TOL)
C INPUTS:
C	FUNC		name of the function
C			should be declared as double precision external function in the
C			calling program, i.e.
C
C			double precision	FUNC
C		 	external		FUNC
C
C	X1,X2		double precision	x-values which bracket the required node,
C						i.e. FUNC(X1)*FUNC(X2) >= 0
C	TOL		double precision	tolerance in position of node
C OUTPUTS:
C	nrZBrentd	double precision	x-value (with accuracy TOL) which makes function value zero
C CALLS:
C	Say, nrZBrak, nrZBrac
C SEE ALSO:
C	nrZBrent
C RESTRICTIONS:
C >	The node must be bracketed by X1 and X2, i.e. FUNC(X1)*FUNC(X2) >= 0.
C PROCEDURE:
C	See Press and Teukolsky, Numerical Recipes
C-
	implicit double precision (A-H,O-Z)
	parameter (ITMAX = 100)
	parameter (EPS   = 3E-8)

	A = X1
	B = X2
	FA = FUNC(A)
	FB = FUNC(B)
	if (FB*FA .gt. 0) call Say('nrZBrentd','E','Stop','root is not bracketed')

	FC = FB
	do ITER=1,ITMAX
	    if (FB*FC .gt. 0.) then
		C = A		! Rename A,B,C and adjust bounding interval
		FC = FA
		D = B-A
		E = D
	    end if
	    if (abs(FC) .lt. abs(FB)) then
		A = B
		B = C
		C = A
		FA = FB
		FB = FC
		FC = FA
	    end if
	    TOL1 = 2.*EPS*abs(B)+.5*TOL		! Convergence check
	    XM = .5*(C-B)
	    if (abs(XM) .le. TOL1 .or. FB .eq. 0.) then
		nrZBrentd = B
		return
	    end if
	    if (abs(E) .ge. TOL1 .and. abs(FA) .gt. abs(FB)) then
		S = FB/FA		! Attempt inverse quadratic interpolation
		if (A .eq. C) then
		    P = 2.*XM*S
		    Q = 1.-S
		else
		    Q = FA/FC
		    R = FB/FC
		    P = S*(2.*XM*Q*(Q-R)-(B-A)*(R-1.))
		    Q = (Q-1.)*(R-1.)*(S-1.)
		end if
		if (P .gt. 0.) Q = -Q	! Check whether in bounds
		P = abs(P)
		if (2.*P .lt. min(3.*XM*Q-abs(TOL1*Q),abs(E*Q))) then
		    E = D		! Accept interpolation
		    D = P/Q
		else
		    D = XM		! Interpolation failed, use bisection
		    E = D
		end if
	    else			! Bounds decreasing too slowly, use bisection
		D = XM
		E = D
	    end if
	    A = B
	    FA = FB
	    if (abs(D) .gt. TOL1) then	! Evaluate new trial root
		B = B+D
	    else
		B = B+sign(TOL1,XM)
	    end if
	    FB = FUNC(B)
	end do

	nrZBrentd = B

	return
	end