      subroutine fpminf0(f,x,xlb,xub,ftol)
c
c (replacing old fpminf0, dmc 10 Apr 1997)
c solve vector f(x)=0 using 2d Newton's method
c   f and x are 2-vectors.
c
c input:
c    f -- function subroutine with call interface of form
c         call f(x,ans,jac)
c           real x(2),ans(2),jac(2,2)
c           x - input argument
c           ans - return value
c           jac - 2x2 jacobian, df/dx
c
c    x -- initial guess **modified on output**
c
c    xlb -- lower bounds for x(1) and x(2)
c      use xlb to bind x(2) away from origin singularity!
c    xub -- upper bounds for x(1) and x(2)
c    ftol -- error tolerance.  Done when
c       max(abs(ans(1)),abs(ans(2))).le.ftol
c
c output:
c    x -- the approximate solution for f(x)=0.
c
c special assumptions:
c    x really indexes polar coordinates (chi,psi) with chi the angle-like
c      cordinate and psi the radius-like coordinate.  The "origin" of the
c      system is at x(2)=xlb(2), where all chi variation disappears.
c
c    so...
c    f is periodic in x(1) with one period given in range (xlb(1),xub(1))
c
c------------------------------------------------------
C============
C idecl:  explicitize implicit INTEGER declarations:
      IMPLICIT NONE
      INTEGER niter,nitmax
C============
C idecl:  explicitize implicit REAL declarations:
      REAL*8 zdamp,xrange,period,zdela,zdmax,zdet,zfac
C============
      external f			! the function subroutine
      REAL*8 x(2)			! init guess & answer returned
      REAL*8 xlb(2)			! x lower bounds
      REAL*8 xub(2)			! x upper bounds
      REAL*8 ftol			! error tolerance
c
      REAL*8 xold(2)                    ! local:  new value of x
      REAL*8 xdel(2)                    ! local:  change vector
c
      REAL*8 ans(2)			! f return value
      REAL*8 jac(2,2)			! df/dx (2x2 jacobian)
c------------------------------------------------------
c
      niter=0
      nitmax=100
      zdamp=1.0			! step damper
c
      xrange=xub(2)-xlb(2)
c
      period=xub(1)-xlb(1)
      zdela=0.45d0*period               ! angle change if origin is crossed
c
c  factor of 0.45 is used instead of 0.5 to prevent stable oscillation
c  across the origin
c
 100  continue
      niter=niter+1
c
      if(niter.eq.nitmax/2) then
         write(6,
     >'('' %fpminf0 warning:  high iteration count'')')
         zdamp=0.6666667d0
      endif
c
      if(niter.gt.nitmax) then
         call errmsg_exit('?too many iterations in fpminf0')
      endif
c
c  evaluate function
c
      call f(x,ans,jac)
c
c  tolerance test
c
      zdmax=max(abs(ans(1)),abs(ans(2)))
      if(zdmax.lt.ftol) go to 900	! success
c
c  ok solve for next guess using Newton method for change.
c  the coordinate system is assumed to be well behaved!
c
      xold(1)=x(1)
      xold(2)=x(2)
c
c  direct 2x2 solution  of  jac*xdel = -ans
c
      zdet=jac(1,1)*jac(2,2)-jac(2,1)*jac(1,2)
c
      xdel(1)=-(jac(2,2)*ans(1)-jac(1,2)*ans(2))/zdet
      xdel(2)=-(jac(1,1)*ans(2)-jac(2,1)*ans(1))/zdet
c
c  new x
c
      x(1)=xold(1)+zdamp*xdel(1)
      x(2)=xold(2)+zdamp*xdel(2)
c
c  constraints on one step change
c
      zfac=0.25d0
      x(1)=max((xold(1)-zfac*period),min((xold(1)+zfac*period),x(1)))
      x(2)=max((xold(2)-zfac*xrange),min((xold(2)+zfac*xrange),x(2)))
c
c  global range checks
c
      if(x(2).lt.xlb(2)) then
c  crossed the origin
         x(2)=xlb(2)+(xlb(2)-x(2))	! cross back
         x(1)=x(1)+zdela              ! correct the angle
      endif
      x(2)=min(x(2),xub(2))
c
c  branch cut check
c
      if(x(1).lt.xlb(1)) then
         x(1)=x(1)+period
      else if(x(1).gt.xub(1)) then
         x(1)=x(1)-period
      endif
c
c  next iteration
c
      go to 100
c
c-------------------------------
c
 900  continue
      return
      end


! 09May2001 fgtok -s r8_precision.sub xpl_r4_remove.inf r8_convert.list "r8con.csh conversion"
! 09May2001 fgtok
