pyspen.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyspen.f

 
C***********************************************************************
 
C...PYSPEN
C...Calculates real and imaginary part of Spence function; see
C...G. 't Hooft and M. Veltman, Nucl. Phys. B153 (1979) 365.
 
      FUNCTION pyspen(XREIN,XIMIN,IREIM)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Commonblocks.
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      SAVE /pydat1/
C...Local array and data.
      dimension b(0:14)
      DATA b/
     &1.000000d+00,        -5.000000d-01,         1.666667d-01,
     &0.000000d+00,        -3.333333d-02,         0.000000d+00,
     &2.380952d-02,         0.000000d+00,        -3.333333d-02,
     &0.000000d+00,         7.575757d-02,         0.000000d+00,
     &-2.531135d-01,         0.000000d+00,         1.166667d+00/
 
      xre=xrein
      xim=ximin
      IF(abs(1d0-xre).LT.1d-6.AND.abs(xim).LT.1d-6) THEN
        IF(ireim.EQ.1) pyspen=paru(1)**2/6d0
        IF(ireim.EQ.2) pyspen=0d0
        RETURN
      ENDIF
 
      xmod=sqrt(xre**2+xim**2)
      IF(xmod.LT.1d-6) THEN
        IF(ireim.EQ.1) pyspen=0d0
        IF(ireim.EQ.2) pyspen=0d0
        RETURN
      ENDIF
 
      xarg=sign(acos(xre/xmod),xim)
      sp0re=0d0
      sp0im=0d0
      sgn=1d0
      IF(xmod.GT.1d0) THEN
        algxre=log(xmod)
        algxim=xarg-sign(paru(1),xarg)
        sp0re=-paru(1)**2/6d0-(algxre**2-algxim**2)/2d0
        sp0im=-algxre*algxim
        sgn=-1d0
        xmod=1d0/xmod
        xarg=-xarg
        xre=xmod*cos(xarg)
        xim=xmod*sin(xarg)
      ENDIF
      IF(xre.GT.0.5d0) THEN
        algxre=log(xmod)
        algxim=xarg
        xre=1d0-xre
        xim=-xim
        xmod=sqrt(xre**2+xim**2)
        xarg=sign(acos(xre/xmod),xim)
        algyre=log(xmod)
        algyim=xarg
        sp0re=sp0re+sgn*(paru(1)**2/6d0-(algxre*algyre-algxim*algyim))
        sp0im=sp0im-sgn*(algxre*algyim+algxim*algyre)
        sgn=-sgn
      ENDIF
 
      xre=1d0-xre
      xim=-xim
      xmod=sqrt(xre**2+xim**2)
      xarg=sign(acos(xre/xmod),xim)
      zre=-log(xmod)
      zim=-xarg
 
      spre=0d0
      spim=0d0
      savere=1d0
      saveim=0d0
      DO 100 i=0,14
        IF(max(abs(savere),abs(saveim)).LT.1d-30) goto 110
        termre=(savere*zre-saveim*zim)/dble(i+1)
        termim=(savere*zim+saveim*zre)/dble(i+1)
        savere=termre
        saveim=termim
        spre=spre+b(i)*termre
        spim=spim+b(i)*termim
  100 CONTINUE
 
  110 IF(ireim.EQ.1) pyspen=sp0re+sgn*spre
      IF(ireim.EQ.2) pyspen=sp0im+sgn*spim
 
      RETURN
      END