135 lines
3.0 KiB
Fortran
135 lines
3.0 KiB
Fortran
c****************************************************************
|
|
|
|
function zbrent(func,x1,x2,tol)
|
|
|
|
c****************************************************************
|
|
c**
|
|
c** FILE NAME: zbrent.f
|
|
c**
|
|
c** DATE WRITTEN: 8/3/90
|
|
c**
|
|
c** PROGRAMMER:Scott Hensley
|
|
c**
|
|
c** FUNCTIONAL DESCRIPTION: The subroutine is a routine taken
|
|
c** from "Numerical Recipes" for finding the zero of a function
|
|
c** to within a given tolerance. This routine requires that
|
|
c** the function and desired tolerance be passed in as arguements
|
|
c** along with two values which bracket the zero (i.e. the
|
|
c** function values at these points must be of opposite sign).
|
|
c**
|
|
c** ROUTINES CALLED:none
|
|
c**
|
|
c** NOTES: none
|
|
c**
|
|
c** UPDATE LOG:
|
|
c**
|
|
c*****************************************************************
|
|
|
|
IMPLICIT real*8 (a-h,o-z)
|
|
|
|
c INPUT VARIABLES:
|
|
real*8 x1 !bracketing parameter
|
|
real*8 x2 !bracketing parameter
|
|
real*8 tol !tolerance
|
|
real*8 zbrent !value of zero
|
|
|
|
c OUTPUT VARIABLES:see input
|
|
|
|
c LOCAL VARIABLES:
|
|
parameter (itmax=100,eps=3.d-8)
|
|
|
|
c PROCESSING STEPS:
|
|
|
|
c for details see "Numerical Recipes"
|
|
|
|
a = x1
|
|
b = x2
|
|
fa = func(a)
|
|
fb = func(b)
|
|
|
|
c write(*,*) 'fa fb ',a,b,fa,fb
|
|
if(fb*fa .gt. 0)then
|
|
zbrent=1.e9 !no zero in the range
|
|
return
|
|
end if
|
|
|
|
fc = fb
|
|
|
|
do iter=1,itmax
|
|
|
|
if(fb*fc .gt. 0)then
|
|
c = a
|
|
fc = fa
|
|
d = b - a
|
|
e = d
|
|
endif
|
|
|
|
if(dabs(fc) .lt. dabs(fb))then
|
|
a = b
|
|
b = c
|
|
c = a
|
|
fa = fb
|
|
fb = fc
|
|
fc = fa
|
|
endif
|
|
|
|
tol1 = 2*eps*dabs(b) + .5*tol
|
|
xm = .5*(c-b)
|
|
|
|
if(dabs(xm) .le. tol1 .or. fb .eq. 0)then
|
|
zbrent = b
|
|
return
|
|
endif
|
|
|
|
if(dabs(e) .ge. tol1 .and. dabs(fa) .gt. dabs(fb))then
|
|
|
|
sa = fb/fa
|
|
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)
|
|
endif
|
|
|
|
if(p .gt. 0) q = -q
|
|
p = dabs(p)
|
|
|
|
if(2*p .lt. dmin1(3.*xm*q - dabs(tol1*q),dabs(e*q)))then
|
|
e = d
|
|
d = p/q
|
|
else
|
|
d = xm
|
|
e = d
|
|
endif
|
|
|
|
else
|
|
|
|
d = xm
|
|
e = d
|
|
|
|
endif
|
|
|
|
a = b
|
|
fa = fb
|
|
|
|
if(dabs(d) .gt. tol1)then
|
|
b = b + d
|
|
else
|
|
b = b + sign(tol1,xm)
|
|
endif
|
|
|
|
fb = func(b)
|
|
|
|
enddo
|
|
|
|
pause 'it ain t gonna work this time'
|
|
|
|
zbrent = b
|
|
|
|
return
|
|
|
|
end
|