double precision function zbrent(func,x1,x2,tol)
      implicit real*8 (a-h,o-z)
      parameter (itmax=200,eps=3.0D-08)
      external func

      a = x1
      b = x2
      fa = func(a)
      fb = func(b)
      zip = 0.0D+00

      if ((fa .gt. zip .and. fb .gt. zip) .or. (fa .lt. zip .and.
     &    fb .lt. zip)) pause 'root must be bracketed'

      c = b
      fc = fb
      do 11 iter = 1,itmax
         if ((fb .gt. zip .and. fc .gt. zip) .or. (fb .lt. zip .and.
     &       fc .lt. zip)) 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.0D+00*eps*dabs(b)+0.5D+00*tol
         xm = 0.5D+00*(c-b)
         if (dabs(xm) .le. tol1 .or. fb .eq. zip) then
            zbrent = b
            return
         endif
         if (dabs(e) .ge. tol1 .and. dabs(fa) .gt. dabs(fb)) then
            s = fb/fa
            if (a .eq. c) then
               p = 2.0D+00*xm*s
               q = 1.0D+00-s
            else
               q = fa/fc
               r = fb/fc
               p = s*(2.0D+00*xm*q*(q-r)-(b-a)*(r-1.0D+00))
               q = (q-1.0D+00)*(r-1.0D+00)*(s-1.0D+00)
            endif
            if (p .gt. zip) q = -q
            p = dabs(p)
            if (2.0D+00*p .lt. dmin1(3.0D+00*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)
 11   continue

      pause 'zbrent exceeding maximum iterations'
      zbrent = b

      return
      end

      double precision function rtsafe(funcd,x1,x2,xacc)
      implicit real*8 (a-h,o-z)
      parameter (maxit=100)
      external funcd

      call funcd(x1,fl,df)
      call funcd(x2,fh,df)
      if ((fl .gt. 0.0D+00 .and. fh .gt. 0.0D+00) .or.
     &    (fl .lt. 0.0D+00 .and. fh .lt. 0.0D+00)) pause
     &    'root must be bracketed in rtsafe'
      if (fl .eq. 0.0D+00) then
         rtsafe = x1
         return
      elseif (fh .eq. 0.0D+00) then
         rtsafe = x2
         return
      elseif (fl .lt. 0.0D+00) then
         xl = x1
         xh = x2
      else
         xh = x1
         xl = x2
      endif
      rtsafe = 0.5D+00*(x1+x2)
      dxold = dabs(x2-x1)
      dx = dxold
      call funcd(rtsafe,f,df)
      do 11 j = 1,maxit
         if (((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f) .gt. 0.0D+00
     &        .or. dabs(2.0D+00*f) .gt. dabs(dxold-df)) then
            dxold = dx
            dx = 0.5D+00*(xh-xl)
            rtsafe = xl+dx
            if (xl .eq. rtsafe) return
         else
            dxold = dx
            dx = f/df
            temp = rtsafe
            rtsafe = rtsafe-dx
            if (temp .eq. rtsafe) return
         endif
         if (dabs(dx) .lt. xacc) return
         call funcd(rtsafe,f,df)
         if (f .lt. 0.0D+00) then
            xl = rtsafe
         else
            xh = rtsafe
         endif
 11   continue
      pause 'rtsafe exceeding maximum iterations'
      return
      end

      subroutine zbrac(func,x1,x2,succes)
	implicit real*8 (a-h,o-z)
	parameter (factor=1.05D+00,ntry=50)
	logical succes
	external func

	if (dabs(x1-x2) .lt. 1.0D-06) pause 'no initial range'
	f1 = func(x1)
	f2 = func(x2)

	succes = .true.

	do 11 j = 1,ntry
         if (f1*f2 .lt. 0.0D+00) return
	   if (dabs(f1) .lt. dabs(f2)) then
	      x1 = x1+factor*(x1-x2)
	      f1 = func(x1)
	   else
	      x2 = x2+factor*(x2-x1)
	      f2 = func(x2)
	   endif
 11   continue

      succes = .false.

	return
	end