Commit 3148ce21 authored by ulrich_y's avatar ulrich_y
Browse files

First public release of McMule

parent 9d6d323f
Pipeline #1223 passed with stages
in 33 minutes and 12 seconds
[submodule "tools/pymule"]
path = tools/pymule
url = ../pymule.git
url = https://gitlab.psi.ch/mcmule/pymule.git
branch = master
......@@ -3,7 +3,7 @@ FROM yulrich/mcmule-pre:1.0.0
LABEL maintainer="yannick.ulrich@psi.ch"
LABEL version="1.0"
LABEL description="The full McMule suite"
RUN pip3 install git+https://"gitlab+deploy-token-2":TToX8HiCR5mV6hLH_UXD@gitlab.psi.ch/mcmule/pymule.git
RUN pip3 install git+https://gitlab.psi.ch/mcmule/pymule.git
COPY . /monte-carlo
RUN cd /monte-carlo && ./configure
RUN cd /monte-carlo && make
This diff is collapsed.
# The McMule internal repository
McMule is a framework for fully differential higher-order QED
calculations of scattering and decay processes involving leptons. It
keeps finite lepton masses, which regularises collinear singularities.
Soft singularities are treated with dimensional regularisation and
using $`{\rm FKS}^\ell`$ subtraction
Please find the manual for McMule
[here](https://gitlab.psi.ch/mcmule/manual) or download it here [as a
pdf](https://gitlab.psi.ch/mcmule/manual/-/jobs/artifacts/master/raw/manual.pdf?job=compile)
If you find McMule useful, please consider citing
[[2007.xxxxx]](https://arxiv.org/abs/2007.xxxxx)
> QED at NNLO with McMule
>
> P. Banerjee, T. Engel, A. Signer, Y. Ulrich
## Getting started
To obtain a copy of McMule we recommend the following approach
```bash
$ git clone --recursive https://gitlab.psi.ch/mcmule/mcmule
````
To build McMule, a Fortran compiler such as `gfortran` and a python
installation is needed. The main executable can be compiled by running
```bash
$ ./configure
$ make mcmule
```
Alternatively, we provide a Docker container for easy deployment and
legacy results. In multi-user environments, udocker can be used
instead. In either case, a pre-compiled copy of the code can be
obtained by calling
```bash
$ docker pull yulrich/mcmule # requires Docker to be installed
$ udocker pull yulrich/mcmule # requires uDocker to be installed
```
Please refer to the manual on how to use McMule.
group=mudec
AUXFILES= mudec_pm2enngg.f95 mudec_pm2ennggav.f95 \
AUXFILES= mudec_pm2enngg.f95 mudec_pm2ennggav.f95 mudec_h2lff.f95 \
mudec_pm2ennglav.flex.opt.f95 mudec_TLEPS0.f95 mudec_TLEPS1.f95
......
......@@ -3,7 +3,7 @@
MODULE mudec_TLEPS0
!!!!!!!!!!!!!!!!!!!!!!
contains
function TLEPS0(EL,GF,MOU2,MIN2,p1p2,p1p4,p2p4,asym,np2,np3, &
function TLEPS0(EL,GF,MOU2,MIN2,p1p2,p1p4,p2p4,asym,np2,np3, &
np4,mu2)
!
use global_def, only: pi, prec
......@@ -30,74 +30,74 @@
integer, parameter :: dd2 = 7
integer, parameter :: dd3 = 10
integer, parameter :: dd00 = 13
external A0, B0i
external B0
external DB0
external C0, C0i
external D0, D0i
!
real(kind=prec) TLEPS0,EL,GF,MOU2,MIN2, &
p1p2,p1p4,p2p4,asym, &
np2,np3,np4,mu2, &
real(kind=prec) TLEPS0,EL,GF,MOU2,MIN2, &
p1p2,p1p4,p2p4,asym, &
np2,np3,np4,mu2, &
aa,bb
!
aa=0._prec
bb=-1._prec
!
!
!
! call setlambda(-1._prec)
! call setdelta(0._prec)
! call setmudim(mu2)
! call setminmass(0._prec)
!
!
!
TLEPS0=real((-32*EL**2*GF**2*(MOU2**2*p1p4* &
(MIN2*p1p4 - 3*sqrt(MIN2)*np4*p2p4 + &
3*p1p4*(-p1p2 + p2p4)) + &
MOU2*(3*MIN2**1.5*np4*(2*p1p4 - p2p4)*p2p4 + &
MIN2**2*(-2*p1p4**2 + p2p4**2) + &
p1p4*(-4*p1p4**3 + 8*p1p4**2*p2p4 + &
2*p1p2*(7*p1p4 - 3*p2p4)*p2p4 - &
11*p1p4*p2p4**2 + 3*p2p4**3 + &
p1p2**2*(-10*p1p4 + 6*p2p4)) + &
sqrt(MIN2)*np4* &
(4*p1p4**3 - 8*p1p4**2*p2p4 + &
3*(p1p2 - p2p4)*p2p4**2 + &
p1p4*p2p4*(-7*p1p2 + 11*p2p4)) + &
MIN2*(-8*p1p4**2*p2p4 + 3*p2p4**3 + &
p1p2*(10*p1p4**2 - 2*p1p4*p2p4 - 3*p2p4**2))) &
+ (p1p4 - p2p4)* &
(-3*MIN2**2.5*np4*p2p4 + MIN2**3*(p1p4 + p2p4) + &
MIN2**1.5*np4* &
(-4*p1p4**2 + 4*p1p4*p2p4 + &
7*(p1p2 - p2p4)*p2p4) - &
4*p1p4*(p1p2 - p2p4)* &
(2*p1p2**2 + 2*p1p4**2 - 2*p1p2*p2p4 - &
2*p1p4*p2p4 + p2p4**2) + &
MIN2*(4*p1p4**3 - 4*p1p4**2*p2p4 + &
5*p1p4*p2p4**2 + 4*p2p4**3 + &
2*p1p2**2*(7*p1p4 + 2*p2p4) - &
2*p1p2*p2p4*(9*p1p4 + 4*p2p4)) + &
MIN2**2*(5*p2p4*(p1p4 + p2p4) - &
p1p2*(7*p1p4 + 5*p2p4)) - &
4*sqrt(MIN2)*np4* &
(p1p2**2*p2p4 + &
p2p4*(2*p1p4**2 - 2*p1p4*p2p4 + p2p4**2) - &
2*p1p2*(p1p4**2 - p1p4*p2p4 + p2p4**2))) + &
sqrt(MIN2)*np2* &
(3*MOU2**2*p1p4**2 + &
(p1p4 - p2p4)*(3*MIN2 - 4*p1p2 + 4*p2p4)* &
(p1p4*(-2*p1p2 + p2p4) + MIN2*(p1p4 + p2p4)) + &
MOU2*(MIN2*(-6*p1p4**2 + 3*p2p4**2) + &
p1p4*(10*p1p2*p1p4 - 6*p1p2*p2p4 - &
7*p1p4*p2p4 + 3*p2p4**2)))))/ &
TLEPS0=real((-32*EL**2*GF**2*(MOU2**2*p1p4* &
(MIN2*p1p4 - 3*sqrt(MIN2)*np4*p2p4 + &
3*p1p4*(-p1p2 + p2p4)) + &
MOU2*(3*MIN2**1.5*np4*(2*p1p4 - p2p4)*p2p4 + &
MIN2**2*(-2*p1p4**2 + p2p4**2) + &
p1p4*(-4*p1p4**3 + 8*p1p4**2*p2p4 + &
2*p1p2*(7*p1p4 - 3*p2p4)*p2p4 - &
11*p1p4*p2p4**2 + 3*p2p4**3 + &
p1p2**2*(-10*p1p4 + 6*p2p4)) + &
sqrt(MIN2)*np4* &
(4*p1p4**3 - 8*p1p4**2*p2p4 + &
3*(p1p2 - p2p4)*p2p4**2 + &
p1p4*p2p4*(-7*p1p2 + 11*p2p4)) + &
MIN2*(-8*p1p4**2*p2p4 + 3*p2p4**3 + &
p1p2*(10*p1p4**2 - 2*p1p4*p2p4 - 3*p2p4**2))) &
+ (p1p4 - p2p4)* &
(-3*MIN2**2.5*np4*p2p4 + MIN2**3*(p1p4 + p2p4) + &
MIN2**1.5*np4* &
(-4*p1p4**2 + 4*p1p4*p2p4 + &
7*(p1p2 - p2p4)*p2p4) - &
4*p1p4*(p1p2 - p2p4)* &
(2*p1p2**2 + 2*p1p4**2 - 2*p1p2*p2p4 - &
2*p1p4*p2p4 + p2p4**2) + &
MIN2*(4*p1p4**3 - 4*p1p4**2*p2p4 + &
5*p1p4*p2p4**2 + 4*p2p4**3 + &
2*p1p2**2*(7*p1p4 + 2*p2p4) - &
2*p1p2*p2p4*(9*p1p4 + 4*p2p4)) + &
MIN2**2*(5*p2p4*(p1p4 + p2p4) - &
p1p2*(7*p1p4 + 5*p2p4)) - &
4*sqrt(MIN2)*np4* &
(p1p2**2*p2p4 + &
p2p4*(2*p1p4**2 - 2*p1p4*p2p4 + p2p4**2) - &
2*p1p2*(p1p4**2 - p1p4*p2p4 + p2p4**2))) + &
sqrt(MIN2)*np2* &
(3*MOU2**2*p1p4**2 + &
(p1p4 - p2p4)*(3*MIN2 - 4*p1p2 + 4*p2p4)* &
(p1p4*(-2*p1p2 + p2p4) + MIN2*(p1p4 + p2p4)) + &
MOU2*(MIN2*(-6*p1p4**2 + 3*p2p4**2) + &
p1p4*(10*p1p2*p1p4 - 6*p1p2*p2p4 - &
7*p1p4*p2p4 + 3*p2p4**2)))))/ &
(3.*p1p4**2*(p1p4 - p2p4)**2))
!
! T1BOXAC=T1BOXAC+real(
......
......@@ -3,7 +3,7 @@
MODULE mudec_TLEPS1
!!!!!!!!!!!!!!!!!!!!!!
contains
function TLEPS1(EL,GF,MOU2,MIN2,p1p2,p1p4,p2p4,asym,np2,np3, &
function TLEPS1(EL,GF,MOU2,MIN2,p1p2,p1p4,p2p4,asym,np2,np3, &
np4,mu2)
!
use global_def, only: pi, prec
......@@ -30,78 +30,78 @@
integer, parameter :: dd2 = 7
integer, parameter :: dd3 = 10
integer, parameter :: dd00 = 13
external A0, B0i
external B0
external DB0
external C0, C0i
external D0, D0i
!
real(kind=prec) TLEPS1,EL,GF,MOU2,MIN2, &
p1p2,p1p4,p2p4,asym, &
np2,np3,np4,mu2, &
real(kind=prec) TLEPS1,EL,GF,MOU2,MIN2, &
p1p2,p1p4,p2p4,asym, &
np2,np3,np4,mu2, &
aa,bb
!
aa=0._prec
bb=-1._prec
!
!
!
! call setlambda(-1._prec)
! call setdelta(0._prec)
! call setmudim(mu2)
! call setminmass(0._prec)
!
!
!
TLEPS1=real((16*EL**2*GF**2*(2*MOU2**2*p1p4* &
(7*MIN2*p1p4 - 9*sqrt(MIN2)*np4*p2p4 + &
9*p1p4*(-p1p2 + p2p4)) + &
2*sqrt(MIN2)*np2* &
(9*MOU2**2*p1p4**2 + &
(p1p4 - p2p4)*(9*MIN2 + 16*(-p1p2 + p2p4))* &
(p1p4*(-2*p1p2 + p2p4) + MIN2*(p1p4 + p2p4)) + &
MOU2*(9*MIN2*(-2*p1p4**2 + p2p4**2) + &
p1p4*(34*p1p2*p1p4 - 18*p1p2*p2p4 - &
25*p1p4*p2p4 + 9*p2p4**2))) + &
MOU2*(18*MIN2**1.5*np4*(2*p1p4 - p2p4)*p2p4 + &
14*MIN2**2*(-2*p1p4**2 + p2p4**2) - &
4*p1p4*(8*p1p4**3 + &
p1p2**2*(17*p1p4 - 9*p2p4) - &
16*p1p4**2*p2p4 + 25*p1p4*p2p4**2 - &
9*p2p4**3 + p1p2*p2p4*(-25*p1p4 + 9*p2p4)) - &
2*MIN2*(32*p1p4**2*p2p4 - 9*p2p4**3 + &
p1p2*(-46*p1p4**2 + 14*p1p4*p2p4 + 9*p2p4**2) &
) + sqrt(MIN2)* &
(3*CMPLX(aa,bb)*asym*p2p4*(-p1p4 + p2p4) + &
2*np4*(16*p1p4**3 - 32*p1p4**2*p2p4 - &
25*p1p4*(p1p2 - 2*p2p4)*p2p4 + &
9*(p1p2 - 2*p2p4)*p2p4**2))) + &
(p1p4 - p2p4)* &
(-18*MIN2**2.5*np4*p2p4 + &
14*MIN2**3*(p1p4 + p2p4) - &
8*p1p4*(p1p2 - p2p4)* &
(8*p1p2**2 + 8*p1p4**2 - 8*p1p2*p2p4 - &
8*p1p4*p2p4 + 7*p2p4**2) + &
4*MIN2*(8*p1p4**3 - 8*p1p4**2*p2p4 + &
13*p1p4*p2p4**2 + 8*p2p4**3 + &
p1p2**2*(31*p1p4 + 8*p2p4) - &
p1p2*p2p4*(39*p1p4 + 16*p2p4)) + &
MIN2**2*(46*p2p4*(p1p4 + p2p4) - &
2*p1p2*(37*p1p4 + 23*p2p4)) + &
MIN2**1.5*(3*CMPLX(aa,bb)*asym*p2p4 + &
np4*(-32*p1p4**2 + 50*p1p2*p2p4 + &
32*p1p4*p2p4 - 68*p2p4**2)) + &
sqrt(MIN2)* &
(3*CMPLX(aa,bb)*asym*p1p4*(-p1p2 + p1p4 + p2p4) - &
8*np4*(-8*p1p2*p1p4**2 + 4*p1p2**2*p2p4 + &
8*p1p2*p1p4*p2p4 + 8*p1p4**2*p2p4 - &
11*p1p2*p2p4**2 - 8*p1p4*p2p4**2 + &
7*p2p4**3)))))/ &
TLEPS1=real((16*EL**2*GF**2*(2*MOU2**2*p1p4* &
(7*MIN2*p1p4 - 9*sqrt(MIN2)*np4*p2p4 + &
9*p1p4*(-p1p2 + p2p4)) + &
2*sqrt(MIN2)*np2* &
(9*MOU2**2*p1p4**2 + &
(p1p4 - p2p4)*(9*MIN2 + 16*(-p1p2 + p2p4))* &
(p1p4*(-2*p1p2 + p2p4) + MIN2*(p1p4 + p2p4)) + &
MOU2*(9*MIN2*(-2*p1p4**2 + p2p4**2) + &
p1p4*(34*p1p2*p1p4 - 18*p1p2*p2p4 - &
25*p1p4*p2p4 + 9*p2p4**2))) + &
MOU2*(18*MIN2**1.5*np4*(2*p1p4 - p2p4)*p2p4 + &
14*MIN2**2*(-2*p1p4**2 + p2p4**2) - &
4*p1p4*(8*p1p4**3 + &
p1p2**2*(17*p1p4 - 9*p2p4) - &
16*p1p4**2*p2p4 + 25*p1p4*p2p4**2 - &
9*p2p4**3 + p1p2*p2p4*(-25*p1p4 + 9*p2p4)) - &
2*MIN2*(32*p1p4**2*p2p4 - 9*p2p4**3 + &
p1p2*(-46*p1p4**2 + 14*p1p4*p2p4 + 9*p2p4**2) &
) + sqrt(MIN2)* &
(3*CMPLX(aa,bb)*asym*p2p4*(-p1p4 + p2p4) + &
2*np4*(16*p1p4**3 - 32*p1p4**2*p2p4 - &
25*p1p4*(p1p2 - 2*p2p4)*p2p4 + &
9*(p1p2 - 2*p2p4)*p2p4**2))) + &
(p1p4 - p2p4)* &
(-18*MIN2**2.5*np4*p2p4 + &
14*MIN2**3*(p1p4 + p2p4) - &
8*p1p4*(p1p2 - p2p4)* &
(8*p1p2**2 + 8*p1p4**2 - 8*p1p2*p2p4 - &
8*p1p4*p2p4 + 7*p2p4**2) + &
4*MIN2*(8*p1p4**3 - 8*p1p4**2*p2p4 + &
13*p1p4*p2p4**2 + 8*p2p4**3 + &
p1p2**2*(31*p1p4 + 8*p2p4) - &
p1p2*p2p4*(39*p1p4 + 16*p2p4)) + &
MIN2**2*(46*p2p4*(p1p4 + p2p4) - &
2*p1p2*(37*p1p4 + 23*p2p4)) + &
MIN2**1.5*(3*CMPLX(aa,bb)*asym*p2p4 + &
np4*(-32*p1p4**2 + 50*p1p2*p2p4 + &
32*p1p4*p2p4 - 68*p2p4**2)) + &
sqrt(MIN2)* &
(3*CMPLX(aa,bb)*asym*p1p4*(-p1p2 + p1p4 + p2p4) - &
8*np4*(-8*p1p2*p1p4**2 + 4*p1p2**2*p2p4 + &
8*p1p2*p1p4*p2p4 + 8*p1p4**2*p2p4 - &
11*p1p2*p2p4**2 - 8*p1p4*p2p4**2 + &
7*p2p4**3)))))/ &
(9.*p1p4**2*(p1p4 - p2p4)**2))
!
! T1BOXAC=T1BOXAC+real(
......
!!!!!!!!!!!!!!!!!!!!!!
MODULE mudec_h2lff
!!!!!!!!!!!!!!!!!!!!!!
use functions
implicit none
real(kind=prec), parameter :: zeta(2:5) = (/1.6449340668482264365, &
1.2020569031595942854, &
1.0823232337111381915, &
1.0369277551433699263 /)
real(kind=8) Hr1(-1:+1)
real(kind=8) Hr2(-1:+1,-1:+1)
real(kind=8) Hr3(-1:+1,-1:+1,-1:+1)
real(kind=8) Hr4(-1:+1,-1:+1,-1:+1,-1:+1)
real(kind=8) Hex(2)
real(kind=prec) :: logz(1:4)
real(kind=prec) :: xpow(-3:5)
real(kind=prec) :: xm1(2:4)
contains
SUBROUTINE INIT_HPLS(x, z)
real(kind=prec), intent(in) :: x, z
integer, parameter :: n1 = -1
integer, parameter :: n2 = +1
integer, parameter :: nw = 4
complex(kind=16) Hc1(-1:+1)
complex(kind=16) Hc2(-1:+1,-1:+1)
complex(kind=16) Hc3(-1:+1,-1:+1,-1:+1)
complex(kind=16) Hc4(-1:+1,-1:+1,-1:+1,-1:+1)
real(kind=8) Hi1(-1:+1)
real(kind=8) Hi2(-1:+1,-1:+1)
real(kind=8) Hi3(-1:+1,-1:+1,-1:+1)
real(kind=8) Hi4(-1:+1,-1:+1,-1:+1,-1:+1)
call hplog(x, nw, Hc1, Hc2, Hc3, Hc4, &
Hr1, Hr2, Hr3, Hr4, &
Hi1, Hi2, Hi3, Hi4, n1, n2)
xpow = x**(/-3,-2,-1,0,1,2,3,4,5/)
xm1 = 1._prec/(x-1)**(/2,3,4/)
logz = log(z)**(/1,2,3,4/)
Hex(1) = log(1-0.5*x)
Hex(2) = 0.5*log(1-x)**2 * log(2-x) + log(1-x) * Li2(x-1) - Li3(x-1) - 3. * zeta(3) / 4.
END SUBROUTINE
FUNCTION FF1A1ES()
! Calculate the alpha^1 e^-1 coeff. of F1
implicit none
real(kind=prec) :: ff1a1es
ff1a1es = logz(1)*(-1)
ff1a1es = ff1a1es -1
ff1a1es = ff1a1es + Hr1(1)*(-1)
END FUNCTION FF1A1ES
FUNCTION FF1A1EF()
! Calculate the alpha^1 e^0 coeff. of F1
implicit none
real(kind=prec) :: ff1a1ef
ff1a1ef = logz(1)*(-0.5)
ff1a1ef = ff1a1ef + logz(2)*(1)
ff1a1ef = ff1a1ef -2
ff1a1ef = ff1a1ef + Hr2(1,1)*(-2)
ff1a1ef = ff1a1ef + Hr1(1)*(-1.5)
ff1a1ef = ff1a1ef + Hr2(0,1)*(-1)
END FUNCTION FF1A1EF
FUNCTION FF1A1EL()
! Calculate the alpha^1 e^1 coeff. of F1
implicit none
real(kind=prec) :: ff1a1el
ff1a1el = Hr3(0,0,1)*(-1)
ff1a1el = ff1a1el + logz(3)*(-(2./3.))
ff1a1el = ff1a1el + logz(2)*(0.5)
ff1a1el = ff1a1el + Hr1(1)*((-7 - 2*zeta(2))/2.)
ff1a1el = ff1a1el + logz(1)*((-3 - 2*zeta(2))/2.)
ff1a1el = ff1a1el -4 - zeta(2)
ff1a1el = ff1a1el + Hr3(1,1,1)*(-4)
ff1a1el = ff1a1el + Hr2(1,1)*(-3)
ff1a1el = ff1a1el + Hr3(0,1,1)*(-2)
ff1a1el = ff1a1el + Hr3(1,0,1)*(-2)
ff1a1el = ff1a1el + Hr2(0,1)*(-1.5)
END FUNCTION FF1A1EL
FUNCTION FF1A2ED()
! Calculate the alpha^2 e^-2 coeff. of F1
implicit none
real(kind=prec) :: ff1a2ed
ff1a2ed = logz(1)*(1)
ff1a2ed = ff1a2ed + Hr1(1)*logz(1)*(1)
ff1a2ed = ff1a2ed + 0.5
ff1a2ed = ff1a2ed + logz(2)*(0.5)
ff1a2ed = ff1a2ed + Hr1(1)*(1)
ff1a2ed = ff1a2ed + Hr2(1,1)*(1)
END FUNCTION FF1A2ED
FUNCTION FF1A2ES()
! Calculate the alpha^2 e^-1 coeff. of F1
implicit none
real(kind=prec) :: ff1a2es
ff1a2es = Hr2(1,1)*(5)
ff1a2es = ff1a2es + Hr3(1,1,1)*(6)
ff1a2es = ff1a2es + 2
ff1a2es = ff1a2es + Hr1(1)*logz(2)*(-1)
ff1a2es = ff1a2es + logz(3)*(-1)
ff1a2es = ff1a2es + logz(2)*(-0.5)
ff1a2es = ff1a2es + Hr2(0,1)*(1)
ff1a2es = ff1a2es + Hr3(1,0,1)*(1)
ff1a2es = ff1a2es + Hr2(0,1)*logz(1)*(1)
ff1a2es = ff1a2es + Hr3(0,1,1)*(2)
ff1a2es = ff1a2es + Hr1(1)*logz(1)*(2)
ff1a2es = ff1a2es + Hr2(1,1)*logz(1)*(2)
ff1a2es = ff1a2es + logz(1)*(2.5)
ff1a2es = ff1a2es + Hr1(1)*(3.5)
END FUNCTION FF1A2ES
FUNCTION FF1A2EF()
! Calculate the alpha^2 e^0 coeff. of F1
implicit none
real(kind=prec) :: ff1a2ef
ff1a2ef = Hr3(0,1,1)*logz(1)*(2)
ff1a2ef = ff1a2ef + Hr3(1,0,1)*logz(1)*(2)
ff1a2ef = ff1a2ef + Hr2(1,1)*logz(1)*(4)
ff1a2ef = ff1a2ef + Hr3(1,1,1)*logz(1)*(4)
ff1a2ef = ff1a2ef + Hr4(1,1,0,1)*(6)
ff1a2ef = ff1a2ef + Hr4(1,0,1,1)*(10)
ff1a2ef = ff1a2ef + Hr4(0,1,1,1)*(12)
ff1a2ef = ff1a2ef + Hr3(1,1,1)*(22)
ff1a2ef = ff1a2ef + Hr4(1,1,1,1)*(28)
ff1a2ef = ff1a2ef + Hr3(-1,0,1)*(2*xm1(3)*(1 + xpow(1))*(9 - 8*xpow(1) + 5*xpow(2)))
ff1a2ef = ff1a2ef + Hr4(0,0,1,1)*(2*xm1(4)*xpow(1)*(-28 - 9*xpow(1) - 6*xpow(2) + xpow(3&
&)))
ff1a2ef = ff1a2ef + Hex(2)*(-(xm1(4)*(-12 + 6*xpow(1) - 5*xpow(2) + 2*xpow(3)))/2.)
ff1a2ef = ff1a2ef + Hr3(1,0,1)*(xm1(3)*(-12 - 16*xpow(1) - 19*xpow(2) + 5*xpow(3)))
ff1a2ef = ff1a2ef + Hr4(0,-1,0,1)*(-4*xm1(4)*(1 + 8*xpow(2) - 4*xpow(3) + xpow(4)))
ff1a2ef = ff1a2ef + Hr4(0,0,0,1)*(-(xm1(4)*(3 + 19*xpow(2) - 6*xpow(3) + 2*xpow(4))))
ff1a2ef = ff1a2ef + Hr4(0,1,0,1)*(xm1(4)*(3 + 16*xpow(1) + 27*xpow(2) - 6*xpow(3) + 2*xp&
&ow(4)))
ff1a2ef = ff1a2ef + Hr3(0,0,1)*(-(xm1(4)*(7 + 14*xpow(1) + 8*xpow(2) - 28*xpow(3) + 8*xp&
&ow(4)))/2.)
ff1a2ef = ff1a2ef + Hr3(0,1,1)*((xm1(4)*(-18 - 122*xpow(1) + 185*xpow(2) - 62*xpow(3) + &
&26*xpow(4)))/2.)
ff1a2ef = ff1a2ef + Hr1(-1)*(xm1(3)*(1 + xpow(1))*(9 - 8*xpow(1) + 5*xpow(2))*zeta(2))
ff1a2ef = ff1a2ef + Hex(1)*((-3*xm1(4)*(-12 + 6*xpow(1) - 5*xpow(2) + 2*xpow(3))*zeta(2)&
&)/2.)
ff1a2ef = ff1a2ef + Hr2(0,-1)*(-2*xm1(4)*(1 + 8*xpow(2) - 4*xpow(3) + xpow(4))*zeta(2))
ff1a2ef = ff1a2ef + logz(2)*((-7 + 8*zeta(2))/8.)
ff1a2ef = ff1a2ef + Hr1(1)*logz(1)*((23 + 8*zeta(2))/4.)
ff1a2ef = ff1a2ef + Hr2(1,1)*((xm1(2)*(-23 - 66*xpow(1) + 14*xpow(2) + 4*zeta(2) - 8*xpo&
&w(1)*zeta(2) + 4*xpow(2)*zeta(2)))/2.)
ff1a2ef = ff1a2ef + Hr2(0,1)*(-(xm1(4)*xpow(-1)*(2 - 11*xpow(1) - 60*xpow(2) + 18*xpow(3&
&) + 28*xpow(4) + 23*xpow(5) - 12*xpow(1)*zeta(2) - 224*xpow(2)*zeta(2) - 148*xpow(3)*zeta(&
&2) - 24*xpow(4)*zeta(2)))/4.)
ff1a2ef = ff1a2ef + logz(1)*((49 + 40*zeta(2) - 48*zeta(3))/8.)
ff1a2ef = ff1a2ef + Hr1(1)*((xm1(3)*xpow(-1)*(-2 - 71*xpow(1) + 229*xpow(2) - 237*xpow(3&
&) + 81*xpow(4) - 208*xpow(1)*zeta(2) - 400*xpow(2)*zeta(2) - 248*xpow(3)*zeta(2) + 40*xpow&
&(4)*zeta(2) + 16*xpow(1)*zeta(3) - 48*xpow(2)*zeta(3) + 48*xpow(3)*zeta(3) - 16*xpow(4)*ze&
&ta(3)))/8.)
ff1a2ef = ff1a2ef -(xm1(4)*(-92 + 368*xpow(1) - 552*xpow(2) + 368*xpow(3) - 92*xpow(4) +&
& 88*zeta(2) - 336*log2*zeta(2) - 676*xpow(1)*zeta(2) + 552*log2*xpow(1)*zeta(2) + 344*xpow&
&(2)*zeta(2) - 348*log2*xpow(2)*zeta(2) + 280*xpow(3)*zeta(2) - 72*log2*xpow(3)*zeta(2) - 3&
&6*xpow(4)*zeta(2) + 96*log2*xpow(4)*zeta(2) + 80*zeta(3) - 238*xpow(1)*zeta(3) + 205*xpow(&
&2)*zeta(3) - 66*xpow(3)*zeta(3) - 8*xpow(4)*zeta(3) + 210*zeta(4) + 376*xpow(1)*zeta(4) + &
&1598*xpow(2)*zeta(4) - 516*xpow(3)*zeta(4) + 156*xpow(4)*zeta(4)))/8.
ff1a2ef = ff1a2ef + Hr4(1,0,0,1)*(-3)
ff1a2ef = ff1a2ef + Hr1(1)*logz(2)*(-2)
ff1a2ef = ff1a2ef + Hr2(1,1)*logz(2)*(-2)
ff1a2ef = ff1a2ef + Hr2(0,1)*logz(2)*(-1)
ff1a2ef = ff1a2ef + logz(3)*(-(1./3.))
ff1a2ef = ff1a2ef + Hr1(1)*logz(3)*((2./3.))
ff1a2ef = ff1a2ef + Hr3(0,0,1)*logz(1)*(1)
ff1a2ef = ff1a2ef + logz(4)*((7./6.))
ff1a2ef = ff1a2ef + Hr2(0,1)*logz(1)*(2)
END FUNCTION FF1A2EF
FUNCTION FF2A1EF()
! Calculate the alpha^1 e^0 coeff. of F2
implicit none
real(kind=prec) :: ff2a1ef
ff2a1ef = Hr1(1)*(-xpow(-1))
END FUNCTION FF2A1EF
FUNCTION FF2A1EL()
! Calculate the alpha^1 e^1 coeff. of F2
implicit none
real(kind=prec) :: ff2a1el
ff2a1el = Hr2(0,