C	subroutine "dispcurve" to compute dispersion curve using grderiv

C	Sven Ivansson 1993






	SUBROUTINE DISPCURVE ( HSTEPFSUBPRE, NFR,CFREQARR, USTART, 
     $	DSFMIN,DSFMAX,  DSFFAKTMAX, DSUMAX, DRELGOALFAKT, RSUG2, 
     $	DLOGREDUCFAKT,DDELT,XYTOL, NFROK,UCURVE,DUDFCURVE,UGRCURVE, 
     $	IGRCALL,IASFAIL, TERM)
C	hstepfsubpre         for hstepfsub in "dpcfsub.inc"
C	nfr                  number of complex frequencies
C	cfreqarr(1:nfr)      complex frequencies 
C	                         ( cfreqarr(1) is the start frequency )
C	ustart               mode-u at cfreqarr(1), it must not be on a 
C	                         branch-cut
C	dsfmin               minimal absolute freq step that is handled 
C	                     (strict limit), always take  dsfmin >= 0.0d0
C	dsfmax               maximal absolute freq step
C	dsffaktmax           maximal factor for change of "planned" freq steps
C	dsumax               maximum allowed absolute "planned" change of u 
C	                     in each freq step
C	drelgoalfakt         for choosing "next freq step", always take  
C	                     drelgoalfakt < 1.0d0
C	rsug2                half-width for u tolerance rectangle to dzeroholo
C	dlogreducfakt        "dlog" of required reduction factor in each step 
C	                     for the absolute value of "the function"  
C	                     ("dzeroholo")
C	ddelt                distance-difference limit for "secant derivative"
C	                     ("dzeroholo")
C	xytol                desired accuracy for the zeros  ("dzeroholo")
C	nfrok               output number of complex frequencies ok
C	ucurve(1:nfrok)     output  u values with ucurve(1):=ustart
C	dudfcurve(1:nfrok)  output  du/df values
C	ugrcurve(1:nfrok)   output  group-u values 
C	igrcall             input/output (OBS)  performed number of function 
C	                                                         evaluations 
C	iasfail             output  0 upon successful result
C	                            1 for inappropriate cfreqarr
C	                            2 for zero "deriv-du" encountered
C	                            3 for "dsf<=dsfmin"
C	term                input .true. for successive ucurve output on 
C	                    terminal

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C	nol <= noldim

	INCLUDE 'INMODEL.INC'
	INCLUDE 'INPARAM.INC'
C	note: "icuttyp,isheettyp" etc must be set from calling program
	INCLUDE 'DPCFSUB.INC'

	INTEGER NFR
	REAL*8 HSTEPFSUBPRE, DSFMIN,DSFMAX, 
     $	DSFFAKTMAX, DSUMAX, DRELGOALFAKT, RSUG2, 
     $	DLOGREDUCFAKT,DDELT,XYTOL
	COMPLEX*16 USTART, 
     $	CFREQARR(*)
	LOGICAL TERM

	INTEGER IGRCALL

	INTEGER NFROK, IASFAIL
	COMPLEX*16 UCURVE(*), DUDFCURVE(*), UGRCURVE(*)

	INTEGER NXTFR, IDZHFAIL
	REAL*8 RLNDET,DLNDO,DLNDU, SFTOT, SF,SFOLD, DSF,DSFOLD, 
     $	RSUG2NOW, PX,PY, X1,X2,Y1,Y2, XYTOLUT, PI, PI2,PI2INV
	COMPLEX*16 CFREQRIKT, CFREQ,CFREQOLD, UP,UPOLD, 
     $	CDET,DERIVDO,DERIVDU, DUDO
	LOGICAL CUTTING

	EXTERNAL GRVALUESUB

	DATA PI /3.141592653589793D0/

	PI2 = 2*PI
	PI2INV = 1/PI2

	IF (U0.NE.(1.0D0,0.0D0)) STOP 3415   ! u,up must be interchangeable; 
C	branch-points and -cuts must be frequency-independent

	IF ( DSFMIN.LT.0.0D0 .OR. DRELGOALFAKT.GE.1.0D0 ) STOP 27  ! safety

	CFREQ = CFREQARR(1)
	OMEG = PI2 * CFREQ
	HSTEPFSUB = HSTEPFSUBPRE / CDABS(CFREQ)      ! obs
	RK0 = OMEG                 ! note that "u0=(1.0d0,0.0d0)"
	OMEG2 = OMEG*OMEG

	IGRCALL = IGRCALL + 1
	CALL GRDERIV ( USTART, RLNDET,CDET, DLNDO,DERIVDO, 
     $	DLNDU,DERIVDU )                         ! note: u,up interchangeable
	IF (DERIVDU.EQ.(0.0D0,0.0D0)) THEN
	  IASFAIL = 2        ! zero "deriv-u" encountered
	  NFROK = 0
	  RETURN
	END IF
	DUDO = - DEXP(DLNDO-DLNDU) * DERIVDO/DERIVDU
	UCURVE(1) = USTART
	DUDFCURVE(1) = PI2*DUDO
	UGRCURVE(1) = USTART + OMEG*DUDO

	NXTFR = 2
	CFREQRIKT = CFREQARR(2) - CFREQ
	SFTOT = CDABS(CFREQRIKT)
	IF (SFTOT.EQ.0.0D0) THEN
	  IASFAIL = 1
	  NFROK = 1
	  RETURN
	END IF
	CFREQRIKT = CFREQRIKT/SFTOT
	SF = 0.0D0

	DSFOLD = DSFMAX
	UP = USTART                             ! note: u,up interchangeable

	DO WHILE (NXTFR.LE.NFR)

	  CFREQOLD = CFREQ
	  SFOLD = SF
	  UPOLD = UP
	  IF (DUDO.NE.(0.0D0,0.0D0)) THEN
	    DSF = DRELGOALFAKT * DMIN1 ( DSFMAX , PI2INV*DSUMAX/CDABS(DUDO) )
	  ELSE
	    DSF = DRELGOALFAKT * DSFMAX
	  END IF
	  IF ( DSF .GT. DSFFAKTMAX*DSFOLD ) DSF = DSFFAKTMAX*DSFOLD
 100	  IF (DSF.LE.DSFMIN) THEN
	    IASFAIL = 3
	    NFROK = NXTFR - 1
	    RETURN
	  END IF
	  DSFOLD = DSF
	  SF = SFOLD + DSF
	  IF (SF.LT.SFTOT) THEN
	    CFREQ = CFREQOLD + DSF*CFREQRIKT
	  ELSE
	    SF = SFTOT
	    CFREQ = CFREQARR(NXTFR)
	  END IF

	  OMEG = PI2 * CFREQ
	  HSTEPFSUB = HSTEPFSUBPRE / CDABS(CFREQ)      ! obs
	  RK0 = OMEG                 ! note that "u0=(1.0d0,0.0d0)"
	  OMEG2 = OMEG*OMEG

	  X1 = DREAL(UP)
	  Y1 = DIMAG(UP)
	  UP = UPOLD + PI2*(CFREQ-CFREQOLD)*DUDO   ! note: u,up interchangeable
	  PX = DREAL(UP)
	  PY = DIMAG(UP)
	  IF ( CUTTING(X1,PX,Y1,PY) ) THEN 
	    SF = SFOLD
	    DSF = 0.5D0 * CDABS(CFREQ-CFREQOLD)
	    GOTO 100
	  END IF
	  PX = 0.01D0*X1 + 0.99D0*PX    ! cutting border-safety
	  PY = 0.01D0*Y1 + 0.99D0*PY    ! cutting border-safety
	  UP = DCMPLX(PX,PY)                    ! note: u,up interchangeable
	  IGRCALL = IGRCALL + 1
	  CALL GRVALUE ( UP, RLNDET,CDET )
	  RSUG2NOW = RSUG2
 200	  X1 = PX - RSUG2NOW
	  X2 = PX + RSUG2NOW
	  Y1 = PY - RSUG2NOW
	  Y2 = PY + RSUG2NOW
	  IF ( CUTTING(X1,X2,Y1,Y2) ) THEN 
	    RSUG2NOW = 0.5D0 * RSUG2NOW
	    GOTO 200
	  END IF
	  CALL DZEROHOLO ( GRVALUESUB, X1,X2,Y1,Y2, DLOGREDUCFAKT,
     $	  DDELT,XYTOL, PX,PY,XYTOLUT,RLNDET,CDET, IGRCALL,IDZHFAIL )
	  IF (IDZHFAIL.NE.0) THEN
	    SF = SFOLD
	    DSF = 0.5D0 * CDABS(CFREQ-CFREQOLD)
	    GOTO 100
	  END IF

	  UP = DCMPLX(PX,PY)                ! note: u,up interchangeable
	  CALL GRDERIVREST ( UP, DLNDO,DERIVDO, DLNDU,DERIVDU )
	  IF (DERIVDU.EQ.(0.0D0,0.0D0)) THEN
	    IASFAIL = 2        ! zero "deriv-u" encountered
	    NFROK = NXTFR - 1
	    RETURN
	  END IF
	  DUDO = - DEXP(DLNDO-DLNDU) * DERIVDO/DERIVDU

	  IF (CFREQ.EQ.CFREQARR(NXTFR)) THEN
	    UCURVE(NXTFR) = UP                  ! note: u,up interchangeable
	    IF (TERM) WRITE(*,*) SNGL(DREAL(CFREQ)),UP
	    DUDFCURVE(NXTFR) = PI2*DUDO         ! note: u,up interchangeable
	    UGRCURVE(NXTFR) = UP + OMEG*DUDO    ! note: u,up interchangeable
	    NXTFR = NXTFR + 1
	    CFREQRIKT = CFREQARR(NXTFR) - CFREQ
	    SFTOT = CDABS(CFREQRIKT)
	    IF (SFTOT.EQ.0.0D0) THEN
	      IASFAIL = 1
	      NFROK = NXTFR - 1
	      RETURN
	    END IF
	    CFREQRIKT = CFREQRIKT / SFTOT
	    SF = 0.0D0
	  END IF

	END DO

	NFROK = NFR

	RETURN
	END




	LOGICAL FUNCTION CUTTING ( X1,X2,Y1,Y2 )
C	to test if the input rectangle intersects the branch-cuts
C	must have  x1 <= x2  and  y1 <= y2

	IMPLICIT NONE

	INCLUDE 'INPARAM.INC'

	REAL*8 X1,X2,Y1,Y2

	INTEGER IGANG, IC, INUTI, ISLASK
	REAL*8 XBP(4), YBP(4), XYHYP(4)
	LOGICAL CUT(4)
	DATA IGANG /0/
	SAVE CUT, XBP,YBP,XYHYP

	IF (IGANG.EQ.0) THEN
	  IF (UPBPA.NE.(0.0D0,0.0D0)) THEN
	    CUT(1) = .TRUE.
	    IF (UPBPB.NE.(0.0D0,0.0D0)) THEN
	      CUT(2) = .TRUE.
	    ELSE
	      CUT(2) = .FALSE.
	    END IF
	  ELSE
	    CUT(1) = .FALSE.
	    CUT(2) = .FALSE.
	  END IF
	  IF (CUT(1)) THEN
	    XBP(1) = DREAL(UPBPA)
	    YBP(1) = DIMAG(UPBPA)
	    XYHYP(1) = XBP(1)*YBP(1)
	  END IF
	  IF (CUT(2)) THEN
	    XBP(2) = DREAL(UPBPB)
	    YBP(2) = DIMAG(UPBPB)
	    XYHYP(2) = XBP(2)*YBP(2)
	  END IF
	  IF (UPBPASURF.NE.(0.0D0,0.0D0)) THEN
	    CUT(3) = .TRUE.
	    IF (UPBPBSURF.NE.(0.0D0,0.0D0)) THEN
	      CUT(4) = .TRUE.
	    ELSE
	      CUT(4) = .FALSE.
	    END IF
	  ELSE
	    CUT(3) = .FALSE.
	    CUT(4) = .FALSE.
	  END IF
	  IF (CUT(3)) THEN
	    XBP(3) = DREAL(UPBPASURF)
	    YBP(3) = DIMAG(UPBPASURF)
	    XYHYP(3) = XBP(3)*YBP(3)
	  END IF
	  IF (CUT(4)) THEN
	    XBP(4) = DREAL(UPBPBSURF)
	    YBP(4) = DIMAG(UPBPBSURF)
	    XYHYP(4) = XBP(4)*YBP(4)
	  END IF
	  IGANG = 1
	END IF

	CUTTING = .FALSE.   ! temporary

	DO IC = 1, 4

	  IF (CUT(IC)) THEN
	    IF (X1.GT.XBP(IC)) THEN
	    ELSE IF (X2.LT.-XBP(IC)) THEN
	    ELSE IF ( -XBP(IC).LT.Y1 .AND. Y2.LT.XBP(IC) ) THEN
	    ELSE
C	      now treating y1
	      IF (DABS(Y1).LT.YBP(IC)) THEN
	        INUTI = 1
	      ELSE
	        IF (ICUTTYP.EQ.1) THEN
C	          now the point "x1,y1"
	          IF ( X1*Y1 .LT. XYHYP(IC) ) THEN
	            INUTI = 1
	          ELSE
	            INUTI = 0
	          END IF
C	          now the point "x2,y1"
	          IF ( X2*Y1 .LT. XYHYP(IC) ) THEN
	            ISLASK = 1
	          ELSE
	            ISLASK = 0
	          END IF
	        ELSE
	          IF (ICUTTYP.NE.2) STOP 45178      ! not implemented
C	          now the point "x1,y1"
	          IF ( X1.LT.0.0D0 .AND. Y1.GT.0.0D0 ) THEN
	            INUTI = 1
	          ELSE IF ( X1.GT.0.0D0 .AND. Y1.LT.0.0D0 ) THEN
	            INUTI = 1
	          ELSE IF ( DABS(X1) .LT. XBP(IC) ) THEN
	            INUTI = 1
	          ELSE
	            INUTI = 0
	          END IF
C	          now the point "x2,y1"
	          IF ( X2.LT.0.0D0 .AND. Y1.GT.0.0D0 ) THEN
	            ISLASK = 1
	          ELSE IF ( X2.GT.0.0D0 .AND. Y1.LT.0.0D0 ) THEN
	            ISLASK = 1
	          ELSE IF ( DABS(X2) .LT. XBP(IC) ) THEN
	            ISLASK = 1
	          ELSE
	            ISLASK = 0
	          END IF
	        END IF
	        IF (ISLASK.NE.INUTI) THEN
	          CUTTING = .TRUE.
	          GOTO 999
	        END IF
	      END IF
C	      now treating y2
	      IF (DABS(Y2).LT.YBP(IC)) THEN
	        IF (INUTI.NE.1) THEN
	          CUTTING = .TRUE.
	          GOTO 999
	        END IF
	      ELSE
	        IF (ICUTTYP.EQ.1) THEN
C	          now the point "x1,y2"
	          IF ( X1*Y2 .LT. XYHYP(IC) ) THEN
	            IF (INUTI.NE.1) THEN
	              CUTTING = .TRUE.
	              GOTO 999
	            END IF
	          ELSE
	            IF (INUTI.NE.0) THEN
	              CUTTING = .TRUE.
	              GOTO 999
	            END IF
	          END IF
C	          now the point "x2,y2"
	          IF ( X2*Y2 .LT. XYHYP(IC) ) THEN
	            IF (INUTI.NE.1) THEN
	              CUTTING = .TRUE.
	              GOTO 999
	            END IF
	          ELSE
	            IF (INUTI.NE.0) THEN
	              CUTTING = .TRUE.
	              GOTO 999
	            END IF
	          END IF
	        ELSE
	          IF (ICUTTYP.NE.2) STOP 45178      ! not implemented
C	          now the point "x1,y2"
	          IF ( X1.LT.0.0D0 .AND. Y2.GT.0.0D0 ) THEN
	            IF (INUTI.NE.1) THEN
	              CUTTING = .TRUE.
	              GOTO 999
	            END IF
	          ELSE IF ( X1.GT.0.0D0 .AND. Y2.LT.0.0D0 ) THEN
	            IF (INUTI.NE.1) THEN
	              CUTTING = .TRUE.
	              GOTO 999
	            END IF
	          ELSE IF ( DABS(X1) .LT. XBP(IC) ) THEN
	            IF (INUTI.NE.1) THEN
	              CUTTING = .TRUE.
	              GOTO 999
	            END IF
	          ELSE
	            IF (INUTI.NE.0) THEN
	              CUTTING = .TRUE.
	              GOTO 999
	            END IF
	          END IF
C	          now the point "x2,y2"
	          IF ( X2.LT.0.0D0 .AND. Y2.GT.0.0D0 ) THEN
	            IF (INUTI.NE.1) THEN
	              CUTTING = .TRUE.
	              GOTO 999
	            END IF
	          ELSE IF ( X2.GT.0.0D0 .AND. Y2.LT.0.0D0 ) THEN
	            IF (INUTI.NE.1) THEN
	              CUTTING = .TRUE.
	              GOTO 999
	            END IF
	          ELSE IF ( DABS(X2) .LT. XBP(IC) ) THEN
	            IF (INUTI.NE.1) THEN
	              CUTTING = .TRUE.
	              GOTO 999
	            END IF
	          ELSE
	            IF (INUTI.NE.0) THEN
	              CUTTING = .TRUE.
	              GOTO 999
	            END IF
	          END IF
	        END IF
	      END IF
	    END IF
	  END IF

	END DO

 999	RETURN
	END