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 IMPLICIT NONE INCLUDE 'RPPARAMETER.INC' C nol <= noldim , nhyd <= nhyddim , li <= lidim INCLUDE 'INMODEL.INC' INCLUDE 'INRANGE.INC' 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