      function fi1(f,y,x1,x2) ! bernd berg, sep 19 2000.
c inverse of the function f.
c result:     fi1=x such that y=f(x). 
c precission: eps=1/10**8 can be changed, see below.
c method:     bisecting interval, initial range [x1,x2] (x1<x2).
c assumption: f(x) monoton (in- or decreasing) in the interval [x1,x2].
      include '../../libs/fortran/implicit.sta'
      include '../../libs/fortran/constants.par'
      external f
      eps=one/ten**8
      itermax=1000 ! maximum number of iterations.
      iter=0
c
      xx1=x1
      xx2=x2
1     y1=f(xx1)
      y2=f(xx2)
c
      if(y1.gt.y2) then
      xx1=x2
      xx2=x1
      go to 1
      end if
c
      if(y.lt.y1 .or. y.gt.y2) go to 2 
c
3     continue
      iter=iter+1
      if(iter.gt.itermax) go to 2
        xx=half*(xx1+xx2)
        ff=f(xx)
        if(ff.le.y) xx1=xx
        if(ff.gt.y) xx2=xx
      if(abs(xx2-xx1).gt.eps ) go to 3
c
      fi1=half*(xx2+xx1) 
c
      return
2     print*,'fi1: no convergence or y out of range !'
      print*,'itermax,eps:',itermax,eps
      print*,'xx1,xx2:    ',xx1,xx2
      print*,'y1,y,y2:    ',y1,y,y2
      stop 'fi1'
      end
