C----------------------------------------------------------------------- C Implementation of Michel spectrum functions. C Radiative corrections in different approximations. C C Andrej Arbuzov C C $Id: fg_nla.f,v 1.3 2003/07/18 23:22:53 andr Exp $ C C----------------------------------------------------------------------- C.... The muon decay differential width reads: C \dd^2\Gamma^{\mu^{-}\to ...}/(\dd x\dd c)=\Gamma_0(FFF(x) + c\xi GGG(x)) C "c" ---> "-c" for \mu^{+} decay C C \Gamma_0 &=& \frac{G_F^2 m_{\mu}^5}{192\pi^3} \\ C &\times& ( 1 + \frac{3m_\mu^2}{5m_W^2} ) C C PARAM(1) = RHOM !(\rho Michel) C PARAM(2) = ETAM !(\eta Michel) C PARAM(3) = CHIM !(\chi Michel) C PARAM(4) = DELTAM !(\delta Michel) C PARAM(5) = (M_E/M_MU)**2 != AME2/AMU2 C ALME = DLOG(AMU2/AME2) C PARAM(6) = ALME C PARAM(7) = ALPHA_QED(0)/2D0/PI C X0 = 2D0*DSQRT(AME2/AMU2)/(1D0+AME2/AMU2) C PARAM(8) = X0 C C ALPHA = ALPHA_QED(0) ---> ALPHA_QED(AMU2) = ALPHAM C DALP = ALPHA/3D0/PI*DLOG(AMU2/AME2) C ALPHAM = ALPHA/(1D0-ALPHA/3D0/PI*DLOG(AMU2/AME2)) C & + ALPHA**3/4D0/PI**2*DLOG(AMU2/AME2) C.... Running alpha_QED should not be used ! C C The best approximation for photonic RC: C FFF = FFM0 + FFM1 + FFL2 + FNL2 + FFL3 + FAHE C GGG = GGM0 + GGM1 + GGL2 + GNL2 + GGL3 + GAHE C C The best approximation for photonic + pair RC: C FFF = FFF + FFVP(X,PARAM) + FFSP(X,ALDP,PARAM) C GGG = GGG + GGVP(X,PARAM) + GGSP(X,ALDP,PARAM) C Note extra parameter in soft-pair contribution: C ALDP = \ln(\Delta), where \Delta is the maximal energy C fraction of the soft pair (E_pair < \Delta m_\mu/2) C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FFF0(X,PARAM) C --------------------------------------- C.... BORN LEVEL F-FUNCTION WITH M_E ---> 0 IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 RHOM,FIS,WAC FFF0 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN RHOM = PARAM(1) C.... EQ.(31) FROM KUNO&OKADA (M_E OMITTED) FIS = X*(1D0-X) + 2D0/9D0*RHOM*( 4D0*X**2 & - 3D0*X ) WAC = X*6D0 FFF0 = FIS*WAC RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GGG0(X,PARAM) C --------------------------------------- C.... BORN LEVEL F-FUNCTION WITH M_E ---> 0 IMPLICIT NONE REAL*8 X, PARAM(11) REAL*8 CHIM, DELTAM,FAS,WAC GGG0 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN CHIM = PARAM(3) DELTAM= PARAM(4) C.... EQ.(32) FROM KUNO&OKADA (M_E OMITTED) FAS = 1D0/3D0*CHIM*X*( & + 1D0 - X + 2D0/3D0*DELTAM*(4D0*X - 3D0) ) WAC = X*6D0 GGG0 = - FAS*WAC RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FFM0(X,PARAM) C --------------------------------------- C.... BORN LEVEL F-FUNCTION WITH EXACT MASS DEPENDENCE IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 AEM2,RHOM,ETAM,X0,FIS,WAC FFM0 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN AEM2 = PARAM(5) RHOM = PARAM(1) ETAM = PARAM(2) X0 = PARAM(8) C.... EQ.(31) FROM KUNO&OKADA FIS = X*(1D0-X) + 2D0/9D0*RHOM*( 4D0*X**2 & - 3D0*X - X0**2 ) + ETAM*X0*(1D0-X) WAC = (1D0+AEM2)**4*DSQRT(X**2-X0**2)*6D0 FFM0 = FIS*WAC RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GGM0(X,PARAM) C --------------------------------------- C.... BORN LEVEL G-FUNCTION WITH EXACT MASS DEPENDENCE IMPLICIT NONE REAL*8 x,PARAM(11) REAL*8 AEM2,CHIM,DELTAM,X0,FAS,WAC GGM0 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN AEM2 = PARAM(5) CHIM = PARAM(3) DELTAM= PARAM(4) X0 = PARAM(8) C.... EQ.(32) FROM KUNO&OKADA FAS = 1D0/3D0*CHIM*DSQRT(X**2-X0**2)*( & + 1D0 - X + 2D0/3D0*DELTAM*(4D0*X - 3D0 & + DSQRT(1D0-X0**2) - 1D0 ) ) WAC = (1D0+AEM2)**4*DSQRT(X**2-X0**2)*6D0 GGM0 = - FAS*WAC RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FFL1(X,PARAM) C --------------------------------------- C.... O(ALPHA*L) F-FUNCTION (M_E ---> 0) IMPLICIT NONE REAL*8 X,PARAM(11) REAl*8 DZ2,ALME,RXL,FLKS,FCTR PARAMETER(DZ2 = 1.64493406684823D0) FFL1 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) C.... LL PART OF EQ.(2.6) FROM K.-S. RXL = 0D0 & + ALME/2D0*( 3D0/2D0 + 2D0*DLOG((1D0-X)/X) ) C.... LL PART OF EQ.(2.4) FROM K.-S. FLKS = (6D0-4D0*X)*RXL & + (1D0-X)/3D0/X**2*( & ( 5D0 + 17D0*X - 34D0*X**2 )*ALME/2D0 & ) FCTR = (1D0-(PARAM(8)/X)**2)*(1D0+PARAM(5))**4 FFL1 = PARAM(7)*FLKS*X**2*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GGL1(X,PARAM) C --------------------------------------- C.... O(ALPHA*L) G-FUNCTION (M_E ---> 0) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,ALME,RXL,GLKS,FCTR PARAMETER(DZ2 = 1.64493406684823D0) GGL1 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) C.... LL PART OF EQ.(2.6) FROM K.-S. RXL = 0D0 & + ALME/2D0*( 3D0/2D0 + 2D0*DLOG((1D0-X)/X) ) C.... LL PART OF EQ.(2.5) FROM K.-S. GLKS = (2D0-4D0*X)*RXL & - (1D0-X)/3D0/X**2*( & ( 1D0 + X + 34D0*X**2 )*ALME/2D0 & ) FCTR = (1D0-(PARAM(8)/X)**2)*(1D0+PARAM(5))**4 GGL1 = PARAM(7)*GLKS*X**2*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FFF1(X,PARAM) C --------------------------------------- C.... O(ALPHA) F-FUNCTION (M_E ---> 0) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,ALME,ALME2,ALX,RX,FFKS,FCTR PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE EXTERNAL SPENCE FFF1 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) ALX = DLOG(X) ALME2 = ALME/2D0 C.... EQ.(2.6) FROM K.-S. RX = 0D0 & + 2D0*SPENCE(X) - 2D0*DZ2 - 2D0 & + ALME2*( 3D0/2D0 + 2D0*DLOG((1D0-X)/X) ) & - ALX*( 2D0*ALX - 1D0 ) & + DLOG(1D0-X)*( 3D0*ALX - 1D0 - 1D0/X ) C.... EQ.(2.4) FROM K.-S. FFKS = (6D0-4D0*X)*RX & + 6D0*(1D0-X)*ALX & + (1D0-X)/3D0/X**2*( & ( 5D0 + 17D0*X - 34D0*X**2 )*( ALME2 + ALX ) & - 22D0*X + 34D0*X**2 & ) FCTR = (1D0-(PARAM(8)/X)**2)*(1D0+PARAM(5))**4 FFF1 = PARAM(7)*FFKS*X**2*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GGG1(X,PARAM) C --------------------------------------- C.... O(ALPHA) G-FUNCTION (M_E ---> 0) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,ALME,ALX,ALME2,RX,FCTR,GGKS PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE EXTERNAL SPENCE GGG1 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) ALX = DLOG(X) ALME2 = ALME/2D0 C.... EQ.(2.6) FROM K.-S. RX = 0D0 & + 2D0*SPENCE(X) - 2D0*DZ2 - 2D0 & + ALME2*( 3D0/2D0 + 2D0*DLOG((1D0-X)/X) ) & - ALX*( 2D0*ALX - 1D0 ) & + DLOG(1D0-X)*( 3D0*ALX - 1D0 - 1D0/X ) C.... EQ.(2.5) FROM K.-S. GGKS = (2D0-4D0*X)*RX & + (2D0-6D0*X)*ALX & - (1D0-X)/3D0/X**2*( & ( 1D0 + X + 34D0*X**2 )*( ALME2 + ALX ) & + 3D0 - 7D0*X - 32D0*X**2 & + 4D0*(1D0-X)**2/X*DLOG(1D0-X) & ) FCTR = (1D0-(PARAM(8)/X)**2)*(1D0+PARAM(5))**4 GGG1 = PARAM(7)*GGKS*X**2*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FFM1(X,PARAM) C --------------------------------------- C.... O(ALPHA) F-FUNCTION WITH EXACT MASS DEPENDENCE IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 AEM2,ALME,Z,BET2,BETA,Q2MU,AAA,F0,FBORN,F1ZZ REAL*8 SPENCE EXTERNAL SPENCE FFM1 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN AEM2 = PARAM(5) ALME = PARAM(6) Z = X*(1D0+AEM2) BET2 = 1D0-4D0*AEM2/Z**2 IF(BET2.LE.0D0) THEN BETA = 0D0 ELSE BETA = DSQRT(BET2) ENDIF C.... BORN-LEVEL F-FUNCTION (WITH Z**2, WITHOUT BETA) F0 = Z**2*(3D0-2D0*Z) + Z**3/4D0*(3D0*Z-4D0)*(1D0-BET2) C.... O(ALPHA) F-FUNCTION (WITHOUT Z**2, WITHOUT BETA) Q2MU = 1D0 - Z + AEM2 C.... EQ.(13) FROM A.A. AAA = ALME*( DLOG(Q2MU) - DLOG(Z) & + DLOG((1D0+BETA)/2D0/BETA) & + DLOG((2D0-Z*(1D0-BETA))/2D0/BETA) ) & + ( DLOG(Q2MU) - 2D0*DLOG(Z) + 2D0*DLOG((1D0+BETA)/2D0) & + 4D0*DLOG((2D0-Z*(1D0-BETA))/2D0/BETA) )*( DLOG(Z) & + DLOG((1D0+BETA)/2D0) ) & + 2D0*SPENCE((1D0-BETA)*(2D0-Z*(1D0+BETA)) & /(1D0+BETA)/(2D0-Z*(1D0-BETA))) & - 2D0*SPENCE((2D0-Z*(1D0+BETA))/(2D0-Z*(1D0-BETA))) C.... EQ.(1) FROM A.A. FBORN = 3D0 - 2D0*Z + Z/4D0*( 3D0*Z - 4D0 )*(1D0-BETA**2) C.... EQ.(12) FROM A.A. F1ZZ = FBORN*( 2D0/BETA*AAA & + ( Z**2*(1D0-BETA**2) - 4D0*(1D0+Z*BETA) )/2D0/Z/BETA & *DLOG(Q2MU) & + ( 4D0 - Z**2*(1D0-BETA**2) )/Z/BETA & *DLOG((2D0-Z*(1D0-BETA))/2D0) & ) & + ( ALME + 2D0*DLOG(Z) + 2D0*DLOG((1D0+BETA)/2D0) )/BETA & *( 5D0*Z**4/384D0*(1D0-BETA**2)**3 - Z**3/4D0*(1D0-BETA**2)**2 & + 3D0*Z**2/32D0*( 3D0 - 12D0*BETA + BETA**2 )*(1D0-BETA**2) & + Z*( 2D0/3D0 + 2D0*BETA + (1D0-BETA**2)*(3D0/2D0+BETA) ) & + 1D0/8D0*( - 20D0 - 12D0*BETA - 19D0*(1D0-BETA**2) ) & + 2D0/Z + 5D0/6D0/Z**2 & ) & + ( DLOG(Z) + DLOG((1D0+BETA)/2D0) )*( & + 9D0/4D0*Z**2*(1D0-BETA**2) & + 2D0*Z*(BETA**2-3D0) + 3D0 & ) & + FBORN*( - 11D0/18D0*Z*(1D0-BETA**2) + 22D0/27D0*BETA**2 & - 2D0/9D0 ) & + Z*( - 22D0/27D0*BETA**4 + BETA**2/2D0 - 11D0/6D0 ) & + 22D0/9D0*(3D0-BETA**2) - 22D0/3D0/Z FFM1 = PARAM(7)*F1ZZ*Z**2*BETA*(1D0+PARAM(5)) RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GGM1(X,PARAM) C --------------------------------------- C.... O(ALPHA) G-FUNCTION WITH EXACT MASS DEPENDENCE IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 AEM2,Z,BET2,BETA,G0,Q2MU,ALME,AAA,GBORN,G1ZZ REAL*8 SPENCE EXTERNAL SPENCE GGM1 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN AEM2 = PARAM(5) Z = X*(1D0+AEM2) BET2 = 1D0-4D0*AEM2/Z**2 IF(BET2.LE.0D0) THEN BETA = 0D0 ELSE BETA = DSQRT(BET2) ENDIF C.... BORN-LEVEL F-FUNCTION (WITH Z**2, WITHOUT BETA) G0 = Z**2*(1D0-2D0*Z)*BETA & + 3D0*Z**4/4D0*(1D0-BET2)*BETA C.... O(ALPHA) G-FUNCTION (WITHOUT Z**2, WITHOUT BETA) Q2MU = 1D0 - Z + AEM2 ALME = PARAM(6) C.... EQ.(13) FROM A.A. AAA = ALME*( DLOG(Q2MU) - DLOG(Z) & + DLOG((1D0+BETA)/2D0/BETA) & + DLOG((2D0-Z*(1D0-BETA))/2D0/BETA) ) & + ( DLOG(Q2MU) - 2D0*DLOG(Z) + 2D0*DLOG((1D0+BETA)/2D0) & + 4D0*DLOG((2D0-Z*(1D0-BETA))/2D0/BETA) )*( DLOG(Z) & + DLOG((1D0+BETA)/2D0) ) & + 2D0*SPENCE((1D0-BETA)*(2D0-Z*(1D0+BETA)) & /(1D0+BETA)/(2D0-Z*(1D0-BETA))) & - 2D0*SPENCE((2D0-Z*(1D0+BETA))/(2D0-Z*(1D0-BETA))) C.... EQ.(1) FROM A.A. GBORN = (1D0-2D0*Z)*BETA + 3D0/4D0*Z**2*(1D0-BETA**2)*BETA C.... EQ.(14) FROM A.A. G1ZZ = GBORN*( 2D0/BETA*AAA & - 4D0*DLOG((2D0-Z*(1D0-BETA))/2D0) & ) & + ( ALME + 2D0*DLOG(Z) + 2D0*DLOG((1D0+BETA)/2D0) )/BETA**2 & *( 5D0*Z**4/384D0*(1D0-BETA**2)**3 & + Z**3/8D0*(1D0-BETA**2)**2*(1D0-3D0*BETA**2) & + 3D0*Z**2/32D0*(1D0-BETA**2)*( - 11D0 + 15D0*BETA**2 & - 12D0*BETA**3 ) & + Z*( 2D0/3D0 + 2D0*BETA + (1D0-BETA**2)*( BETA**2/2D0 & - 2D0*BETA + 3D0/2D0 ) ) & - 7D0/2D0 - BETA/2D0 & + (1D0-BETA**2)*( 17D0/8D0 + BETA/2D0 ) - 1D0/6D0/Z**2 ) & + ( DLOG(Z) + DLOG((1D0+BETA)/2D0) )*BETA*( & + 9D0/4D0*Z**2*(1D0-BETA**2) - 4D0*Z + 1D0 ) & + 1D0/BETA**2*( DLOG(Q2MU) & - 2D0*DLOG((2D0-Z*(1D0-BETA))/2D0) )*( & - Z**3/48D0*(1D0-BETA**2)**2*(1D0-19D0*BETA**2) & + Z**2*(1D0-BETA**2)*( - 3D0/2D0*BETA**3 & - 5D0/4D0*BETA**2 + 1D0/4D0 ) & + Z*( 4D0*BETA + (1D0-BETA**2)*( - 3D0/4D0*BETA**2 & - 4D0*BETA - 5D0/4D0 ) ) & + 16D0/3D0 - 2D0*BETA - (1D0-BETA**2)*2D0*(1D0-BETA) & + 1D0/Z*( - 6D0 + 1D0 - BETA**2 ) + 4D0/Z**2 & - 4D0/3D0/Z**3 ) & + GBORN*( - 5D0/144D0/BETA**2*Z**2*(1D0-BETA**2)**2 & - 10D0/27D0/BETA**2*Z*(1D0-BETA**2) - 55D0/54D0 & + 203D0/162D0/BETA**2 ) & + Z/81D0*( 17D0/BETA - 195D0*BETA ) & - 1D0/324D0*( 595D0/BETA - 1923D0*BETA ) & + 10D0/3D0/Z/BETA - 1D0/Z**2/BETA GGM1 = PARAM(7)*G1ZZ*Z**2*BETA*(1D0+PARAM(5)) RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FFL2(X,PARAM) C --------------------------------------- C.... O(ALPHA^2L^2) F-FUNCTION (M_E ---> 0) C (PURE PHOTONIC CORRECTIONS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,ALME,PHI,FF2,FCTR PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE EXTERNAL SPENCE FFL2 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) C.... EQ.(14) FROM ACG PHI = DLOG(X)**2/2D0 + DLOG(1D0-X)**2 & - 2D0*DLOG(X)*DLOG(1D0-X) - SPENCE(1D0-X) - DZ2 FF2 = 1D0/2D0*ALME**2*( & + 4D0*X**2*( 3D0 - 2D0*X )*PHI & + ( 10D0/3D0 + 8D0*X - 16D0*X**2 + 32D0/3D0*X**3 )*DLOG(1D0-X) & + ( - 5D0/6D0 - 2D0*X + 8D0*X**2 - 32D0/3D0*X**3 )*DLOG(X) & + 11D0/36D0 + 17D0/6D0*X + 8D0/3D0*X**2 - 32D0/9D0*X**3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) FFL2 = PARAM(7)**2*FF2*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GGL2(X,PARAM) C --------------------------------------- C.... O(ALPHA^2L^2) G-FUNCTION (M_E ---> 0) C (PURE PHOTONIC CORRECTIONS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2, ALME,PHI,GG2,FCTR PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE EXTERNAL SPENCE GGL2 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) C.... EQ.(15) FROM ACG PHI = DLOG(X)**2/2D0 + DLOG(1D0-X)**2 & - 2D0*DLOG(X)*DLOG(1D0-X) - SPENCE(1D0-X) - DZ2 GG2 = 1D0/2D0*ALME**2*( & + 4D0*X**2*( 1D0 - 2D0*X )*PHI & + ( - 2D0/3D0 - 16D0*X**2 + 32D0/3D0*X**3 )*DLOG(1D0-X) & + ( 1D0/6D0 + 8D0*X**2 - 32D0/3D0*X**3 )*DLOG(X) & - 7D0/36D0 - 7D0/6D0*X + 8D0/3D0*X**2 - 32D0/9D0*X**3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) GGL2 = PARAM(7)**2*GG2*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FNL2(X,PARAM) C --------------------------------------- C.... O(ALPHA^2L^1) F-FUNCTION (M_E ---> 0) C (PURE PHOTONIC CORRECTIONS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,DL1X,SP1X,TRIX,S12X,FN2,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE,TRILOG,S12 EXTERNAL SPENCE,TRILOG,S12 FNL2 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) DL1X = DLOG(1D0-X) SP1X = SPENCE(1D0-X) TRIX = TRILOG(1D0-X) S12X = S12(1D0-X) C.... EQ.(26) FROM AM C SPX = SPENCE(X) C TRX = TRILOG(X) C S1X = S12(X) C FN2AM = ALME*( 0D0 C & + 2D0*X**2*(3D0-2D0*X)*( - 2D0*TRX - 2D0*S1X C & + 2D0*SPX*DL1X + 2D0*SPX*DLX + 5D0*DLX*DL1X**2 C & - 5D0*DLX**2*DL1X + 2D0*DLX**3 - 2D0*DZ2*DL1X C & - 2D0*DZ2*DLX + 7D0*DZ3 ) C & + SPX*( 10D0/3D0 + 14D0*X - 40D0*X**2 + 92D0/3D0*X**3 ) C & + DLX*DL1X * ( 25D0/3D0 + 32D0*X - 54D0*X**2 + 92D0/3D0*X**3 ) C & + DL1X**2 * ( - 12D0*X - 4D0*X**2 + 8D0*X**3 ) C & + DLX**2 * ( - 25D0/12D0 - 5D0*X + 22D0*X**2 C & - 70D0/3D0*X**3 ) C & + DL1X * ( - 17D0/3D0 - 53D0/3D0*X + 64D0/3D0*X**2 C & - 12D0*X**3 ) C & + DLX * ( - 3D0/4D0 + 37D0/6D0*X + 4D0/3D0*X**2 C & + 44D0/9D0*X**3 ) C & + DZ2*( - 10D0/3D0 - 2D0*X + 35D0*X**2 - 98D0/3D0*X**3 ) C & + 211D0/216D0 - 287D0/12D0*X + 83D0/3D0*X**2 C & - 559D0/54D0*X**3 C & ) C.... EQUIVALENT TO: FN2 = ALME*( 0D0 & + TRIX * ( 12D0*X**2 - 8D0*X**3 ) & + S12X * ( 12D0*X**2 - 8D0*X**3 ) & + DZ2*DLX * ( - 12D0*X**2 + 8D0*X**3 ) & + DZ2 * ( 12D0*X - 5D0*X**2 - 2D0*X**3 ) & + DZ3 * ( 18D0*X**2 - 12D0*X**3 ) & + SP1X*DL1X * ( - 24D0*X**2 + 16D0*X**3 ) & + DLX*DL1X**2 * ( 12D0*X**2 - 8D0*X**3 ) & + DLX**2*DL1X * ( - 36D0*X**2 + 24D0*X**3 ) & + DLX**3 * ( 12D0*X**2 - 8D0*X**3 ) & + SP1X * ( - 10D0/3D0 - 14D0*X + 40*X**2 - 92D0/3D0*X**3 ) & + DLX*DL1X * ( 5D0 + 18D0*X - 14D0*X**2 ) & + DL1X**2 * ( - 12D0*X - 4D0*X**2 + 8D0*X**3 ) & + DLX**2 * ( - 25D0/12D0 - 5D0*X + 22D0*X**2 & - 70D0/3D0*X**3 ) & + DLX * ( - 3D0/4D0 + 37D0/6D0*X + 4D0/3D0*X**2 & + 44D0/9D0*X**3 ) & + DL1X * ( - 17D0/3D0 - 53D0/3D0*X + 64D0/3D0*X**2 & - 12D0*X**3 ) & + 211D0/216D0 - 287D0/12D0*X + 83D0/3D0*X**2 & - 559D0/54D0*X**3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) FNL2 = PARAM(7)**2*FN2*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GNL2(X,PARAM) C --------------------------------------- C.... O(ALPHA^2L^1) G-FUNCTION (M_E ---> 0) C (PURE PHOTONIC CORRECTIONS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,DL1X,SP1X,TRIX,S12X,GN2,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE,TRILOG,S12 EXTERNAL SPENCE,TRILOG,S12 GNL2 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) DL1X = DLOG(1D0-X) SP1X = SPENCE(1D0-X) TRIX = TRILOG(1D0-X) S12X = S12(1D0-X) C.... EQ.(5.3) FROM AA (G_2^{(1,\gamma)}) GN2 = ALME*( 0D0 & + TRIX * ( 4D0*X**2 - 8D0*X**3 ) & + S12X * ( 4D0*X**2 - 8D0*X**3 ) & + SP1X*DL1X * ( - 8D0*X**2 + 16D0*X**3 ) & + DLX*DL1X**2 * ( 4D0*X**2 - 8D0*X**3 ) & + DLX**2*DL1X * ( - 12D0*X**2 + 24D0*X**3 ) & + DLX**3 * ( 4D0*X**2 - 8D0*X**3 ) & + DZ2*DLX * ( - 4D0*X**2 + 8D0*X**3 ) & + DZ3 * ( 6D0*X**2 - 12D0*X**3 ) & + SP1X * ( 14D0/3D0 - 8D0/3D0/X - 6D0*X + 24D0*X**2 & - 92D0/3D0*X**3 ) & + DLX*DL1X * ( - 5D0 + 6D0*X - 86D0/3D0*X**2 ) & + DL1X**2 * ( 8D0 - 8D0/3D0/X - 12D0*X + 20D0/3D0*X**2 & + 8D0*X**3 ) & + DLX**2 * ( 5D0/12D0 + 18D0*X**2 - 70D0/3D0*X**3 ) & + DL1X * ( - 13D0/3D0 + 37D0/3D0*X + 50D0/3D0*X**2 & - 32D0/3D0*X**3 ) & + DLX * ( 25D0/12D0 - 59D0/6D0*X + 6D0*X**2 + 32D0/9D0*X**3 ) & + DZ2 * ( - 8D0 + 8D0/3D0/X + 12D0*X - 29D0/3D0*X**2 & - 2D0*X**3 ) & + 817D0/216D0 - 91D0/12D0*X + 62D0/3D0*X**2 - 607D0/54D0*X**3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) GNL2 = PARAM(7)**2*GN2*FCTR RETURN END C----------------------------------------------------------------------- C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FFL3(X,PARAM) C --------------------------------------- C.... O(ALPHA^3L^3) F-FUNCTION (M_E ---> 0) C (PURE PHOTONIC CORRECTIONS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,PSI,FF3,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE,TRILOG,S12 EXTERNAL SPENCE,TRILOG,S12 FFL3 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) C.... EQ.(XX) FROM A.A. PSI = 3D0*TRILOG(1D0-X) - 2D0*S12(1D0-X) & + DLOG(1D0-X)**3 - DLOG(X)**3/6D0 & + 3D0/2D0*DLOG(X)**2*DLOG(1D0-X) & - 3D0*DLOG(X)*DLOG(1D0-X)**2 & - 3D0*SPENCE(1D0-X)*DLOG(1D0-X) & + 2D0*DZ3 - 3D0*DZ2*( DLOG(1D0-X) - DLOG(X) ) FF3 = 1D0/6D0*(ALME-1D0)**3*( & + 8D0*X**2*(3D0-2D0*X)*PSI & + ( 10D0 + 24D0*X - 48D0*X**2 + 32D0*X**3 )*DLOG(1D0-X)**2 & + ( 5D0/12D0 + X - 8D0*X**2 + 16D0*X**3 )*DLOG(X)**2 & + ( - 5D0 - 12D0*X + 48D0*X**2 - 64D0*X**3 )*DLOG(X)*DLOG(1D0-X) & + ( 5D0 + 12D0*X - 32D0*X**3 )*SPENCE(1D0-X) & + ( - 10D0 - 24D0*X + 48D0*X**2 - 32D0*X**3 )*DZ2 & + ( - 13D0/18D0 - 21D0/2D0*X + 64D0/3D0*X**3 )*DLOG(X) & + ( 11D0/6D0 + 17D0*X + 16D0*X**2 - 64D0/3D0*X**3 )*DLOG(1D0-X) & + 569D0/216D0 + 4D0/3D0*X - 16D0/3D0*X**2 + 128D0/27D0*X**3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) FFL3 = PARAM(7)**3*FF3*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GGL3(X,PARAM) C --------------------------------------- C.... O(ALPHA^3L^3) G-FUNCTION (M_E ---> 0) C (PURE PHOTONIC CORRECTIONS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,PSI,GG3,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE,TRILOG,S12 EXTERNAL SPENCE,TRILOG,S12 GGL3 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) C.... EQ.(XX) FROM A.A. PSI = 3D0*TRILOG(1D0-X) - 2D0*S12(1D0-X) & + DLOG(1D0-X)**3 - DLOG(X)**3/6D0 & + 3D0/2D0*DLOG(X)**2*DLOG(1D0-X) & - 3D0*DLOG(X)*DLOG(1D0-X)**2 & - 3D0*SPENCE(1D0-X)*DLOG(1D0-X) & + 2D0*DZ3 - 3D0*DZ2*( DLOG(1D0-X) - DLOG(X) ) GG3 = 1D0/6D0*(ALME-1D0)**3*( & + 8D0*X**2*(1D0-2D0*X)*PSI & + ( - 2D0 - 48D0*X**2 + 32D0*X**3 )*DLOG(1D0-X)**2 & + ( - 1D0/12D0 - 8D0*X**2 + 16D0*X**3 )*DLOG(X)**2 & + ( 1D0 + 48D0*X**2 - 64D0*X**3 )*DLOG(X)*DLOG(1D0-X) & + ( - 1D0 - 32D0*X**3 )*SPENCE(1D0-X) & + ( 2D0 + 48D0*X**2 - 32D0*X**3 )*DZ2 & + ( 5D0/18D0 + 5D0/2D0*X + 64D0/3D0*X**3 )*DLOG(X) & + ( - 7D0/6D0 - 7D0*X + 16D0*X**2 - 64D0/3D0*X**3 )*DLOG(1D0-X) & - 133D0/216D0 - 13D0/6D0*X - 16D0/3D0*X**2 + 128D0/27D0*X**3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) GGL3 = PARAM(7)**3*GG3*FCTR RETURN END C----------------------------------------------------------------------- C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FNS0(X,PARAM) C --------------------------------------- C.... O(ALPHA^2L^2) F-FUNCTION (NS PAIRS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,DL1X,FNSLL,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) FNS0 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) DL1X = DLOG(1D0-X) C.... EQ.(12,18) FROM ACG FNSLL = 1D0/3D0*ALME**2*( 0D0 & + 5D0/6D0 + 2D0*X - 4D0*X**2 + 8D0/3D0*X**3 & + 2D0*X**2*(3D0-2D0*X)*( DL1X - DLX ) & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) FNS0 = PARAM(7)**2*FNSLL*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GNS0(X,PARAM) C --------------------------------------- C.... O(ALPHA^2L^2) G-FUNCTION (NS PAIRS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,DL1X,SP1X,TRIX,S12X,GNSLL,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE,TRILOG,S12 EXTERNAL SPENCE,TRILOG,S12 GNS0 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) DL1X = DLOG(1D0-X) SP1X = SPENCE(1D0-X) TRIX = TRILOG(1D0-X) S12X = S12(1D0-X) C.... EQ.(13,18) FROM ACG GNSLL = 1D0/3D0*ALME**2*( 0D0 & - 1D0/6D0 - 4D0*X**2 + 8D0/3D0*X**3 & + 2D0*X**2*(1D0-2D0*X)*( DL1X - DLX ) & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) GNS0 = PARAM(7)**2*GNSLL*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FFS0(X,PARAM) C --------------------------------------- C.... O(ALPHA^2L^2) F-FUNCTION (S PAIRS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,FSLL,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) FFS0 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) C.... EQ.(19) FROM ACG FSLL = 1D0/2D0*ALME**2*( 0D0 & + 17D0/9D0 + 2D0/3D0/X + 3D0*X - 14D0/3D0*X**2 & - 8D0/9D0*X**3 + ( 5D0/3D0 + 4D0*X + 4D0*X**2 )*DLX & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) FFS0 = PARAM(7)**2*FSLL*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GGS0(X,PARAM) C --------------------------------------- C.... O(ALPHA^2L^2) G-FUNCTION (S PAIRS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,GSLL,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) GGS0 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) C.... EQ.(19) FROM ACG GSLL = 1D0/2D0*ALME**2*( 0D0 & - 1D0/9D0 - 2D0/9D0/X + X + 2D0/9D0*X**2 & - 8D0/9D0*X**3 + ( - 1D0/3D0 + 4D0/3D0*X**2 )*DLX & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) GGS0 = PARAM(7)**2*GSLL*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FNS1(X,PARAM) C --------------------------------------- C.... O(ALPHA^2L) F-FUNCTION (NS PAIRS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,DL1X,SP1X,TRIX,S12X,FNSNL,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE,TRILOG,S12 EXTERNAL SPENCE,TRILOG,S12 FNS1 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) DL1X = DLOG(1D0-X) SP1X = SPENCE(1D0-X) TRIX = TRILOG(1D0-X) S12X = S12(1D0-X) C.... EQ.(30) FROM AM FNSNL = ALME*( 0D0 & + 2D0*X**2*(3D0-2D0*X)*( - 2D0*SP1X - 2D0/3D0*DLX*DL1X & + 2D0/3D0*DL1X**2 - DLX**2 - 2D0/3D0*DZ2 ) & + DL1X*( 10D0/9D0 - 4D0/3D0*X - 46D0/3D0*X**2 + 12D0*X**3 ) & + DLX*( 5D0/9D0 + 4D0/3D0*X + 8D0*X**2 - 76D0/9D0*X**3 ) & - 11D0/6D0 - 19D0/3D0*X + 100D0/9D0*X**2 - 64D0/9D0*X**3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) FNS1 = PARAM(7)**2*FNSNL*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GNS1(X,PARAM) C --------------------------------------- C.... O(ALPHA^2L) G-FUNCTION (NS PAIRS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,DL1X,SP1X,TRIX,S12X,GNSNL,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE,TRILOG,S12 EXTERNAL SPENCE,TRILOG,S12 GNS1 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) DL1X = DLOG(1D0-X) SP1X = SPENCE(1D0-X) TRIX = TRILOG(1D0-X) S12X = S12(1D0-X) C.... EQ.(XX) FROM ABA GNSNL = ALME*( 0D0 & + 4D0*X**2*(1D0-2D0*X)*( - SP1X - 1D0/3D0*DLX*DL1X & + 1D0/3D0*DL1X**2 - DLX**2/2D0 - 1D0/3D0*DZ2 ) & + DL1X*( 22D0/9D0 - 8D0/9D0/X - 4D0*X - 6D0*X**2 + 12D0*X**3 ) & + DLX*( - 1D0/9D0 + 8D0/9D0*X**2 - 76D0/9D0*X**3 ) & - 7D0/18D0 + 5D0/3D0*X + 86D0/9D0*X**2 - 20D0/3D0*X**3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) GNS1 = PARAM(7)**2*GNSNL*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FFS1(X,PARAM) C --------------------------------------- C.... O(ALPHA^2L) F-FUNCTION (S PAIRS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,DL1X,SP1X,TRIX,S12X,FFSNL,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE,TRILOG,S12 EXTERNAL SPENCE,TRILOG,S12 FFS1 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) DL1X = DLOG(1D0-X) SP1X = SPENCE(1D0-X) TRIX = TRILOG(1D0-X) S12X = S12(1D0-X) C.... EQ.(33) FROM AM FFSNL = ALME*( 0D0 & + ( SP1X + DLX*DL1X ) * ( 5D0/3D0 + 4D0*X + 4D0*X**2 ) & + DLX**2*( 5D0/2D0 + 6D0*X + 4D0*X**2 ) & + DL1X*( 17D0/9D0 + 2D0/3D0/X + 3D0*X - 14D0/3D0*X**2 & - 8D0/9D0*X**3 ) & + DLX*( 8D0/9D0 + 4D0/3D0/X - 5D0/6D0*X - 19D0/3D0*X**2 ) & - 67D0/9D0 - 1D0/3D0/X + 43D0/18D0*X + 77D0/18D0*X**2 & + 10D0/9D0*X**3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) FFS1 = PARAM(7)**2*FFSNL*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GGS1(X,PARAM) C --------------------------------------- C.... O(ALPHA^2L) G-FUNCTION (S PAIRS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,DL1X,SP1X,TRIX,S12X,GGSNL,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE,TRILOG,S12 EXTERNAL SPENCE,TRILOG,S12 GGS1 = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) DL1X = DLOG(1D0-X) SP1X = SPENCE(1D0-X) TRIX = TRILOG(1D0-X) S12X = S12(1D0-X) C.... EQ.(XX) FROM ABA GGSNL = ALME*( 0D0 & + ( 4D0/3D0*X**2 - 1D0/3D0 )*( SP1X + DLX*DL1X ) & + DLX**2*( - 1D0/2D0 + 4D0/3D0*X**2 ) & + DL1X*( - 1D0/9D0 - 2D0/9D0/X + X + 2D0/9D0*X**2 & - 8D0/9D0*X**3 ) & + DLX*( 5D0/9D0 - 4D0/9D0/X + 5D0/2D0*X + 5D0/9D0*X**2 ) & + 4D0/3D0 + 1D0/3D0/X - 7D0/18D0*X - 43D0/18D0*X**2 & + 10D0/9D0*X**3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) GGS1 = PARAM(7)**2*GGSNL*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FFFI(X,PARAM) C --------------------------------------- C.... O(ALPHA^2L) F-FUNCTION (S*NS PAIR INTERFERENCE ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,DL1X,SP1X,TRIX,S12X,FFSNS,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE,TRILOG,S12 EXTERNAL SPENCE,TRILOG,S12 FFFI = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) DL1X = DLOG(1D0-X) SP1X = SPENCE(1D0-X) TRIX = TRILOG(1D0-X) S12X = S12(1D0-X) C.... EQ.(35) FROM AM FFSNS = ALME*( 0D0 & + 4D0*X**2*(3D0-2D0*X)*( TRIX - 2D0*S12X - DLX*SP1X ) & + SP1X*( 5D0/3D0 + 4D0*X - 26D0*X**2 + 52D0/3D0*X**3 ) & + DLX**2 * ( - 9D0*X**2 + 26D0/3D0*X**3 ) & + DLX * ( - 5D0/3D0 - 5D0/3D0*X - 28D0/3D0*X**2 ) & - 62D0/9D0 + 41D0/3D0*X - 55D0/3D0*X**2 + 104D0/9D0*X**3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) FFFI = PARAM(7)**2*FFSNS*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GGGI(X,PARAM) C --------------------------------------- C.... O(ALPHA^2L) G-FUNCTION (S*NS PAIR INTERFERENCE ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,DL1X,SP1X,TRIX,S12X,GGSNS,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE,TRILOG,S12 EXTERNAL SPENCE,TRILOG,S12 GGGI = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) DL1X = DLOG(1D0-X) SP1X = SPENCE(1D0-X) TRIX = TRILOG(1D0-X) S12X = S12(1D0-X) C.... EQ.(XX) FROM ABA GGSNS = ALME*( 0D0 & + 4D0*X**2*(1D0-2D0*X)*( TRIX - 2D0*S12X - DLX*SP1X ) & + SP1X*( - 1D0/3D0 - 14D0*X**2 + 52D0/3D0*X**3 ) & + DLX**2 * ( - 3D0*X**2 + 26D0/3D0*X**3 ) & + DLX * ( 1D0/3D0 + 1D0/3D0*X - 28D0/3D0*X**2 ) & + 10D0/9D0 - 1D0/3D0*X - 37D0/3D0*X**2 + 104D0/9D0*X**3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) GGGI = PARAM(7)**2*GGSNS*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION F3NS(X,PARAM) C --------------------------------------- C.... O(ALPHA^3L^3) F-FUNCTION (NS PAIRS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,DL1X,SP1X,F3N,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE,TRILOG,S12 EXTERNAL SPENCE,TRILOG,S12 F3NS = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) DL1X = DLOG(1D0-X) SP1X = SPENCE(1D0-X) C.... EQ.(5.10) FROM ABA F3N = 1D0/6D0*ALME**3*( 0D0 & + 8D0*X**2*(3D0-2D0*X)*( DLX**2/2D0 + DL1X**2 & - 2D0*DLX*DL1X - SP1X - DZ2 ) & + ( 20D0/3D0 + 16D0*X - 80D0/3D0*X**2 + 160D0/9D0*X**3 )*DL1X & + ( - 5D0/3D0 - 4D0*X + 32D0/3D0*X**2 - 160D0/9D0*X**3 )*DLX & + 73D0/54D0 + 67D0/9D0*X + 16D0/9D0*X**2 - 128D0/27D0*X**3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) F3NS = PARAM(7)**3*F3N*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION G3NS(X,PARAM) C --------------------------------------- C.... O(ALPHA^3L^3) G-FUNCTION (NS PAIRS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,DL1X,SP1X,G3N,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE,TRILOG,S12 EXTERNAL SPENCE,TRILOG,S12 G3NS = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) DL1X = DLOG(1D0-X) SP1X = SPENCE(1D0-X) C.... EQ.(5.11) FROM ABA G3N = 1D0/6D0*ALME**3*( 0D0 & + 8D0*X**2*(1D0-2D0*X)*( DLX**2/2D0 + DL1X**2 & - 2D0*DLX*DL1X - SP1X - DZ2 ) & + ( - 4D0/3D0 - 272D0/9D0*X**2 + 160D0/9D0*X**3 )*DL1X & + ( 1D0/3D0 + 128D0/9D0*X**2 - 160D0/9D0*X**3 )*DLX & - 29D0/54D0 - 7D0/3D0*X + 16D0/9D0*X**2 - 128D0/27D0*X**3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) G3NS = PARAM(7)**3*G3N*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION F3SS(X,PARAM) C --------------------------------------- C.... O(ALPHA^3L^3) F-FUNCTION (S PAIRS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,DL1X,SP1X,F3S,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE,TRILOG,S12 EXTERNAL SPENCE,TRILOG,S12 F3SS = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) DL1X = DLOG(1D0-X) SP1X = SPENCE(1D0-X) C.... EQ.(5.12) FROM ABA F3S = 1D0/6D0*ALME**3*( 0D0 & + ( 5D0/3D0 + 4D0*X + 4D0*X**2 )*( 4D0*SP1X & + 4D0*DLX*DL1X - DLX**2 ) - 4D0*X**2*DLX**2 & + ( 68D0/9D0 + 8D0/3D0/X + 12D0*X - 56D0/3D0*X**2 & - 32D0/9D0*X**3 )*DL1X & + ( - 29D0/9D0 - 14D0/3D0*X + 16D0*X**2 + 32D0/9D0*X**3 )*DLX & - 287D0/27D0 - 4D0/9D0/X - 13D0/9D0*X + 86D0/9D0*X**2 & + 80D0/27D0*X**3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) F3SS = PARAM(7)**3*F3S*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION G3SS(X,PARAM) C --------------------------------------- C.... O(ALPHA^3L^3) G-FUNCTION (S PAIRS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,DL1X,SP1X,G3S,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE,TRILOG,S12 EXTERNAL SPENCE,TRILOG,S12 G3SS = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) DL1X = DLOG(1D0-X) SP1X = SPENCE(1D0-X) C.... EQ.(5.13) FROM ABA G3S = 1D0/6D0*ALME**3*( 0D0 & + ( 4D0/3D0*X**2 - 1D0/3D0 )*( 4D0*SP1X & + 4D0*DLX*DL1X - DLX**2 ) - 4D0/3D0*X**2*DLX**2 & + ( - 4D0/9D0 - 8D0/9D0/X + 4D0*X + 8D0/9D0*X**2 )*DL1X & + ( 1D0/9D0 - 2D0*X - 16D0/9D0*X**2 + 32D0/9D0*X**3 )*DLX & + 31D0/27D0 + 4D0/27D0/X - 35D0/9D0*X - 10D0/27D0*X**2 & + 80D0/27D0*X**3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) G3SS = PARAM(7)**3*G3S*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FFSP(X,ALDP,PARAM) C -------------------------------------------- C.... O(ALPHA^2) F-FUNCTION (SOFT PAIRS ONLY) IMPLICIT NONE REAL*8 X,ALDP,PARAM(11) REAL*8 DZ2,DZ3,ALME,AEM2,AEM,XXP,F0,ALA,FFFSP,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) FFSP = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) AEM2 = PARAM(5) AEM = DSQRT(AEM2) C.... MAXIMAL ELECTRON ENERGY FRACTION FOR FIVE-LEPTONIC DECAY XXP = ( 1D0 - 4D0*AEM + 5D0*AEM2 ) & / ( 1D0 - 2D0*AEM ) / ( 1D0 + AEM2 ) IF(X.GE.XXP) RETURN C.... EQ.(15) FROM DRAFT F0 = X**2*(3D0-2D0*X) ALA = ALME + 2D0*ALDP FFFSP = F0*4D0/3D0*( & + ALA**3/12D0 & - 2D0/3D0*ALA**2 & + ALA*( 61D0/18D0 - DZ2 ) & - 223D0/27D0 + 8D0/3D0*DZ2 + 2D0*DZ3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) FFSP = PARAM(7)**2*FFFSP*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GGSP(X,ALDP,PARAM) C -------------------------------------------- C.... O(ALPHA^2) G-FUNCTION (SOFT PAIRS ONLY) IMPLICIT NONE REAL*8 X,ALDP,PARAM(11) REAL*8 DZ2,DZ3,ALME,AEM2,AEM,XXP,G0,ALA,GGGSP,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) GGSP = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) AEM2 = PARAM(5) AEM = DSQRT(AEM2) C.... MAXIMAL ELECTRON ENERGY FRACTION FOR FIVE-LEPTONIC DECAY XXP = ( 1D0 - 4D0*AEM + 5D0*AEM2 ) & / ( 1D0 - 2D0*AEM ) / ( 1D0 + AEM2 ) IF(X.GE.XXP) RETURN C.... EQ.(15) FROM DRAFT G0 = X**2*(1D0-2D0*X) ALA = ALME + 2D0*ALDP GGGSP = G0*4D0/3D0*( & + ALA**3/12D0 & - 2D0/3D0*ALA**2 & + ALA*( 61D0/18D0 - DZ2 ) & - 223D0/27D0 + 8D0/3D0*DZ2 + 2D0*DZ3 & ) FCTR = DSQRT(1D0-(PARAM(8)/X)**2) GGSP = PARAM(7)**2*GGGSP*FCTR RETURN END C--------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FFVP(X,PARAM) C --------------------------------------- C.... O(ALPHA^2) F-FUNCTION (VIRTUAL PAIRS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,DL1X,SP1X,TRIX,S12X,F0,WX,FFFVP,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE,TRILOG,S12 EXTERNAL SPENCE,TRILOG,S12 FFVP = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) DL1X = DLOG(1D0-X) SP1X = SPENCE(1D0-X) TRIX = TRILOG(1D0-X) S12X = S12(1D0-X) C.... EQ.(??) FROM DRAFT F0 = X**2*(3D0-2D0*X) WX = - ALME**3/9D0 & + ( 25D0/18D0 - 2D0/3D0*DLX )*ALME**2 & + ( - 397D0/54D0 - 4D0/3D0*DZ2 + 38D0/9D0*DLX & - 4D0/3D0*DLX**2 - 4D0/3D0*SP1X )*ALME & + 517D0/27D0 - 8D0/3D0*DZ2*DLX + 22D0/9D0*DZ2 & + 4D0/3D0*DZ3 - 8D0/3D0*DLX*SP1X - 265D0/27D0*DLX & + 38D0/9D0*DLX**2 - 8D0/9D0*DLX**3 & + 38D0/9D0*SP1X - 8D0/3D0*S12X + 4D0/3D0*TRIX FFFVP = F0*WX & - 2D0*X**2*DLX*ALME & - 2D0*X**2*DLX**2 & - 2D0*X**2*SP1X - 2D0/3D0/(1D0-X)*DLX & + 2D0/3D0*X*DLX + 7D0*X**2*DLX + 2D0/3D0*DLX FCTR = DSQRT(1D0-(PARAM(8)/X)**2) FFVP = PARAM(7)**2*FFFVP*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GGVP(X,PARAM) C --------------------------------------- C.... O(ALPHA^2) G-FUNCTION (VIRTUAL PAIRS ONLY) IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 DZ2,DZ3,ALME,DLX,DL1X,SP1X,TRIX,S12X,G0,WX,GGGVP,FCTR PARAMETER(DZ3 = 1.20205690315959D0) PARAMETER(DZ2 = 1.64493406684823D0) REAL*8 SPENCE,TRILOG,S12 EXTERNAL SPENCE,TRILOG,S12 GGVP = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ALME = PARAM(6) DLX = DLOG(X) DL1X = DLOG(1D0-X) SP1X = SPENCE(1D0-X) TRIX = TRILOG(1D0-X) S12X = S12(1D0-X) C.... EQ.(15) FROM DRAFT G0 = X**2*(1D0-2D0*X) WX = - ALME**3/9D0 & + ( 25D0/18D0 - 2D0/3D0*DLX )*ALME**2 & + ( - 397D0/54D0 - 4D0/3D0*DZ2 + 38D0/9D0*DLX & - 4D0/3D0*DLX**2 - 4D0/3D0*SP1X )*ALME & + 517D0/27D0 - 8D0/3D0*DZ2*DLX + 22D0/9D0*DZ2 & + 4D0/3D0*DZ3 - 8D0/3D0*DLX*SP1X - 265D0/27D0*DLX & + 38D0/9D0*DLX**2 - 8D0/9D0*DLX**3 & + 38D0/9D0*SP1X - 8D0/3D0*S12X + 4D0/3D0*TRIX GGGVP = G0*WX & - 2D0/3D0*X**2*DLX*ALME & - 2D0/3D0*X**2*DLX**2 & - 2D0/3D0*X**2*SP1X + 2D0/3D0/(1D0-X)*DLX & - 2D0/3D0*X*DLX + 13D0/9D0*X**2*DLX - 2D0/3D0*DLX FCTR = DSQRT(1D0-(PARAM(8)/X)**2) GGVP = PARAM(7)**2*GGGVP*FCTR RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FAHE(X,PARAM) C --------------------------------------- C.... AD HOC EXPONENTIATION IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 ARG REAL*8 FFF0 EXTERNAL FFF0 FAHE = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ARG = 2D0*PARAM(7)*(PARAM(6)-2D0)*DLOG(1D0-X) FAHE = FFF0(X,PARAM)*( DEXP(ARG) & - 1D0 - ARG & - (2D0*PARAM(7)*DLOG(1D0-X))**2/2D0*( PARAM(6)**2 & - 4D0*PARAM(6) ) & - (2D0*PARAM(7)*DLOG(1D0-X))**3/6D0*PARAM(6)**3 & ) RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GAHE(X,PARAM) C --------------------------------------- C.... AD HOC EXPONENTIATION IMPLICIT NONE REAL*8 X,PARAM(11) REAL*8 ARG REAL*8 GGG0 EXTERNAL GGG0 GAHE = 0D0 IF(X.LE.PARAM(8).OR.X.GE.1D0) RETURN ARG = 2D0*PARAM(7)*(PARAM(6)-2D0)*DLOG(1D0-X) GAHE = GGG0(X,PARAM)*( DEXP(ARG) & - 1D0 - ARG & - (2D0*PARAM(7)*DLOG(1D0-X))**2/2D0*( PARAM(6)**2 & - 4D0*PARAM(6) ) & - (2D0*PARAM(7)*DLOG(1D0-X))**3/6D0*PARAM(6)**3 & ) RETURN END C---------------------------------------------------------------------- C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION DDILOG(X) * ====== ========= ======== ========= C IMPLICIT NONE DOUBLE PRECISION X,Y,T,S,A,PI3,PI6,ZERO,ONE,HALF,MALF,MONE,MTWO DOUBLE PRECISION C(0:18),H,ALFA,B0,B1,B2 INTEGER I C DATA ZERO /0.0D0/, ONE /1.0D0/ DATA HALF /0.5D0/, MALF /-0.5D0/, MONE /-1.0D0/, MTWO /-2.0D0/ DATA PI3 /3.28986 81336 96453D0/, PI6 /1.64493 40668 48226D0/ C DATA C( 0) / 0.42996 69356 08137 0D0/ DATA C( 1) / 0.40975 98753 30771 1D0/ DATA C( 2) /-0.01858 84366 50146 0D0/ DATA C( 3) / 0.00145 75108 40622 7D0/ DATA C( 4) /-0.00014 30418 44423 4D0/ DATA C( 5) / 0.00001 58841 55418 8D0/ DATA C( 6) /-0.00000 19078 49593 9D0/ DATA C( 7) / 0.00000 02419 51808 5D0/ DATA C( 8) /-0.00000 00319 33412 7D0/ DATA C( 9) / 0.00000 00043 45450 6D0/ DATA C(10) /-0.00000 00006 05784 8D0/ DATA C(11) / 0.00000 00000 86121 0D0/ DATA C(12) /-0.00000 00000 12443 3D0/ DATA C(13) / 0.00000 00000 01822 6D0/ DATA C(14) /-0.00000 00000 00270 1D0/ DATA C(15) / 0.00000 00000 00040 4D0/ DATA C(16) /-0.00000 00000 00006 1D0/ DATA C(17) / 0.00000 00000 00000 9D0/ DATA C(18) /-0.00000 00000 00000 1D0/ IF(X .EQ. ONE) THEN DDILOG=PI6 RETURN ELSE IF(X .EQ. MONE) THEN DDILOG=MALF*PI6 RETURN END IF T=-X IF(T .LE. MTWO) THEN Y=MONE/(ONE+T) S=ONE A=-PI3+HALF*(LOG(-T)**2-LOG(ONE+ONE/T)**2) ELSE IF(T .LT. MONE) THEN Y=MONE-T S=MONE A=LOG(-T) A=-PI6+A*(A+LOG(ONE+ONE/T)) ELSE IF(T .LE. MALF) THEN Y=(MONE-T)/T S=ONE A=LOG(-T) A=-PI6+A*(MALF*A+LOG(ONE+T)) ELSE IF(T .LT. ZERO) THEN Y=-T/(ONE+T) S=MONE A=HALF*LOG(ONE+T)**2 ELSE IF(T .LE. ONE) THEN Y=T S=ONE A=ZERO ELSE Y=ONE/T S=MONE A=PI6+HALF*LOG(T)**2 END IF H=Y+Y-ONE ALFA=H+H B1=ZERO B2=ZERO DO 1 I = 18,0,-1 B0=C(I)+ALFA*B1-B2 B2=B1 1 B1=B0 DDILOG=-(S*(B0-H*B2)+A) RETURN END DOUBLE PRECISION FUNCTION TRILOG(X) * ====== ========= ======== ========= C Tsuyoshi Matsuura 1987 C C TRILOG: Li3 between 0 and 1 !!! C IMPLICIT NONE DOUBLE PRECISION X,S(0:10),L(0:20),U,Z,HELP,Z3 DOUBLE PRECISION DDILOG INTEGER I C C S: COEFFICIENTS OF S1,2 C L: COEFFICIENTS OF Li3 C DATA S/0.5000000000000000D+00,2.0833333333333333D-02, 1 -2.3148148148148148D-04,4.1335978835978837D-06, 2 -8.2671957671957673D-08,1.7397297489890083D-09, 3 -3.7744215276339238D-11,8.3640853316779243D-13, 4 -1.8831557201792127D-14,4.2930310281389223D-16, 5 -9.8857668116275541D-18/ DATA L/1.0000000000000000D+00,-0.3750000000000000D+00, 1 7.8703703703703705E-02,-8.6805555555555557E-03, 2 1.2962962962962963E-04, 8.1018518518518520E-05, 3 -3.4193571608537598E-06,-1.3286564625850340E-06, 4 8.6608717561098512E-08, 2.5260875955320400E-08, 5 -2.1446944683640648E-09,-5.1401106220129790E-10, 6 5.2495821146008296E-11, 1.0887754406636318E-11, 7 -1.2779396094493696E-12,-2.3698241773087452E-13, 8 3.1043578879654624E-14, 5.2617586299125060E-15, 9 -7.5384795499492655E-16,-1.1862322577752285E-16, 1 1.8316979965491384E-17/ DATA Z3/1.2020569031595943D+00/ C IF(X .LT. 0D0 .OR. X .GT. 1D0) THEN CCC WRITE(*,*)' **************************************' CCC WRITE(*,*)' Li3 called with X = ',X CCC WRITE(*,*)' This should lie between 0 and 1 !!!' CCC STOP' The program stops right here !!!' CCC CCC fortan writes do not work nicely with C++ code CCC Do whatever a standard math funcion would do with invalid arg. CCC TRILOG = DLOG(-1.D0) TRILOG = DLOG(1.D0) ENDIF IF (X.LT.0.5D0) THEN IF (X.GT.0.0D0) THEN U=-DLOG(1.0D0-X) HELP=U*L(20)+L(19) DO 10 I=18,0,-1 HELP=U*HELP+L(I) 10 CONTINUE TRILOG=U*HELP ELSE TRILOG=0.0D0 ENDIF ELSE IF (X.LT.1.0D0) THEN U=-DLOG(X) Z=U*U HELP=Z*S(10)+S(9) DO 20 I=8,0,-1 HELP=Z*HELP+S(I) 20 CONTINUE HELP=1.0D0/2.0D0*(Z*HELP-Z*U/6.0D0) C LI3=-HELP+DLOG(X)*DILOG(X)+0.5D0*DLOG(1D0-X)*DLOG(X)**2+Z3 TRILOG=-HELP+DLOG(X)*DDILOG(X)+0.5D0*DLOG(1D0-X)*DLOG(X)**2+Z3 ELSE TRILOG=Z3 ENDIF ENDIF END DOUBLE PRECISION FUNCTION S12(X) * ====== ========= ======== ====== C Tsuyoshi Matsuura 1987 C C S1,2 between 0 and 1 !!! C IMPLICIT NONE REAL*8 X,S(0:10),L(0:20),U,Z,HELP,Z3 REAL*8 DDILOG INTEGER I C C S: COEFFICIENTS OF S1,2 C L: COEFFICIENTS OF Li3 C DATA S/0.5000000000000000D+00,2.0833333333333333D-02, 1 -2.3148148148148148D-04,4.1335978835978837D-06, 2 -8.2671957671957673D-08,1.7397297489890083D-09, 3 -3.7744215276339238D-11,8.3640853316779243D-13, 4 -1.8831557201792127D-14,4.2930310281389223D-16, 5 -9.8857668116275541D-18/ DATA L/1.0000000000000000D+00,-0.3750000000000000D+00, 1 7.8703703703703705E-02,-8.6805555555555557E-03, 2 1.2962962962962963E-04, 8.1018518518518520E-05, 3 -3.4193571608537598E-06,-1.3286564625850340E-06, 4 8.6608717561098512E-08, 2.5260875955320400E-08, 5 -2.1446944683640648E-09,-5.1401106220129790E-10, 6 5.2495821146008296E-11, 1.0887754406636318E-11, 7 -1.2779396094493696E-12,-2.3698241773087452E-13, 8 3.1043578879654624E-14, 5.2617586299125060E-15, 9 -7.5384795499492655E-16,-1.1862322577752285E-16, 1 1.8316979965491384E-17/ DATA Z3/1.2020569031595943D+00/ C IF(X .LT. 0D0 .OR. X .GT. 1D0) THEN CCC WRITE(*,*)' **************************************' CCC WRITE(*,*)' S12 called with X = ',X CCC WRITE(*,*)' This should lie between 0 and 1 !!!' CCC STOP' The program stops right here !!!' CCC CCC fortan writes do not work nicely with C++ code CCC Do whatever a standard math funcion would do with invalid arg. CCC S12 = DLOG(-1.D0) S12 = DLOG(1.D0) ENDIF IF (X.LT.0.5D0) THEN IF (X.GT.0.0D0) THEN U=-DLOG(1.0D0-X) Z=U*U HELP=Z*S(10)+S(9) DO 10 I=8,0,-1 HELP=HELP*Z+S(I) 10 CONTINUE S12=1.0D0/2.0D0*(Z*HELP-Z*U/6.0D0) ELSE S12=0.0D0 ENDIF ELSE IF (X.LT.1.0D0) THEN U=-DLOG(X) HELP=U*L(20)+L(19) DO 20 I=18,0,-1 HELP=HELP*U+L(I) 20 CONTINUE HELP=U*HELP S12=-HELP+DLOG(1.0D0-X)*DDILOG(1.0D0-X)+ + 0.5D0*DLOG(X)*DLOG(1.0D0-X)**2+Z3 ELSE S12=Z3 ENDIF ENDIF END FUNCTION SPENCE(X) * IMPLICIT REAL*8(A-H,O-Z) PARAMETER (F1=1.64493406684822618D0) * IF(X)8,1,1 1 IF(X-.5D0)2,2,3 2 SPENCE=FSPENS(X) RETURN 3 IF(X-1.D0)4,4,5 4 SPENCE=F1-LOG(X)*LOG(1D0-X+1D-15)-FSPENS(1D0-X) RETURN 5 IF(X-2.D0)6,6,7 6 SPENCE=F1-.5D0*LOG(X)*LOG((X-1D0)**2/X)+FSPENS(1D0-1D0/X) RETURN 7 SPENCE=2D0*F1-.5D0*LOG(X)**2-FSPENS(1D0/X) RETURN 8 IF(X+1.D0)10,9,9 9 SPENCE=-.5D0*LOG(1D0-X)**2-FSPENS(X/(X-1D0)) RETURN 10 SPENCE=-.5D0*LOG(1D0-X)*LOG(X**2/(1D0-X))-F1+FSPENS(1D0/(1D0-X)) END FUNCTION FSPENS(X) * IMPLICIT REAL*8(A-H,O-Z) * A=1D0 F=0D0 AN=0D0 TCH=1D-16 1 AN=AN+1D0 A=A*X B=A/AN**2 F=F+B IF(B-TCH)2,2,1 2 FSPENS=F END C-----------------------------------------------------------------------