C R's  optimize() :   function	fmin(ax,bx,f,tol)
C =    ==========		~~~~~~~~~~~~~~~~~
c
c      an approximation  x  to the point where  f  attains a minimum  on
c  the interval  (ax,bx)  is determined.
c
c  input..
c
c  ax    left endpoint of initial interval
c  bx    right endpoint of initial interval
c  f     function subprogram which evaluates  f(x)  for any  x
c        in the interval  (ax,bx)
c  tol   desired length of the interval of uncertainty of the final
c        result (.ge.0.)
c
c  output..
c
c  fmin  abcissa approximating the point where  f  attains a
c        minimum
c
c      the method used is a combination of  golden  section  search  and
c  successive parabolic interpolation.  convergence is never much slower
c  than  that  for  a  fibonacci search.  if  f  has a continuous second
c  derivative which is positive at the minimum (which is not  at  ax  or
c  bx),  then  convergence  is  superlinear, and usually of the order of
c  about  1.324....
c      the function  f  is never evaluated at two points closer together
c  than  eps*abs(fmin)+(tol/3), where eps is  approximately  the  square
c  root  of  the  relative  machine  precision.   if   f   is a unimodal
c  function and the computed values of   f   are  always  unimodal  when
c  separated  by  at least  eps*abs(x)+(tol/3), then  fmin  approximates
c  the abcissa of the global minimum of  f  on the interval  ax,bx  with
c  an error less than  3*eps*abs(fmin)+tol.  if   f   is  not  unimodal,
c  then fmin may approximate a local, but perhaps non-global, minimum to
c  the same accuracy.
c      this function subprogram is a slightly modified  version  of  the
c  algol  60 procedure  localmin  given in richard brent, algorithms for
c  minimization without derivatives, prentice-hall, inc. (1973).
c
c
      double precision function fmin(ax,bx,f,tol)

      implicit none
      double precision ax,bx,f,tol
c
      EXTERNAL  d1mach
      double precision  dabs,dsqrt,d1mach

      double precision  a,b,c,d,e,eps,xm,p,q,r,tol1,t2,u,v,w,fu,fv,fw,
     2    fx,x,tol3

c
c  c is the squared inverse of the golden ratio
      c=0.5d0*(3.0d0-dsqrt(5.0d0))
c
c  eps is approximately the square root of the relative machine
c  precision.
c
      eps=d1mach(4)
      tol1=eps+1.0d0
      eps=dsqrt(eps)
c
      a=ax
      b=bx
      v=a+c*(b-a)
      w=v
      x=v
c
c  -Wall indicates that d may be used before being assigned
c
      d=0.0d0
      e=0.0d0
      fx=f(x)
      fv=fx
      fw=fx
      tol3=tol/3.0d0
c
c  main loop starts here -----------------------------------
c
   20 xm=0.5d0*(a+b)
      tol1=eps*dabs(x)+tol3
      t2=2.0d0*tol1
c
c  check stopping criterion
c
      if (dabs(x-xm).le.(t2-0.5d0*(b-a))) go to 900
      p=0.0d0
      q=0.0d0
      r=0.0d0
      if (dabs(e).gt.tol1) then
c     
c     fit parabola
c     
         r=(x-w)*(fx-fv)
         q=(x-v)*(fx-fw)
         p=(x-v)*q-(x-w)*r
         q=2.0d0*(q-r)
         if (q.gt.0.0d0) then
            p=-p
         else
            q=-q
         endif
         r=e
         e=d
      endif

      if ((dabs(p).ge.dabs(0.5d0*q*r)).or.(p.le.q*(a-x))
     2                                .or.(p.ge.q*(b-x))) then

c
c     a golden-section step
c     
         if (x.lt.xm) then
            e=b-x
         else
            e=a-x
         endif
         d=c*e
         
      else
c
c     a parabolic-interpolation step
c     
         d=p/q
         u=x+d
c     
c     f must not be evaluated too close to ax or bx
c     
         if (((u-a).ge.t2).and.((b-u).ge.t2)) go to 90
         d=tol1
         if (x.ge.xm) d=-d
      endif
         
c     
c f must not be evaluated too close to x
c     
 90   if (dabs(d).ge.tol1) then
         u=x+d
      else
         if (d.gt.0.0d0) then
            u=x+tol1
         else
            u=x-tol1
         endif
      endif
      fu=f(u)

c
c  update  a, b, v, w, and x
c
      if (fx.le.fu) then
         if (u.lt.x) then
            a=u
         else
            b=u
         endif
      endif

      if (fu.le.fx) then
         
         if (u.lt.x) then
            b=x
         else
            a=x
         endif
         v=w
         fv=fw
         w=x
         fw=fx
         x=u
         fx=fu
         
      else

 170     if ((fu.le.fw).or.(w.eq.x)) then
            v=w
            fv=fw
            w=u
            fw=fu
         else
            if ((fu.le.fv).or.(v.eq.x).or.(v.eq.w)) then
               v=u
               fv=fu
            endif
         endif
      endif
      
      go to 20
c
c  end of main loop
c
  900 fmin=x
      return
      end