
C	subroutine for rpress.f, branch-cut integral(s)

C	Sven Ivansson 1990




	SUBROUTINE RPBRANCHINT ( FSOURC,PM0GOR, IUROT, IFILONPRE, 
     $	NPU,AIMSTA,AIMEND, NHYDLI, 
     $	ICWORKLIM0,ICWORKDIM,CPRESSW )
C	for the "up-hyperbolic" branch-cut path(s)
C	   npu =  -1   --->   P-cut integral
C	   npu =  -2   --->   S-cut integral
C	   npu = -12   --->   P-cut and S-cut integrals
C	hyperbola parameterization with   t = dimag(up) - dreal(up)
C	       tps(1) = t(upbpa)    hps(1) = dreal(upbpa)**2 - dimag(upbpa)**2
C	       tps(2) = t(upbpb)    hps(2) = dreal(upbpb)**2 - dimag(upbpb)**2
C	  P-cut equation   dreal(up) * dimag(up) = psprod(1)  ,  t >= tps(1)
C	     let    slaskrot = dsqrt ( t*t + 4*psprod(1) )
C	            yslask = 0.5d0 * ( t + slaskrot )
C	     then   up = dcmplx ( psprod(1)/yslask , yslask )
C	        dup/dt = (-1/slaskrot) * dconjg(up)
C	  S-cut equation   dreal(up) * dimag(up) = psprod(2)  ,  t >= tps(2)
C	     let    slaskrot = dsqrt ( t*t + 4*psprod(2) )
C	            yslask = 0.5d0 * ( t + slaskrot )
C	     then   up = dcmplx ( psprod(2)/yslask , yslask )
C	        dup/dt = (-1/slaskrot) * dconjg(up)
C	hyperbola parameterization with   s = cdsqrt ( upbp**2 - up**2 ) = 
C                               dsqrt ( hps - dreal(up)**2 + dimag(up)**2 ) = 
C	                        dsqrt ( hps + t * dsqrt(t*t+4*psprod) )
C	  P-cut equation   dreal(up) * dimag(up) = psprod(1)  ,  s >= 0.0d0
C	     let    sq = hps(1) - s*s
C	            slaskrot = dsqrt ( sq*sq + 4*psprod(1)*psprod(1) )
C	            xslask = dsqrt ( 0.5d0 * ( sq + slaskrot ) )
C	     then   up = dcmplx ( xslask , psprod(1)/xslask )
C	        dup/ds = (-s/slaskrot) * dconjg(up)
C	  S-cut equation   dreal(up) * dimag(up) = psprod(2)  ,  s >= 0.0d0
C	     let    sq = hps(2) - s*s
C	            slaskrot = dsqrt ( sq*sq + 4*psprod(2)*psprod(2) )
C	            xslask = dsqrt ( 0.5d0 * ( sq + slaskrot ) )
C	     then   up = dcmplx ( xslask , psprod(2)/xslask )
C	        dup/ds = (-s/slaskrot) * dconjg(up)
C	modelled after rpressint

	USE INMODELa
	USE INRANGEa
	USE INTWORKa

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C       nol <= noldim , nhyd <= nhyddim , li <= lidim

	INCLUDE 'INPARAM.INC'
	INCLUDE 'INPINT.INC'
	INCLUDE 'INPINTX.INC'
	INCLUDE 'INTWORK.INC'
	INCLUDE 'BRANCHPS.INC'

	INTEGER IUROT, IFILONPRE, NPU, NHYDLI
	REAL*8 PM0GOR(*), AIMSTA,AIMEND
	COMPLEX*16 FSOURC

	INTEGER ICWORKLIM0,ICWORKDIM
	COMPLEX*16 CPRESSW(*)

	INTEGER IR, IFILONLOK, NKVAR, NFCAL
	REAL*8 TPSCORN, TSLASK, EPSMAGNLOK, PI, 
     $	TSTAPS(2),TENDPS(2)
	COMPLEX*16 CSLASK, 
     $	GAMMA(lidim), TZ(2), SZ(2)
	LOGICAL DOPS(2)
	EXTERNAL DHBRANCHT, RBRANCHT,RBRANCHPOST , 
     $	         DHBRANCHS, RBRANCHS,RBRANCHPOSS

	DATA PI /3.141592653589793D0/

C	---------------------------------------------------------------------

	IF (NPU.NE.-2) THEN
	  TPS(1) = DIMAG(UPBPA) - DREAL(UPBPA)
	  IF (TPS(1).GE.0.0D0) STOP 560           ! this case not impl here yet
	  HPS(1) = -TPS(1) * ( DIMAG(UPBPA) + DREAL(UPBPA) )
	  DOPS(1) = .TRUE.
	ELSE
	  DOPS(1) = .FALSE.
	END IF
	IF (NPU.NE.-1) THEN
	  TPS(2) = DIMAG(UPBPB) - DREAL(UPBPB)
	  IF (TPS(2).GE.0.0D0) STOP 560           ! this case not impl here yet
	  HPS(2) = -TPS(2) * ( DIMAG(UPBPB) + DREAL(UPBPB) )
	  DOPS(2) = .TRUE.
	ELSE
	  DOPS(2) = .FALSE.
	END IF

	IF (RHO(NOL).GT.0.0D0) THEN
	  PSPROD(1) = DREAL(UPBPA) * DIMAG(UPBPA)
	  IF (PSPROD(1).LE.0.0D0) STOP 560 ! this case not implemented here yet
	  IF (DOPS(1)) THEN
	    IF (AIMSTA.GT.DREAL(UPBPA)) THEN
	      TSTAPS(1) = AIMSTA - PSPROD(1)/AIMSTA
	    ELSE
	      TSTAPS(1) = TPS(1)
	    END IF
	    IF (AIMEND.GT.DREAL(UPBPA)) THEN
	      TENDPS(1) = AIMEND - PSPROD(1)/AIMEND
	    ELSE
	      TENDPS(1) = TPS(1)
	    END IF
	  END IF
	  IF (ISOLID(NOL).NE.0) THEN
	    PSPROD(2) = DREAL(UPBPB) * DIMAG(UPBPB)
	    IF (PSPROD(2).LE.PSPROD(1)) STOP 560  ! this case not impl here yet
	    IF (DOPS(2)) THEN
	      IF (AIMSTA.GT.DREAL(UPBPA)) THEN
	        TSTAPS(2) = AIMSTA - PSPROD(2)/AIMSTA
	      ELSE
	        TSTAPS(2) = TPS(2)
	      END IF
	      IF (AIMEND.GT.DREAL(UPBPA)) THEN
	        TENDPS(2) = AIMEND - PSPROD(2)/AIMEND
	      ELSE
	        TENDPS(2) = TPS(2)
	      END IF
	    END IF
	  ELSE
	    IF (NPU.NE.-1) STOP 560
	  END IF
	ELSE
	  STOP 560
	END IF

	DO IPS = 1, 2

	  IF (DOPS(IPS)) THEN

	    IF ( IFILONPRE.NE.0 .AND. IFILONPRE.NE.-3 ) THEN
C	      adapt to large t
	      DO IR = 1, LI
	        GAMMA(IR) = -RK0*OR(IR)
	      END DO
	    END IF

	    TPSCORN = 0.5D0 * TPS(IPS)         ! t-parameter at "change t/s"
	    IF (TPSCORN.LE.TPS(IPS)) STOP 770

	    TZ(1) = DCMPLX ( DMAX1(TSTAPS(IPS),TPSCORN) , 0.0D0 )
	    TZ(2) = DCMPLX ( DMAX1(TENDPS(IPS),TPSCORN) , 0.0D0 )
	    IF (IADINTXFI.EQ.0) THEN
	      CALL TRAPEZFI (NHYD,LI, DHBRANCHT,RBRANCHT, 
     $	      IFILON,RBRANCHPOST,GAMMA, 2,TZ, HMAX, NHYDLI,CPRESSW, 
     $	      NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,CPRESSW)
	      NFCALLS(3) = NFCALLS(3) + NFCAL
	    ELSE
	      CALL ADINTXFI (NHYD,LI, DHBRANCHT,RBRANCHT, 
     $	      IFILON,RBRANCHPOST,GAMMA, 2,TZ, HMAX,SMAXPRE,EPSH01Q,EPS,
     $	      NELMAXPRE,NBACKMAX, IXTYP,NRUS,IRUS,JRUS, ARGDIFLIM1,
     $	      ARGDIFLIM0,ARGDIFLIM2, 0.0D0,0,NKVAR, IADINTXFI,TEXTLOG, 
     $	      NHYDLI,CPRESSW, EPSMAGNLOK, NFCAL, DHFAKT,RFAKT, 
     $	      ICWORKLIM0,ICWORKDIM,CPRESSW)
	      NFCALLS(3) = NFCALLS(3) + NFCAL
	      IF (EPSMAGNLOK.GT.SUPEPSMAGN) SUPEPSMAGN = EPSMAGNLOK
	    END IF

	    IF ( IFILONPRE.NE.0 .AND. IFILONPRE.NE.-3 ) THEN
C	      adapt to small s
	      IFILONLOK = 0
	    ELSE
	      IFILONLOK = IFILONPRE
	    END IF

	    TSLASK = DMIN1 ( DMAX1(TSTAPS(IPS),TPS(IPS)) , TPSCORN )
	    IF (TSLASK.NE.TPS(IPS)) THEN
	      SZ(1) = DCMPLX ( DSQRT ( HPS(IPS) + 
     $	      TSLASK*DSQRT(TSLASK*TSLASK+4*PSPROD(IPS)) ) , 0.0D0 )
	    ELSE
	      SZ(1) = (0.0D0,0.0D0)
	    END IF
	    TSLASK = DMIN1 ( DMAX1(TENDPS(IPS),TPS(IPS)) , TPSCORN )
	    IF (TSLASK.NE.TPS(IPS)) THEN
	      SZ(2) = DCMPLX ( DSQRT ( HPS(IPS) + 
     $	      TSLASK*DSQRT(TSLASK*TSLASK+4*PSPROD(IPS)) ) , 0.0D0 )
	    ELSE
	      SZ(2) = (0.0D0,0.0D0)
	    END IF
	    IF (IADINTXFI.EQ.0) THEN
	      CALL TRAPEZFI (NHYD,LI, DHBRANCHS,RBRANCHS, 
     $	      IFILONLOK,RBRANCHPOSS,GAMMA, 2,SZ, HMAX, NHYDLI,CPRESSW, 
     $	      NFCAL, DHFAKT,RFAKT, ICWORKLIM0,ICWORKDIM,CPRESSW)
	      NFCALLS(3) = NFCALLS(3) + NFCAL
	    ELSE
	      CALL ADINTXFI (NHYD,LI, DHBRANCHS,RBRANCHS, 
     $	      IFILONLOK,RBRANCHPOSS,GAMMA, 2,SZ, HMAX,SMAXPRE,EPSH01Q,
     $	      EPS,NELMAXPRE,NBACKMAX, IXTYP,NRUS,IRUS,JRUS, ARGDIFLIM1,
     $	      ARGDIFLIM0,ARGDIFLIM2, 0.0D0,0,NKVAR, IADINTXFI,TEXTLOG, 
     $	      NHYDLI,CPRESSW, EPSMAGNLOK, NFCAL, DHFAKT,RFAKT, 
     $	      ICWORKLIM0,ICWORKDIM,CPRESSW)
	      NFCALLS(3) = NFCALLS(3) + NFCAL
	      IF (EPSMAGNLOK.GT.SUPEPSMAGN) SUPEPSMAGN = EPSMAGNLOK
	    END IF

	  END IF

	END DO

	IF ( IFILONPRE.NE.0 .AND. IFILONPRE.NE.-3 ) THEN
C	  reset
	  DO IR = 1, LI
	    CSLASK = RK0*OR(IR)
	    GAMMA(IR) = DCMPLX(-DIMAG(CSLASK),DREAL(CSLASK))    ! i*k*r
	  END DO
	END IF

	RETURN
	END
