*
* ..File: xppf2p.f   
*
* ..Parametrisation of the O(a_em a_s^2) MS(bar) photon-gluon and
*    (flavour singlet and non-singlet) photon-quark splitting functions
*    P_pf^(2), p = g, ns, q,  at mu_r = mu_f.  The expansion parameters
*    are normalized as  a_em = alpha_em/(4 pi),  a_s = alpha_s/(4 pi).
*
* ..The charge prefactors, e.g.,  delta_s = 3 N_f <e^2>,  have to be
*    supplemented by the calling routines.
*
* ..The two sets providing the error estimate are called via  IMOD = 1
*    and  IMOD = 2.  Any other value of IMOD invokes their average.
*    The results are based on the lowest six even-integer moments. 
*
* ..Reference: S. Moch, J. Vermaseren and A. Vogt, 
*              hep-ph/0110331 = Nucl. Phys. B621 (2002) 413 
*
* 
* =====================================================================
*
*
* ..This is the non-singlet photon->quark splitting function P_nsf^(2).
*
       FUNCTION P2NSFA (Y, NF, IMOD)
*
       IMPLICIT REAL*8 (A-Z)
       INTEGER IMOD, NF
       DL  = LOG (Y)
       DL1 = LOG (1.-Y)
*
       P2NSFA1 =   3.8636 * DL1**3 + 97.512 * DL1**2 - 1319.749 * Y**2
     1          + 511.199 * Y + 84.808 * DL - 22.878 * DL**3
       P2NSFA2 = - 5.4658 * DL1**3 - 295.331 * DL1 - 1803.989 * Y
     1          + 740.532 - 259.036 * DL**2 + 27.110 * DL**4
*
       IF (IMOD .EQ. 1) THEN
         P2NSFA = P2NSFA1
       ELSE IF (IMOD .EQ. 2) THEN
         P2NSFA = P2NSFA2
       ELSE
         P2NSFA = 0.5 * (P2NSFA1 + P2NSFA2)
       END IF
*
       P2NSFA = P2NSFA + 128./27.D0 * DL1**4 
     1         + NF * (- 0.2949 * DL1**3 + 34.854 * DL1 + 157.995 * Y
     2                 - 73.672 - 33.059 * DL + 2.887 * DL**3 )
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..This is the pure-singlet photon->quark splitting function P_psf^(2).
*    P_qf^(2) is obtained by adding the non-singlet quantity P_nsf^(2).
*
       FUNCTION P2PSFA (Y, NF, IMOD)
*
       IMPLICIT REAL*8 (A-Z)
       DIMENSION CL3(2)
       INTEGER IMOD, NF
       DL  = LOG (Y)
       DL1 = LOG (1.-Y)
*
       P2PSFA1 = (- 110.167 * Y**2 + 876.629 * Y + 23.605 / Y) * (1.-Y)
     1          + 668.725 * Y*DL + 387.125 * Y*DL**2 + 121.403 * DL  
       P2PSFA2 = (- 34.853 * Y + 305.244 + 50.950 / Y) * (1.-Y)
     1          + 101.246 * Y*DL + 220.083 * DL - 5.0738 * DL**4 
*
       IF (IMOD .EQ. 1) THEN
         P2PSFA = P2PSFA1
       ELSE IF (IMOD .EQ. 2) THEN
         P2PSFA = P2PSFA2
       ELSE
         P2PSFA = 0.5 * (P2PSFA1 + P2PSFA2)
       END IF
*
       P2PSFA = NF * P2PSFA
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..This is the photon->gluon splitting function P_gf^(2).
*
       FUNCTION P2GFA (Y, NF, IMOD)
*
       IMPLICIT REAL*8 (A-Z)
       INTEGER IMOD, NF
       DL  = LOG (Y)
       DL1 = LOG (1.-Y)
*
       P2GFA1 = (769.616 * DL1 + 1329.961 * Y**2 - 391.569 * Y)* (1.-Y) 
     1          + 317.267 * DL**2 + 60.519 * DL**3 + 15.018 * DL/Y 
       P2GFA2 = (-105.632 * DL1**2 - 415.549 * Y**2 - 429.907) * (1.-Y)
     1          - 357.604 * DL - 146.286 * DL**3 + 64.666 * DL/Y
*
       IF (IMOD .EQ. 1) THEN
         P2GFA = P2GFA1
       ELSE IF (IMOD .EQ. 2) THEN
         P2GFA = P2GFA2
       ELSE
         P2GFA = 0.5 * (P2GFA1 + P2GFA2)
       END IF
*
       P2GFA = P2GFA 
     1         + NF * ( (1.-Y) * (43.748 * DL1 + 61.028 * Y**2 - 70.910 
     2                  - 33.913 /Y) - 105.172 * DL + 1.3972 * DL**2 )
*
       RETURN
       END
*
* =================================================================av==
