C	subroutine to refine a "zero within a rectangle" of an 
C       analytic function

C	sven ivansson   1990




	SUBROUTINE DZEROHOLO (CDETSUB, X1,X2,Y1,Y2, DLOGREDUCFAKT,DDELT,
     $	XYTOL, PX,PY,XYTOLUT,RLNDET,CDET, IDZHFCALL,IDZHFAIL)
C	cdetsub          subroutine cdetsub (px,py, rlndet,cdet)  gives 
C	                 "the function" as  dexp(rlndet)*cdet  , "the 
C	                 function" must be analytic inside the rectangle 
C	                 and continuous up to the boundary and zero-free 
C	                 on the boundary  ,  up = px + i*py
C	x1,x2            horizontal limits of the rectangle
C	y1,y2            vertical limits of the rectangle
C	dlogreducfakt    "dlog" of required reduction factor in each step for 
C	                 the absolute value of "the function"
C	ddelt            distance-difference limit for "secant derivative"
C	xytol            desired accuracy for the zero (close to the borders 
C	                 of the rectangle it will be better, see "xytolut")
C	px,py          input/output (OBS)  real & imaginary part of the zero
C	                         should be "within the rectangle" on entry 
C	                         (otherwise the "midpoint" is used)
C	xytolut        output  estimated accuracy for the zero so that 
C	                 * xytolut is an upper bound for the Euclidean 
C	                   distance from (px,py) to the zero
C	                 * the "rectangle (px+-xytolut;py+-xytolut)" is 
C	                   within the "rectangle (x1,x2;y1,y2)"
C	rlndet,cdet    input/output (OBS)  give "function value" at (px,py)
C	idzhfcall      input/output (OBS)  total number of cdetsub calls
C	idzhfail       output  "return status":  0 for successful result
C	                         1 for "unsuccessful, try a smaller rectangle"
C	                         2 for "zero found on the boundary"

	IMPLICIT NONE

	REAL*8 X1,X2,Y1,Y2, DLOGREDUCFAKT,DDELT,XYTOL

	INTEGER IDZHFCALL
	REAL*8 PX,PY, RLNDET
	COMPLEX*16 CDET

	INTEGER IDZHFAIL
	REAL*8 XYTOLUT

	REAL*8 XSIDG2,YSIDG2, DDELTLOK, PXD,PYD, PXNEW,PYNEW, RLNDETD,
     $	RLNDETNEW, XDIFF,YDIFF, SLASK
	COMPLEX*16 CDETD,CDETNEW, ZDIFFD, ZDIFF, CSLASK

        EXTERNAL CDETSUB

	XSIDG2 = 0.5D0 * (X2-X1)
	YSIDG2 = 0.5D0 * (Y2-Y1)
	DDELTLOK = DMIN1 ( 0.5D0*DMIN1(XSIDG2,YSIDG2) , DDELT )

	IF (PX.LE.X1 .OR. X2.LE.PX .OR. PY.LE.Y1 .OR. Y2.LE.PY) THEN
	  PX = 0.5D0 * (X1+X2)
	  PY = 0.5D0 * (Y1+Y2)
	  IDZHFCALL = IDZHFCALL + 1
	  CALL CDETSUB (PX,PY, RLNDET,CDET)
	END IF
	IF (CDET.EQ.(0.0D0,0.0D0)) THEN
	  XYTOLUT = 0.0D0
	  IDZHFAIL = 0
	  RETURN
	END IF
	PXD = PX + DDELTLOK
	IF (PXD.GT.X2) PXD = PX - DDELTLOK
	PYD = PY
	ZDIFFD = DCMPLX(PX-PXD,0.0D0)

 100	IDZHFCALL = IDZHFCALL + 1
	CALL CDETSUB (PXD,PYD, RLNDETD,CDETD)
	IF (CDETD.EQ.(0.0D0,0.0D0)) THEN
	  PX = PXD
	  PY = PYD
	  XYTOLUT = 0.0D0
	  RLNDET = RLNDETD
	  CDET = CDETD
	  IF (X1.LT.PX .AND. PX.LT.X2 .AND. Y1.LT.PY .AND. PY.LT.Y2) THEN 
	    IDZHFAIL = 0
	  ELSE
	    IDZHFAIL = 2
	  END IF
	  RETURN
	END IF

 200	SLASK = RLNDETD - RLNDET
	IF (DABS(SLASK).GE.88.0D0) THEN     !  1.7*(10**38) "=" exp(88.2)
	  IDZHFAIL = 1
	  RETURN
	END IF
	CSLASK = 1 - DEXP(SLASK)*(CDETD/CDET)
	IF (CSLASK.EQ.(0.0D0,0.0D0)) THEN
	  IDZHFAIL = 1
	  RETURN
	END IF
	ZDIFF = - ZDIFFD / CSLASK          ! in general the "secant rule"

	XDIFF = DREAL(ZDIFF)
	YDIFF = DIMAG(ZDIFF)
	PXNEW = PX + XDIFF
	PYNEW = PY + YDIFF
	IF (PXNEW.LT.X1) THEN
	  PYNEW = PY + (X1-PX)*(YDIFF/XDIFF)
	  PXNEW = X1
	ELSE IF (PXNEW.GT.X2) THEN
	  PYNEW = PY + (X2-PX)*(YDIFF/XDIFF)
	  PXNEW = X2
	END IF
	IF (PYNEW.LT.Y1) THEN
	  PXNEW = PX + (Y1-PY)*(XDIFF/YDIFF)
	  PYNEW = Y1
	ELSE IF (PYNEW.GT.Y2) THEN
	  PXNEW = PX + (Y2-PY)*(XDIFF/YDIFF)
	  PXNEW = Y2
	END IF
	IF (PXNEW.LT.X1) PXNEW = X1
	IF (PXNEW.GT.X2) PXNEW = X2
	IF (PYNEW.LT.Y1) PYNEW = Y1
	IF (PYNEW.GT.Y2) PYNEW = Y2
	ZDIFF = DCMPLX(PXNEW-PX,PYNEW-PY)
	IF (ZDIFF.EQ.(0.0D0,0.0D0)) THEN
	  IDZHFAIL = 1           ! this is extremely cautious
	  RETURN
	END IF

	IDZHFCALL = IDZHFCALL + 1
	CALL CDETSUB (PXNEW,PYNEW, RLNDETNEW,CDETNEW)
	IF (CDETNEW.EQ.(0.0D0,0.0D0)) THEN
	  PX = PXNEW
	  PY = PYNEW
	  XYTOLUT = 0.0D0
	  RLNDET = RLNDETNEW
	  CDET = CDETNEW
	  IF (X1.LT.PX .AND. PX.LT.X2 .AND. Y1.LT.PY .AND. PY.LT.Y2) THEN 
	    IDZHFAIL = 0
	  ELSE
	    IDZHFAIL = 2
	  END IF
	  RETURN
	END IF

	SLASK = RLNDETNEW - RLNDET + DLOG(CDABS(CDETNEW/CDET))
	IF (SLASK.GE.DLOGREDUCFAKT) THEN
	  IDZHFAIL = 1
	  RETURN
	END IF

	SLASK = 1.414213562373095D0 * DMAX1 ( DABS(DREAL(ZDIFF)), 
     $	DABS(DIMAG(ZDIFF)) )   ! "cheap" upper bound for Euclidean norm
	IF (SLASK.GE.XYTOL) THEN
	ELSE IF ( (PXNEW-SLASK) .LE. X1 ) THEN
	ELSE IF ( (PXNEW+SLASK) .GE. X2 ) THEN
	ELSE IF ( (PYNEW-SLASK) .LE. Y1 ) THEN
	ELSE IF ( (PYNEW+SLASK) .GE. Y2 ) THEN
	ELSE
	  PX = PXNEW
	  PY = PYNEW
	  XYTOLUT = SLASK
	  RLNDET = RLNDETNEW
	  CDET = CDETNEW
	  IDZHFAIL = 0
	  RETURN
	END IF
	IF (SLASK.LE.DDELTLOK) THEN
	  PXD = PX
	  PYD = PY
	  RLNDETD = RLNDET
	  CDETD = CDET
	  PX = PXNEW
	  PY = PYNEW
	  RLNDET = RLNDETNEW
	  CDET = CDETNEW
	  ZDIFFD = DCMPLX(PX-PXD,PY-PYD)
	  GOTO 200
	ELSE
	  PXD = PXNEW + DDELTLOK
	  IF (PXD.GT.X2) PXD = PXNEW - DDELTLOK
	  PYD = PYNEW
	  PX = PXNEW
	  PY = PYNEW
	  RLNDET = RLNDETNEW
	  CDET = CDETNEW
	  ZDIFFD = DCMPLX(PX-PXD,0.0D0)
	  GOTO 100
	END IF

	END