*
* ..File: xpij1te.f  
*
*
* ..The 2-loop MS(bar) flavour-singlet splitting functions P_ij^(0,1)T
*    for the evolution of unpolarized parton fragmentation densities 
*    at mu_r = mu_f.  The expansion parameter is alpha_s/(4 pi).
*
* ..The distributions (in the mathematical sense) are given as in eq.
*    (B.26) of Floratos, Kounnas, Lacaze: Nucl. Phys. B192 (1981) 417.
*    The name-endings A, B, and C of the functions below correspond to 
*    the kernel superscripts [2], [3], and [1] in that equation.
* 
* ..The code uses the package of Gehrmann and Remiddi for the harmonic
*    polylogarithms published in hep-ph/0107173 = CPC 141 (2001) 296.
*
* ..These are well-known functions, which were re-calulated during the 
*    determinations of their NNLO counterparts presented in 
*
*    S. Moch and A. Vogt,         PL B659 (2008) 290, arXiv:0708.3899 
*    A. Almasy, S. Moch, A. Vogt, NP B854 (2012) 133, arXiv:1107.2263 
*
* =====================================================================
*
*
* ..The 2-loop pure-singlet quark-quark splitting functions P_ps^(1)T.
*   P_qq^(1)T is obtained by adding the non-singlet quantity P_NS^(1)T+
*
       FUNCTION X1PSTA (X, NF)
*
       IMPLICIT REAL*8 (A - Z)
       COMPLEX*16 HC1, HC2, HC3, HC4
       INTEGER NF, NF2, N1, N2, NW, I1, I2, I3, N
       PARAMETER ( N1 = -1, N2 = 1, NW = 2 )
       DIMENSION HC1(N1:N2),HC2(N1:N2,N1:N2),HC3(N1:N2,N1:N2,N1:N2),
     ,           HC4(N1:N2,N1:N2,N1:N2,N1:N2)
       DIMENSION HR1(N1:N2),HR2(N1:N2,N1:N2),HR3(N1:N2,N1:N2,N1:N2),
     ,           HR4(N1:N2,N1:N2,N1:N2,N1:N2)
       DIMENSION HI1(N1:N2),HI2(N1:N2,N1:N2),HI3(N1:N2,N1:N2,N1:N2),
     ,           HI4(N1:N2,N1:N2,N1:N2,N1:N2)
*
* ...Colour factor and an abbreviation
*
       CF = 4./3.D0
       DX = 1.D0/X
*
* ...The harmonic polylogs up to weight 2 by Gehrmann and Remiddi
*
       CALL HPLOG (X, NW, HC1,HC2,HC3,HC4, HR1,HR2,HR3,HR4,
     ,            HI1,HI2,HI3,HI4, N1, N2)
*
* ...The splitting function in terms of the harmonic polylogs
*
      PTqqps1 =
     &  + nf*cf * (  - 32.D0 + 16.D0*x + 224.D0/9.D0*x**2 - 80.D0/9.D0*
     &    dx - 20.D0*Hr1(0) - 36.D0*Hr1(0)*x - 32.D0/3.D0*Hr1(0)*x**2
     &     + 8.D0*Hr2(0,0) + 8.D0*Hr2(0,0)*x )
*
       X1PSTA = PTqqps1
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The 2-loop quark-gluon splitting functions P_qg^(1)T
*
       FUNCTION X1QGTA (X, NF)
*
       IMPLICIT REAL*8 (A - Z)
       COMPLEX*16 HC1, HC2, HC3, HC4
       INTEGER NF, NF2, N1, N2, NW, I1, I2, I3, N
       PARAMETER ( N1 = -1, N2 = 1, NW = 2 )
       DIMENSION HC1(N1:N2),HC2(N1:N2,N1:N2),HC3(N1:N2,N1:N2,N1:N2),
     ,           HC4(N1:N2,N1:N2,N1:N2,N1:N2)
       DIMENSION HR1(N1:N2),HR2(N1:N2,N1:N2),HR3(N1:N2,N1:N2,N1:N2),
     ,           HR4(N1:N2,N1:N2,N1:N2,N1:N2)
       DIMENSION HI1(N1:N2),HI2(N1:N2,N1:N2),HI3(N1:N2,N1:N2,N1:N2),
     ,           HI4(N1:N2,N1:N2,N1:N2,N1:N2)
       PARAMETER ( Z2 = 1.6449 34066 84822 64365 D0 )
*
* ...Colour factors and an abbreviation
*
       CF = 4./3.D0
       CA = 3.D0
       DX = 1.D0/X
*
* ...The harmonic polylogs up to weight 2
*
       CALL HPLOG (X, NW, HC1,HC2,HC3,HC4, HR1,HR2,HR3,HR4,
     ,            HI1,HI2,HI3,HI4, N1, N2)
*
* ...The splitting function in terms of the harmonic polylogs
*
      PTqg1 =
     &  + nf*ca * ( 52.D0/9.D0 - 380.D0/9.D0*x + 712.D0/9.D0*x**2 - 80.D
     &    0/9.D0*dx - 16.D0*z2*x - 16.D0/3.D0*Hr1(0) - 136.D0/3.D0*Hr1(
     &    0)*x - 16.D0/3.D0*Hr1(0)*x**2 - 44.D0/3.D0*Hr1(1) + 40.D0/3.D0
     &    *Hr1(1)*x - 40.D0/3.D0*Hr1(1)*x**2 - 8.D0*Hr2(-1,0) - 16.D0*
     &    Hr2(-1,0)*x - 16.D0*Hr2(-1,0)*x**2 + 8.D0*Hr2(0,0) + 48.D0*
     &    Hr2(0,0)*x - 16.D0*Hr2(0,1) + 32.D0*Hr2(0,1)*x - 32.D0*Hr2(0,
     &    1)*x**2 + 16.D0*Hr2(1,0) - 32.D0*Hr2(1,0)*x + 32.D0*Hr2(1,0)*
     &    x**2 + 8.D0*Hr2(1,1) - 16.D0*Hr2(1,1)*x + 16.D0*Hr2(1,1)*x**2
     &     )
      PTqg1 = PTqg1 + nf*cf * (  - 24.D0 + 46.D0*x - 40.D0*x**2 - 8.D0*
     &    z2 + 16.D0*z2*x - 16.D0*z2*x**2 - 10.D0*Hr1(0) + 8.D0*Hr1(0)*
     &    x + 8.D0*Hr1(0)*x**2 + 12.D0*Hr1(1) - 8.D0*Hr1(1)*x + 8.D0*
     &    Hr1(1)*x**2 - 4.D0*Hr2(0,0) + 8.D0*Hr2(0,0)*x - 16.D0*Hr2(0,0
     &    )*x**2 + 8.D0*Hr2(0,1) - 16.D0*Hr2(0,1)*x + 16.D0*Hr2(0,1)*
     &    x**2 - 24.D0*Hr2(1,0) + 48.D0*Hr2(1,0)*x - 48.D0*Hr2(1,0)*
     &    x**2 - 8.D0*Hr2(1,1) + 16.D0*Hr2(1,1)*x - 16.D0*Hr2(1,1)*x**2
     &     )
      PTqg1 = PTqg1 + nf**2 * (  - 40.D0/9.D0 + 32.D0/9.D0*x - 32.D0/9.D
     &    0*x**2 - 8.D0/3.D0*Hr1(0) + 16.D0/3.D0*Hr1(0)*x - 16.D0/3.D0*
     &    Hr1(0)*x**2 + 8.D0/3.D0*Hr1(1) - 16.D0/3.D0*Hr1(1)*x + 16.D0/
     &    3.D0*Hr1(1)*x**2 )
*
       X1QGTA = PTqg1
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The 2-loop gluon-quark splitting functions P_gq^(1)T
*
       FUNCTION X1GQTA (X, NF)
*
       IMPLICIT REAL*8 (A - Z)
       COMPLEX*16 HC1, HC2, HC3, HC4 
       INTEGER NF, NF2, N1, N2, NW, I1, I2, I3, N
       PARAMETER ( N1 = -1, N2 = 1, NW = 2 ) 
       DIMENSION HC1(N1:N2),HC2(N1:N2,N1:N2),HC3(N1:N2,N1:N2,N1:N2), 
     ,           HC4(N1:N2,N1:N2,N1:N2,N1:N2) 
       DIMENSION HR1(N1:N2),HR2(N1:N2,N1:N2),HR3(N1:N2,N1:N2,N1:N2), 
     ,           HR4(N1:N2,N1:N2,N1:N2,N1:N2) 
       DIMENSION HI1(N1:N2),HI2(N1:N2,N1:N2),HI3(N1:N2,N1:N2,N1:N2), 
     ,           HI4(N1:N2,N1:N2,N1:N2,N1:N2) 
       PARAMETER ( Z2 = 1.6449 34066 84822 64365 D0 )
*
* ...Colour factors and abbreviation
*
       CF = 4./3.D0
       CA = 3.D0
       DX = 1.D0/X
*
* ...The harmonic polylogs up to weight 2
*
       CALL HPLOG (X, NW, HC1,HC2,HC3,HC4, HR1,HR2,HR3,HR4,
     ,            HI1,HI2,HI3,HI4, N1, N2) 
*
* ...The splitting function in terms of the harmonic polylogs
*
      PTgq1 =
     &  + cf*ca * ( 20.D0 - 4.D0*x - 176.D0/9.D0*x**2 + 68.D0/9.D0*dx
     &     + 16.D0*z2 + 32.D0*Hr1(0) + 36.D0*Hr1(0)*x + 32.D0/3.D0*Hr1(
     &    0)*x**2 - 24.D0*Hr1(0)*dx + 8.D0*Hr1(1)*x + 16.D0*Hr2(-1,0)
     &     + 8.D0*Hr2(-1,0)*x + 16.D0*Hr2(-1,0)*dx - 16.D0*Hr2(0,0) -
     &    24.D0*Hr2(0,0)*x - 32.D0*Hr2(0,0)*dx - 16.D0*Hr2(0,1) + 8.D0*
     &    Hr2(0,1)*x + 16.D0*Hr2(0,1)*dx + 48.D0*Hr2(1,0) - 24.D0*Hr2(1
     &    ,0)*x - 48.D0*Hr2(1,0)*dx + 16.D0*Hr2(1,1) - 8.D0*Hr2(1,1)*x
     &     - 16.D0*Hr2(1,1)*dx )
      PTgq1 = PTgq1 + cf**2 * (  - 2.D0 + 18.D0*x - 32.D0*Hr1(0) + 2.D0
     &    *Hr1(0)*x - 8.D0*Hr1(1)*x + 8.D0*Hr2(0,0) - 4.D0*Hr2(0,0)*x
     &     + 32.D0*Hr2(0,1) - 16.D0*Hr2(0,1)*x - 32.D0*Hr2(0,1)*dx - 32.
     &    D0*Hr2(1,0) + 16.D0*Hr2(1,0)*x + 32.D0*Hr2(1,0)*dx - 16.D0*
     &    Hr2(1,1) + 8.D0*Hr2(1,1)*x + 16.D0*Hr2(1,1)*dx )
*   
       X1GQTA = PTgq1
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The 2-loop gluon-gluon splitting function P_gg^(1)T, regular piece 
*
       FUNCTION X1GGTA (X, NF)
*
       IMPLICIT REAL*8 (A - Z)
       COMPLEX*16 HC1, HC2, HC3, HC4 
       INTEGER NF, NF2, N1, N2, NW, I1, I2, I3, N
       PARAMETER ( N1 = -1, N2 = 1, NW = 2 ) 
       DIMENSION HC1(N1:N2),HC2(N1:N2,N1:N2),HC3(N1:N2,N1:N2,N1:N2), 
     ,           HC4(N1:N2,N1:N2,N1:N2,N1:N2) 
       DIMENSION HR1(N1:N2),HR2(N1:N2,N1:N2),HR3(N1:N2,N1:N2,N1:N2), 
     ,           HR4(N1:N2,N1:N2,N1:N2,N1:N2) 
       DIMENSION HI1(N1:N2),HI2(N1:N2,N1:N2),HI3(N1:N2,N1:N2,N1:N2), 
     ,           HI4(N1:N2,N1:N2,N1:N2,N1:N2) 
       PARAMETER ( Z2 = 1.6449 34066 84822 64365 D0,
     ,             Z3 = 1.2020 56903 15959 42854 D0 )
*
* ..The soft coefficient for use in X1GGB and X1GGC
*
       COMMON / P1GSOFT / A2G
*
* ...Colour factors
*
       CF  = 4./3.D0
       CA  = 3.D0
*
* ...Some abbreviations
*
       DX = 1.D0/X
       DM = 1.D0/(1.D0-X)
       DP = 1.D0/(1.D0+X)
*
* ...The harmonic polylogs up to weight 2
*
       CALL HPLOG (X, NW, HC1,HC2,HC3,HC4, HR1,HR2,HR3,HR4,
     ,            HI1,HI2,HI3,HI4, N1, N2) 
*
* ...The splitting function in terms of the harmonic polylogs
*    (without the delta(1-x) part, but with the soft contribution)
*
      PTgg1 =
     &  + ca**2 * (  - 50.D0/9.D0 - 218.D0/9.D0*x + 268.D0/9.D0*dm + 32.
     &    D0*z2 + 16.D0*z2*x**2 - 8.D0*z2*dp - 8.D0*z2*dm - 44.D0*Hr1(0
     &    ) - 4.D0*Hr1(0)*x - 88.D0/3.D0*Hr1(0)*x**2 - 88.D0/3.D0*Hr1(0
     &    )*dx + 88.D0/3.D0*Hr1(0)*dm + 32.D0*Hr2(-1,0) + 16.D0*Hr2(-1,
     &    0)*x + 16.D0*Hr2(-1,0)*x**2 + 16.D0*Hr2(-1,0)*dx - 16.D0*Hr2(
     &    -1,0)*dp - 64.D0*Hr2(0,0)*x + 16.D0*Hr2(0,0)*x**2 - 32.D0*
     &    Hr2(0,0)*dx + 8.D0*Hr2(0,0)*dp - 24.D0*Hr2(0,0)*dm + 32.D0*
     &    Hr2(0,1) - 16.D0*Hr2(0,1)*x + 16.D0*Hr2(0,1)*x**2 - 16.D0*
     &    Hr2(0,1)*dx - 16.D0*Hr2(0,1)*dm + 32.D0*Hr2(1,0) - 16.D0*Hr2(
     &    1,0)*x + 16.D0*Hr2(1,0)*x**2 - 16.D0*Hr2(1,0)*dx - 16.D0*Hr2(
     &    1,0)*dm )
      PTgg1 = PTgg1 + nf*ca * ( 116.D0/9.D0 - 76.D0/9.D0*x + 92.D0/9.D0
     &    *x**2 - 92.D0/9.D0*dx - 40.D0/9.D0*dm + 8.D0*Hr1(0) - 8.D0*
     &    Hr1(0)*x + 16.D0/3.D0*Hr1(0)*x**2 - 16.D0/3.D0*Hr1(0)*dx - 16.
     &    D0/3.D0*Hr1(0)*dm )
      PTgg1 = PTgg1 + nf*cf * (  - 8.D0 + 24.D0*x - 328.D0/9.D0*x**2 +
     &    184.D0/9.D0*dx + 20.D0*Hr1(0) + 28.D0*Hr1(0)*x + 32.D0/3.D0*
     &    Hr1(0)*x**2 + 32.D0/3.D0*Hr1(0)*dx + 8.D0*Hr2(0,0) + 8.D0*
     &    Hr2(0,0)*x )
*
* ...The soft (`+'-distribution) part of the splitting function
*
       A2G   = - 40.D0/9.D0*ca*nf + 268.D0/9.D0*ca**2 - 8.D0*z2*ca**2
       GGG1L = DM * A2G
*
* ...The regular piece of the coefficient function
*
       X1GGTA = PTgg1 - GGG1L
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The singular (soft) piece of P_gg^(1)T  (as in the spacelike case)
*
       FUNCTION X1GGTB (Y, NF)
       IMPLICIT REAL*8 (A - Z)
       INTEGER NF
*
       COMMON / P1GSOFT / A2G
*
       X1GGTB  = A2G/(1.D0-Y)
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The 'local' piece of P_gg^(1)T  (as in the spacelike case)
*
       FUNCTION X1GGTC (Y, NF)
*
       IMPLICIT REAL*8 (A - Z)
       INTEGER NF
       PARAMETER ( Z2 = 1.6449 34066 84822 64365 D0,
     ,             Z3 = 1.2020 56903 15959 42854 D0 )
*
       COMMON / P1GSOFT / A2G
*
* ...Colour factors
*
       CF  = 4./3.D0
       CA  = 3.D0
*
* ...The coefficient of delta(1-x)
*
       P1DELT = 
     ,    - 2.D0*cf*nf
     ,    - 8.D0/3.D0*ca*nf
     ,    + 32.D0/3.D0*ca**2
     ,    + 12.D0*z3*ca**2
*
       X1GGTC = DLOG (1.D0-Y) * A2G + P1DELT
*
       RETURN
       END
*
* =================================================================av==
