diff --git a/src/fmin.f b/src/fmin.f
new file mode 100644
index 0000000..078a9ee
--- /dev/null
+++ b/src/fmin.f
@@ -0,0 +1,159 @@
+      double precision function fmin(ax,bx,f,tol)
+      double precision ax,bx,f,tol
+c
+c      an approximation  x  to the point where  f  attains a minimum  on
+c  the interval  (ax,bx)  is determined.
+c
+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.0d0)
+c
+c
+c  output..
+c
+c  fmin  abcissa approximating the point where  f  attains a minimum
+c
+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  a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w
+      double precision  fu,fv,fw,fx,x
+      double precision  dabs,dsqrt,dsign
+c
+c  c is the squared inverse of the golden ratio
+c
+      c = 0.5d0*(3. - dsqrt(5.0d0))
+c
+c  eps is approximately the square root of the relative machine
+c  precision.
+c
+      eps = 1.0d00
+   10 eps = eps/2.0d00
+      tol1 = 1.0d0 + eps
+      if (tol1 .gt. 1.0d00) go to 10
+      eps = dsqrt(eps)
+c
+c  initialization
+c
+      a = ax
+      b = bx
+      v = a + c*(b - a)
+      w = v
+      x = v
+      e = 0.0d0
+      fx = f(x)
+      fv = fx
+      fw = fx
+c
+c  main loop starts here
+c
+   20 xm = 0.5d0*(a + b)
+      tol1 = eps*dabs(x) + tol/3.0d0
+      tol2 = 2.0d0*tol1
+c
+c  check stopping criterion
+c
+      if (dabs(x - xm) .le. (tol2 - 0.5d0*(b - a))) go to 90
+c
+c is golden-section necessary
+c
+      if (dabs(e) .le. tol1) go to 40
+c
+c  fit parabola
+c
+      r = (x - w)*(fx - fv)
+      q = (x - v)*(fx - fw)
+      p = (x - v)*q - (x - w)*r
+      q = 2.0d00*(q - r)
+      if (q .gt. 0.0d0) p = -p
+      q =  dabs(q)
+      r = e
+      e = d
+c
+c  is parabola acceptable
+c
+   30 if (dabs(p) .ge. dabs(0.5d0*q*r)) go to 40
+      if (p .le. q*(a - x)) go to 40
+      if (p .ge. q*(b - x)) go to 40
+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) .lt. tol2) d = dsign(tol1, xm - x)
+      if ((b - u) .lt. tol2) d = dsign(tol1, xm - x)
+      go to 50
+c
+c  a golden-section step
+c
+   40 if (x .ge. xm) e = a - x
+      if (x .lt. xm) e = b - x
+      d = c*e
+c
+c  f must not be evaluated too close to x
+c
+   50 if (dabs(d) .ge. tol1) u = x + d
+      if (dabs(d) .lt. tol1) u = x + dsign(tol1, d)
+      fu = f(u)
+c
+c  update  a, b, v, w, and x
+c
+      if (fu .gt. fx) go to 60
+      if (u .ge. x) a = x
+      if (u .lt. x) b = x
+      v = w
+      fv = fw
+      w = x
+      fw = fx
+      x = u
+      fx = fu
+      go to 20
+   60 if (u .lt. x) a = u
+      if (u .ge. x) b = u
+      if (fu .le. fw) go to 70
+      if (w .eq. x) go to 70
+      if (fu .le. fv) go to 80
+      if (v .eq. x) go to 80
+      if (v .eq. w) go to 80
+      go to 20
+   70 v = w
+      fv = fw
+      w = u
+      fw = fu
+      go to 20
+   80 v = u
+      fv = fu
+      go to 20
+c
+c  end of main loop
+c
+   90 fmin = x
+      return
+      end