      function fmin(f,a,b,xmin,iter) ! bernd berg, dec 11 2007.
c calculates minimum of an external function f.
c result:     fmin = f(xmin).
c a,b:        a<b intial range needed on input. the function f has to 
c             be monotonically increasing on both sides of its minimum.
c method:     bisection.
      include '../../libs/fortran/implicit.sta'
      include '../../libs/fortran/constants.par'
c precission: eps=1/10**10 can be changed below.
c itermax, maximum number of iterations can be changed below.
      parameter(eps=one/ten**10,itermax=1000)
      dimension fx(3),x(3)
      external f
      ltest=.true.
c initialization:
      iter=0
      x(1)=a
      x(2)=half*(a+b)
      x(3)=b
      fx(1)=f(x(1))
      fx(2)=f(x(2))
      fx(3)=f(x(3))
      print*,"x:",x
      print*,"fx:",fx
      if(fx(2).ge.fx(1).and.fx(2).ge.fx(3)) stop "fmin error 1."
      if(ltest) write(6,'(/,"iter, f, x:",/)')
c
1     iter=iter+1
      if((x(2)-x(1)).gt.(x(3)-x(2))) then
        xtest=half*(x(2)+x(1))
        ftest=f(xtest)
        if(ftest.lt.fx(2)) then
          fx(3)=fx(2)
          x(3)=x(2)
          fx(2)=ftest
          x(2)=xtest
        else
          if(ftest.ge.fx(1)) stop "fmin error 2."
          fx(1)=ftest
          x(1)=xtest
        endif
      else
        xtest=half*(x(3)+x(2))
        ftest=f(xtest)
        if(ftest.lt.fx(2)) then
          fx(1)=fx(2)
          x(1)=x(2)
          fx(2)=ftest
          x(2)=xtest
        else
          if(ftest.ge.fx(3)) stop "fmin error 3."
          fx(3)=ftest
          x(3)=xtest
        endif
      endif
      del=max(abs(fx(1)-fx(2)),abs(fx(3)-fx(2)))
      if(ltest) write(6,'(i6,6g12.3)') iter,fx,x
      if(del.lt.eps) then
        fmin=fx(2)
        xmin=x(2)
        print*,"fmin,xmin:",fmin,xmin
      return
      endif
      if(iter.lt.itermax) goto 1
c
2     print*,"fmin: no convergence in itermax steps!"
      print*,"itermax,eps:",itermax,eps
      print*,"fx:",fx
      print*,"x:",x
      stop 'fmin'
      end
