*
* ..Moments of the timelike (even-N based) splitting functions,
*   using an old code prepared for the GRV fragmentation paper
*
        program NsplitT
        implicit complex (a - z)
        integer nf, k1
*
        nf = 5
*
        do 1 k1 = 4, 4
*
        XN = cmplx(2.e0*k1, 0.e0)
        CALL FRCALC (QQI, QGF, GQI, GGI, GGF, NS1PI, NS1F, 
     1            QQ1F, QG1F, GQ1I, GQ1F, GG1I, GG1F, QPF, QP1F, 
     2            GP1F, C2QI, C1QI, C2GF, C1GF, C2PI, C1PI, XN)
*
* ..Conversion from the 1992 to the 2011 notation
*
        Pqq0N = -QQI/2.e0
        Pqg0N = -GQI/2.e0 * (2.*nf)
        Pgq0N = -QGF*nf/2.e0 / (2.*nf)
        Pgg0N = -(GGI + GGF*nf)/2.e0
*
        Pqq1N = -(NS1PI + (NS1F+QQ1F)*nf)/2.e0
        Pqg1N = -(GQ1I + GQ1F*nf)/2.e0 * (2.*nf)
        Pgq1N = -(QG1F*nf)/2.e0 / (2.*nf)
        Pgg1N = -(GG1I+GG1F*nf)/2.e0
*
* ..Output, bare bones
*
        write(6,*) XN, Pqq0N
        write(6,*) XN, Pqg0N
        write(6,*) XN, Pgq0N
        write(6,*) XN, Pgg0N
*
        write(6,*) XN, Pqq1N
        write(6,*) XN, Pqg1N
        write(6,*) XN, Pgq1N
        write(6,*) XN, Pgg1N
*
   1    continue
        stop
        end
*
*
* ..Original 1992 code, but with the improved parametrization of 
*   Sp(x)/(1+x) in N-space as used in the QCD-Pegasus package
C                                                  
C ..CALCULATION OF ANOMALOUS DIMENSIONS (AND WILSON COEFFICIENTS)      
C ..UP TO THEIR DEPENDENCE OF THE NUMBER OF ACTIVE FLAVOURS F :        
       SUBROUTINE FRCALC (QQI, QGF, GQI, GGI, GGF,        NS1PI, NS1F, 
     1               QQ1F, QG1F, GQ1I, GQ1F, GG1I, GG1F, QPF, QP1F,    
     2               GP1F, C2QI, C1QI, C2GF, C1GF, C2PI, C1PI, N)      
       IMPLICIT COMPLEX (A - Z)                                       
       REAL ZETA2, ZETA3                                             
       NS = N * N                                                    
       N1 = N + 1.                                                   
       N2 = N + 2.                                                   
       NM = N - 1.                                                   
C ..LEADING ORDER :                                                  
       S1  = PSIFN (N1) + 0.577216                                   
       QQI = (8./3.) * (-3.- 2./(N * N1) + 4.* S1)                   
       GQI = -2.* (NS + N +2.) / (N * N1 * N2)                      
       QGF = -(32./3.) * (NS + N + 2.) / (N * N1 * NM)              
       GGI = -22.- 24./(N * NM) - 24./(N1 * N2) + 24.* S1           
       GGF = 4./3.                                                  
C ..NEXT TO LEADING ORDER :                                         
       NT  = NS * N                                                 
       NFO = NT * N                                                 
       N1S = N1 * N1                                                
       N1T = N1S * N1                                               
C ..ANALYTIC CONTINUATIONS OF N-SUMS AS GIVEN IN GLUECK ET AL. (1990) :
       ZETA2 = 1.644934                                                
       ZETA3 = 1.202057                                               
       S11 = S1 + 1./N1                                              
       S2  = ZETA2 - PSIFN1 (N1)                                    
       S21 = S2 + 1./N1S                                           
       G1  = 0.5 * (PSIFN1 (N1/2.) - PSIFN1 (N/2.))               
       G11 = 0.5 * (PSIFN1 (N2/2.) - PSIFN1 (N1/2.))             
*      SPMOM = 1.01/N1 - 0.846/N2 + 1.155/(N+3.) - 1.074/(N+4.) +   
*    1         0.55/(N+5.)                                          
*
* ..New: use an improved parametrization of the moments of Sp(x)/(1+x)
*        for higher accuracy
*
       N3 = N + 3.
       N4 = N + 4.
       N5 = N + 5.
       N6 = N + 6.
       S12 = S11 + 1./N2
       S13 = S12 + 1./N3
       S14 = S13 + 1./N4
       S15 = S14 + 1./N5
       S16 = S15 + 1./N6
       SPMOM = 1.0000D0 * (ZETA2 - S1 / N ) / N  -
     1         0.9992D0 * (ZETA2 - S11/ N1) / N1 +
     2         0.9851D0 * (ZETA2 - S12/ N2) / N2 -
     3         0.9005D0 * (ZETA2 - S13/ N3) / N3 +
     4         0.6621D0 * (ZETA2 - S14/ N4) / N4 -
     5         0.3174D0 * (ZETA2 - S15/ N5) / N5 +
     6         0.0699D0 * (ZETA2 - S16/ N6) / N6
* ..End new block
*
       SLC   = -5./8.* ZETA3                                       
       SLV   = - ZETA2/2.* (PSIFN (N1/2.) - PSIFN (N/2.))          
     1         + S1/NS + SPMOM                                     
c      SSCHLM = SLC - SLV                                         
c      SSTR2M = ZETA2 - PSIFN1 (N1/2.)                            
c      SSTR3M = 0.5 * PSIFN2 (N1/2.) + ZETA3                       
       SSCHLP = SLC + SLV                                          
       SSTR2P = ZETA2 - PSIFN1 (N2/2.)                             
       SSTR3P = 0.5 * PSIFN2 (N2/2.) + ZETA3                       
C...NON-SINGLET PIECES AS GIVEN IN CURCI ET AL. (1980),            
c      NS1MA = 16.* S1 * (2.* N + 1.) / (NS * N1S) +               
c    1         16.* (2.* S1 - 1./(N * N1)) * ( S2 - SSTR2M ) +     
c    2         64.* SSCHLM + 24.* S2 - 3. - 8.* SSTR3M -          
c    3         8.* (3.* NT + NS -1.) / (NT * N1T) +              
c    4         16.* (2.* NS + 2.* N +1.) / (NT * N1T)           
       NS1PA = 16.* S1 * (2.* N + 1.) / (NS * N1S) +           
     1         16.* (2.* S1 - 1./(N * N1)) * ( S2 - SSTR2P ) +        
     2         64.* SSCHLP + 24.* S2 - 3. - 8.* SSTR3P -              
     3         8.* (3.* NT + NS -1.) / (NT * N1T) -                  
     4         16.* (2.* NS + 2.* N +1.) / (NT * N1T)               
       NS1B  = S1 * (536./9. + 8.* (2.* N + 1.) / (NS * N1S)) -    
     1         (16.* S1 + 52./3.- 8./(N * N1)) * S2 - 43./6. -    
     2         (151.* NFO + 263.* NT + 97.* NS + 3.* N + 9.) *   
     3         4./ (9.* NT * N1T)                               
       NS1C  = -160./9.* S1 + 32./3.* S2 + 4./3. +             
     1         16.* (11.* NS + 5.* N - 3.) / (9.* NS * N1S)   
       DNS1I = - 8.* S1 * S2 + 6.* S2 + 4.* (S2 - ZETA2) / (N * N1)   
     1         + 8.* ZETA2 * S1 + 4.* S1 * (2.* N + 1.) / (NS * N1S)  
     2         - (6.* NT + 9.* NS + 7.* N + 2.) / (NT * N1T)          
     3         - 6.* ZETA2                                           
c      NS1MI = -2./9.* NS1MA + 4.* NS1B                             
       NS1PI = -2./9.* NS1PA + 4.* NS1B - 128./9.* DNS1I           
       NS1F  = 2./3. * NS1C                                       
C ..SPACELIKE SINGLET PIECES AS GIVEN IN FLORATOS ET AL. (1981), 
C ..DIFFERENCE OF THE TIMELIKE PIECES AS INVERTED FROM FURMANSKI AND 
C ..PETRONZIO (1980) BY GRV (1992) : 
       NFI = NFO * N                                              
       NSI = NFI * N                                                
       NSE = NSI * N                                             
       NE  = NSE * N                                             
       NN  = NE * N                                              
       NMS = NM * NM                                             
       NMT = NMS * NM                                                 
       N2S = N2 * N2                                                 
       N2T = N2S * N2                                                 
       QQ1FS = (5.* NFI + 32.* NFO + 49.* NT + 38.* NS + 28.* N     
     1           + 8.) / (NM * NT * N1T * N2S) * (-32./3.)         
       DQQ1F = - 80./9./ NM - 12./N + 12./NS + 8./NT - 4./N1      
C .......ERROR IN THE REPRINT (NOT IN THE PREPRINT) OF FP :      
C ..........-8./NS (WRONG) INSTEAD OF 12./NS (CORRECT).         
     1         + 28./N1S + 8./N1T + 224./9./ N2 + 32./3./ N2S  
       QQ1F  = QQ1FS - 16./3.* DQQ1F                          
       GQ1A = 8./3.* S11 * (NS + N + 2.) / (N * N1 * N2)     
     1        + 8./3. * (- 5./(3.* N) + 1./NS - 1./(N * N1)  
     2          + 4./(3.* N1) - 2./N1S - 4./(3.* N2) + 4./N2S)        
       GQ1B = - (2.* S11*S11 - 2.* S11 - 10.* S21) * (NS + N + 2.)   
     1          / (N * N1 * N2)                                     
     2        + 4.* S11 * (1./N - 1./NS + 1./(N * N1) + 2./N1S     
     3          - 4./N2S)                                         
     4        - 12./N + 5./NS - 6./(N * N1) - 12./(N * N1S)      
     5        + 4./(NS * N1) - 2./NT + 23./N1 - 4./N1S + 4./N1T 
     6        - 20./N2                                          
       GQ1C  = (2.* S11 * S11 - 10./3.* S11 - 6.* S21 + 2.* G11      
     1           -6.* ZETA2) * (NS + N + 2.) / (N * N1 * N2)        
     2         - 4.* S11 * (1./N - 2./NS + 1./(N * N1) + 4./N1S     
     3           - 6./N2S)                                          
     4         - 40./9./ NM + 26./9./ N + 8./3./ NS - 190./9./ N1  
     5         + 22./3./(N * N1) + 68./3./ N1S + 16./N1T + 4./NT  
     6         + 8./(N1S * N2) + 356/9./ N2 - 4./N2S             
     7         - 8./(NS * N1S)                                  
C ..GQ1 OF FP HAS TO BE DIVIDED BY  2 * F                      
       GQ1F = - 2.* GQ1A / 2.                                 
       GQ1I = (- 16./3. * GQ1B - 12.* GQ1C) / 2.             
       QG1A  = (S1 * S1 - 3.* S2 - 4.* ZETA2) * (2./NM - 2./N + 1./N1) 
     1         + 2.* S1 * (-2./(NM * N) -1./N1 + 3./N1S + 4./NMS      
     2           - 4./NS)                                            
     3         + 8./(N * NM) - 1./2./ N + 9./2./ N1 - 5./2./ N1S    
     4         + 1./N1T - 8./(NMS * N) + 2./NT                     
       QG1B  = - (S1 * S1 - 5.* S2 + G1 - ZETA2) * (NS + N + 2.)  
     1           / (NM * N * N1)                                 
     2         + 2.* S1 * (NFI - NFO + NT - 9.* NS - 2.* N + 2.)
     3           / (NMS * NS * N1S)                            
     4         - 8./NMT + 6./NMS + 4./(NMS * N) + 17./(9.* NM)
     5         - 12./(NM * NS) + 5./N - 8./NS - 2./(NS * N1)           
     6         - 1./N1 - 7./N1S - 2./N1T - 44./9./ N2 - 8./3./ N2S    
       QG1F  = (- 128./9.* QG1A - 32.* QG1B) * 2.                    
C ..QG1 OF FP HAS TO BE MULTIPLIED BY 2 * F                         
       GG1AS = 16./9.* (38.* NFO + 76.* NT + 94.* NS + 56.* N + 12.)   
     1           / (NM * NS * N1S * N2)   -   160./9.* S1 + 32./3.    
       DGG1A = S2 - 1./NMS + 1./NS - 1./N1S + 1./N2S - ZETA2         
       GG1A  = GG1AS + 64./3.* DGG1A                                
       GG1BS = (2.* NSI + 4.* NFI + NFO - 10.* NT - 5.* NS - 4.* N 
     1          - 4.) * 16. / (NM * NT * N1T * N2)   +   8.       
       DGG1B = 80./9./ NM - 16./3./ NMS + 12./N - 16./NS + 8./NT 
     1         + 4./N1 - 24./N1S + 8./N1T - 224./9./ N2         
     2         - 16./3./ N2S                                   
       GG1B  = GG1BS - 8.* DGG1B                                       
       GG1CS = (2.* NFI + 5.* NFO + 8.* NT + 7.* NS - 2.* N - 2.)     
     1           * 64.* S1 / (NMS * NS * N1S * N2S)                  
     2         + 536./9.* S1 - 64./3.                               
     3         + 32.* SSTR2P * (NS + N + 1.) / (NM * N * N1 * N2)  
     4         - 16.* S1 * SSTR2P + 32.* SSCHLP - 4.* SSTR3P      
     5         - 4.* (457.* NN + 2742.* NE + 6040.* NSE + 6098.* NSI  
     6          + 1567.* NFI - 2344.* NFO - 1632.* NT + 560.* NS     
     7          + 1488.* N + 576.) / (9.* NMS * NT * N1T * N2T)     
       DGG1C = - S1 * S2 + S1 * (ZETA2 + 1./NMS - 1./NS + 1./N1S   
     1         - 1./N2S)                                               
     2         + (S2 - ZETA2) * (11/12.+ 1./NM - 1./N + 1./N1 - 1./N2)
     3         - 1./NMT + 11./12./ NMS - 1./(NMS * N)                
     4         - 1./(NM * NS) - 7./12./ NS - 1./NT + 7./12./ N1S    
     5         - 1./ N1T - 1./(N1S * N2) - 1./(N1 * N2S)           
     6         - 11./12./ N2S - 1./N2T                            
       GG1C  = GG1CS - 64.* DGG1C                                
       GG1I  = 9.* GG1C                                         
       GG1F  = 3./2.* GG1A  + 2./3.* GG1B                      
C ..PHOTON PIECES AS GIVEN IN BARDEEN AND BURAS (1979) :      
       QPF  = -3./4.* QGF                                    
       QP1F = 64./3.* QG1A                                  
       GP1F = 4.* (-4./N + 12./N1 - 164./9./ N2 + 92./9./ NM - 10./NS 
     1          -14./N1S - 16./3./ N2S - 16./3./ NMS + 4./NT + 4./N1T)
C ..TIMELIKE WILSON COEFFICIENTS :                                    
       C2QI = 4./3.* (2.* S1 * S1 - 2.* S2 + 3.* S1 - 9.             
     1       - 2.* S1 / (N * N1) + 3./ N + 4./ N1 + 2./ NS)         
       DC2Q = 8./3.* (6.* S2 + 3./N1S - 3./NS - 7./2./ N - 7./2./ N1)  
       C2QI = C2QI + DC2Q                                             
       C1QI = C2QI + 16./3./ N                                       
       C2GF = (-S1 * (NS + N + 2.) / (NM * N * N1) - 4./(NM * N)    
     1         + 4.* (1.- 2.* N) / (NMS * NS) - 3./ N1S) * 16./3.  
       C1GF = C2GF + 16./3.* (4./NM - 4./N)                       
       C2PI = (- S1 * (NS + N + 2.) / (NM * N * N1) - 4./(NM * N)  
     1         + 4.* (1.- 2.* N) / (NMS * NS) - 3./ N1S) * 4.     
       C1PI = C2PI + 4.* (4./NM - 4./N)                          
       RETURN                                                   
       END                                                     
C                                                             
C ..PSI - FUNCTION FOR COMPLEX ARGUMENT                      
       COMPLEX FUNCTION PSIFN (Z)                           
       COMPLEX Z, ZZ, RZ, DZ, SUB, CLOG                    
       SUB = CMPLX (0.,0.)                                
       ZZ = Z                                            
  1    CONTINUE                                         
       IF (REAL (ZZ) .LT. 10.) THEN                    
         SUB = SUB - 1./ ZZ                           
         ZZ = ZZ + 1.                                
         GOTO 1                                     
       END IF                                                         
       RZ = 1./ ZZ                                                   
       DZ = RZ * RZ                                                 
       PSIFN = SUB + CLOG(ZZ) - RZ/2.- DZ/2520. * ( 210.+ DZ * (-21.+ 
     1         10.*DZ ))                                              
       RETURN                                                        
       END                                                          
C                                                                  
C ..FIRST DERIVATIVE OF THE PSI - FUNCTION FOR COMPLEX ARGUMENT : 
       COMPLEX FUNCTION PSIFN1 (Z)                               
       COMPLEX Z, ZZ, RZ, DZ, SUB                               
       SUB = CMPLX (0.,0.)                                     
       ZZ = Z                                                 
  1    CONTINUE                                              
       IF (REAL (ZZ) .LT. 10.) THEN                         
         SUB = SUB + 1./ (ZZ * ZZ)                         
         ZZ = ZZ + 1.                                     
         GOTO 1                                          
       END IF                                           
       RZ = 1./ ZZ                                     
       DZ = RZ * RZ                                   
       PSIFN1 = SUB + RZ + DZ/2. * ( 1 + RZ/630. * ( 210.- DZ * ( 42.- 
     1         DZ * ( 30.- 42.*DZ ))))                                 
       RETURN                                                         
       END                                                           
C                                                                   
C ..SECOND DERIVATIVE OF THE PSI - FUNCTION FOR COMPLEX ARGUMENT : 
       COMPLEX FUNCTION PSIFN2 (Z)                                
       COMPLEX Z, ZZ, RZ, DZ, SUB                                
       SUB = CMPLX (0.,0.)                                      
       ZZ = Z                                                  
  1    CONTINUE                                               
       IF (REAL (ZZ) .LT. 10.) THEN                          
         SUB = SUB - 2./ (ZZ * ZZ * ZZ)                     
         ZZ = ZZ + 1.                                      
         GOTO 1                                           
       END IF                                            
       RZ = 1./ ZZ                                      
       DZ = RZ * RZ                                    
       PSIFN2 = SUB - DZ/60. * ( 60.+ RZ * ( 60.+ RZ * ( 30.- DZ * 
     1         ( 10.- DZ * ( 10.- DZ * ( 18.- 50.* DZ ))))))      
       RETURN                                                    
       END  
