pythiaeRHIC/pythia/pyspen.f

0001
0002 C***********************************************************************
0003
0004 C...PYSPEN
0005 C...Calculates real and imaginary part of Spence function; see
0006 C...G. 't Hooft and M. Veltman, Nucl. Phys. B153 (1979) 365.
0007
0008       FUNCTION PYSPEN(XREIN,XIMIN,IREIM)
0009
0010 C...Double precision and integer declarations.
0011       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0012       IMPLICIT INTEGER(I-N)
0013       INTEGER PYK,PYCHGE,PYCOMP
0014 C...Commonblocks.
0015       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0016       SAVE /PYDAT1/
0017 C...Local array and data.
0018       DIMENSION B(0:14)
0019       DATA B/
0020      &1.000000D+00,        -5.000000D-01,         1.666667D-01,
0021      &0.000000D+00,        -3.333333D-02,         0.000000D+00,
0022      &2.380952D-02,         0.000000D+00,        -3.333333D-02,
0023      &0.000000D+00,         7.575757D-02,         0.000000D+00,
0024      &-2.531135D-01,         0.000000D+00,         1.166667D+00/
0025
0026       XRE=XREIN
0027       XIM=XIMIN
0028       IF(ABS(1D0-XRE).LT.1D-6.AND.ABS(XIM).LT.1D-6) THEN
0029         IF(IREIM.EQ.1) PYSPEN=PARU(1)**2/6D0
0030         IF(IREIM.EQ.2) PYSPEN=0D0
0031         RETURN
0032       ENDIF
0033
0034       XMOD=SQRT(XRE**2+XIM**2)
0035       IF(XMOD.LT.1D-6) THEN
0036         IF(IREIM.EQ.1) PYSPEN=0D0
0037         IF(IREIM.EQ.2) PYSPEN=0D0
0038         RETURN
0039       ENDIF
0040
0041       XARG=SIGN(ACOS(XRE/XMOD),XIM)
0042       SP0RE=0D0
0043       SP0IM=0D0
0044       SGN=1D0
0045       IF(XMOD.GT.1D0) THEN
0046         ALGXRE=LOG(XMOD)
0047         ALGXIM=XARG-SIGN(PARU(1),XARG)
0048         SP0RE=-PARU(1)**2/6D0-(ALGXRE**2-ALGXIM**2)/2D0
0049         SP0IM=-ALGXRE*ALGXIM
0050         SGN=-1D0
0051         XMOD=1D0/XMOD
0052         XARG=-XARG
0053         XRE=XMOD*COS(XARG)
0054         XIM=XMOD*SIN(XARG)
0055       ENDIF
0056       IF(XRE.GT.0.5D0) THEN
0057         ALGXRE=LOG(XMOD)
0058         ALGXIM=XARG
0059         XRE=1D0-XRE
0060         XIM=-XIM
0061         XMOD=SQRT(XRE**2+XIM**2)
0062         XARG=SIGN(ACOS(XRE/XMOD),XIM)
0063         ALGYRE=LOG(XMOD)
0064         ALGYIM=XARG
0065         SP0RE=SP0RE+SGN*(PARU(1)**2/6D0-(ALGXRE*ALGYRE-ALGXIM*ALGYIM))
0066         SP0IM=SP0IM-SGN*(ALGXRE*ALGYIM+ALGXIM*ALGYRE)
0067         SGN=-SGN
0068       ENDIF
0069
0070       XRE=1D0-XRE
0071       XIM=-XIM
0072       XMOD=SQRT(XRE**2+XIM**2)
0073       XARG=SIGN(ACOS(XRE/XMOD),XIM)
0074       ZRE=-LOG(XMOD)
0075       ZIM=-XARG
0076
0077       SPRE=0D0
0078       SPIM=0D0
0079       SAVERE=1D0
0080       SAVEIM=0D0
0081       DO 100 I=0,14
0082         IF(MAX(ABS(SAVERE),ABS(SAVEIM)).LT.1D-30) GOTO 110
0083         TERMRE=(SAVERE*ZRE-SAVEIM*ZIM)/DBLE(I+1)
0084         TERMIM=(SAVERE*ZIM+SAVEIM*ZRE)/DBLE(I+1)
0085         SAVERE=TERMRE
0086         SAVEIM=TERMIM
0087         SPRE=SPRE+B(I)*TERMRE
0088         SPIM=SPIM+B(I)*TERMIM
0089   100 CONTINUE
0090
0091   110 IF(IREIM.EQ.1) PYSPEN=SP0RE+SGN*SPRE
0092       IF(IREIM.EQ.2) PYSPEN=SP0IM+SGN*SPIM
0093
0094       RETURN
0095       END