(* ::Package:: *)

(*
 anQCDv2: a Mathematica package for analytic running of the
 strong coupling in FAPT, 2danQCD and MPT models.

 Copyright: C. Ayala and G.Cvetic (Aug. 2014), Univ.Tecnica Federico Santa Maria
 Permission is granted for anyone to use this program for any purpose 
 on any computer. The program can be freely distributed. 
 The program has been checked to work in Mathematica 9.0.1, 10.0.1 and 10.0.2 *)

(* ************************************************************ *)

Print["anQCDv2: a Mathematica package for analytic running of the"];
Print["strong coupling in FAPT, 2danQCD and MPT models."];
Print["by C. Ayala and G.Cvetic (Sept. 2014)"];
Print["For details, see Comput.Phys.Commun.190(2015)182-199."]
Print["The software has been checked to work in:"]
Print["  Mathematica 9, 10, 11"];
Print["Adjustment to make it work in Mathematica 11 made in Aug. 2017"]
Print["The software can be freely distributed."];
(* ************************************************************ *)
(* ************************************************************ *)

BeginPackage["anQCD`"];



(* ::Section:: *)
(*Definitions*)

(*----------------------------------------------------------*)
(* pQCD and pQCD global: *)
a1l::usage = "a1l[Nf, Q2, L2, fi] gives the 1-loop QCD coupling with Nf active quarks at squared momentum Q2 and with squared Lambda QCD scale L2={Bar Lambda}^2_{Nf} (both Q2 and L2 in GeV^2).";
a2l::usage = "a2l[Nf, Q2, L2, fi] gives the 2-loop QCD coupling with Nf active quarks at squared momentum Q2 and with squared Lambda QCD scale L2={Bar Lambda}^2_{Nf} (both Q2 and L2 in GeV^2).";
a3l::usage = "a3l[Nf, Q2, L2, fi] gives the 3-loop QCD coupling in MSbar scheme with Nf active quarks at squared momentum Q2 and with square Lambda QCD scale L2={Bar Lambda}^2_{Nf} (both Q2 and L2 in GeV^2).";
a4l::usage = "a4l[Nf, Q2, L2, fi] gives the 4-loop QCD coupling in MSbar scheme with Nf active quarks at squared momentum Q2 and with square Lambda QCD scale L2={Bar Lambda}^2_{Nf} (both Q2 and L2 in GeV^2).";

a1lglob::usage = "a1lglob[Q2, L2Nf3, fi] gives the 1-loop global QCD coupling at Q2=Q^2 and L21lNf3 squared Lambda scale at Nf=3 active flavors to 1-loop, where L2Nf3={Bar Lambda}^2_{Nf=3} (all scales in GeV^2)";
a2lglob::usage = "a2lglob[Q2, L2Nf3, fi] gives the 2-loop global QCD coupling at Q2=Q^2 and L22lNf3 squared Lambda scale at Nf=3 active flavors to 2-loop, where L2Nf3={Bar Lambda}^2_{Nf=3} (all scales in GeV^2).";
a3lglob::usage = "a3lglob[Q2, L2Nf3, fi] gives the 3-loop global QCD coupling at Q2=Q^2 and L23lNf3 squared Lambda scale at Nf=3 active flavors to 3-loop, where L2Nf3={Bar Lambda}^2_{Nf=3} (all scales in GeV^2).";
a4lglob::usage = "a4lglob[Q2, L2Nf3, fi] gives the 4-loop global QCD coupling at Q2=Q^2 and L24lNf3 squared Lambda scale at Nf=3 active flavors to 4-loop, where L2Nf3={Bar Lambda}^2_{Nf=3} (all scales in GeV^2).";

tr1l::usage = "tr1l[Nf, nu, k, sig, L2] gives the 1-loop dimensionless spectral density with the inclusion of logarithm of the coupling in the form Im[a1l^nu Log[a1l]^k] at Q^2=-sig - i*epsilon, at Nf active flavors; nu real; k=0,1,2,...; L2={Bar Lambda}^2_{Nf} (all scales in GeV^2).";
tr2l::usage = "tr2l[Nf, nu, k, sig, L2] gives the 2-loop dimensionless spectral density with the inclusion of logarithm in the form Im[a2l^nu Log[a2l]^k] at Q^2=-sig - i*epsilon, at Nf active flavors; nu real; k=0,1,2,...; L2={Bar Lambda}^2_{Nf} (all scales in GeV^2).";
tr3l::usage = "tr3l[Nf, nu, k, sig, L2] gives the 3-loop dimensionless spectral density with the inclusion of logarithm in the form Im[a3l^nu Log[a3l]^k] at Q^2=-sig - i*epsilon, at Nf active flavors; nu real; k=0,1,2,...; L2={Bar Lambda}^2_{Nf} (all scales in GeV^2).";
tr4l::usage = "tr4l[Nf, nu, k, sig, L2] gives the 4-loop dimensionless spectral density with the inclusion of logarithm in the form Im[a4l^nu Log[a4l]^k] at Q^2=-sig - i*epsilon, at Nf active flavors; nu real; k=0,1,2,...; L2={Bar Lambda}^2_{Nf} (all scales in GeV^2).";

tr1lglob::usage = "tr1lglob[nu, k, sig, L2Nf3] gives the 1-loop dimensionless global spectral density with the inclusion of logarithm in the form Im[a1lglob^nu Log[a1lglob]^k] at Q^2=-sig - i*epsilon; nu real; k=0,1,2,...; L2Nf3={Bar Lambda}^2_{Nf=3} (all scales in GeV^2).";
tr2lglob::usage = "tr2lglob[nu, k, sig, L2Nf3] gives the 2-loop dimensionless global spectral density with the inclusion of logarithm in the form Im[a2lglob^nu Log[a2lglob]^k] at Q^2=-sig - i*epsilon; nu real; k=0,1,2,...; L2Nf3={Bar Lambda}^2_{Nf=3} (all scales in GeV^2).";
tr3lglob::usage = "tr3lglob[nu, k, sig, L2Nf3] gives the 3-loop dimensionless global spectral density with the inclusion of logarithm in the form Im[a3lglob^nu Log[a3lglob]^k] at Q^2=-sig - i*epsilon; nu real; k=0,1,2,...; L2Nf3={Bar Lambda}^2_{Nf=3} (all scales in GeV^2).";
tr4lglob::usage = "tr4lglob[nu, k, sig, L2Nf3] gives the 4-loop dimensionless global spectral density with the inclusion of logarithm in the form Im[a4lglob^nu Log[a4lglob]^k] at Q^2=-sig - i*epsilon; nu real; k=0,1,2,...; L2Nf3={Bar Lambda}^2_{Nf=3} (all scales in GeV^2).";

  (*------------------------------------------------------------*)
  (* Fractional Analytic Perturbation Theory (FAPT): *)
AFAPT1l::usage = "AFAPT1l[Nf, nu, k, Q2, L2, Fi] gives the 1-loop FAPT coupling with the inclusion of logarithm of the coupling in the form [a1l^nu Log[a1l]^k]_{an.FAPT}, where [...]_{an.FAPT} is the FAPT-analytization of the term in [...]; Nf is the number of active flavors; nu real index; k=0,1,2,...; Q2=|Q^2|; L2={Bar Lambda}^2_{Nf} (both Q2 and L2 in GeV^2); Fi=Arg[Q^2], i.e., Q^2 = Q2 Exp[I*Fi].";
AFAPT2l::usage = "AFAPT2l[Nf, nu, k, Q2, L2, Fi] gives the 2-loop FAPT coupling with the inclusion of logarithm of the coupling in the form [a2l^nu Log[a2l]^k]_{an.FAPT}, where [...]_{an.FAPT} is the FAPT-analytization of the term in [...]; Nf is the number of active flavors; nu real index; k=0,1,2,...; Q2=|Q^2|; L2={Bar Lambda}^2_{Nf} (both Q2 and L2 in GeV^2); Fi=Arg[Q^2], i.e., Q^2 = Q2 Exp[I*Fi].";
AFAPT3l::usage = "AFAPT3l[Nf, nu, k, Q2, L2, Fi] gives the 3-loop FAPT coupling with the inclusion of logarithm of the coupling in the form [a3l^nu Log[a3l]^k]_{an.FAPT}, where [...]_{an.FAPT} is the FAPT-analytization of the term in [...]; Nf is the number of active flavors; nu real index; k=0,1,2,...; Q2=|Q^2|; L2={Bar Lambda}^2_{Nf} (both Q2 and L2 in GeV^2); Fi=Arg[Q^2], i.e., Q^2 = Q2 Exp[I*Fi].";
AFAPT4l::usage = "AFAPT4l[Nf, nu, k, Q2, L2, Fi] gives the 4-loop FAPT coupling with the inclusion of logarithm of the coupling in the form [a4l^nu Log[a4l]^k]_{an.FAPT}, where [...]_{an.FAPT} is the FAPT-analytization of the term in [...]; Nf is the number of active flavors; nu real index; k=0,1,2,...; Q2=|Q^2|; L2={Bar Lambda}^2_{Nf} (both Q2 and L2 in GeV^2); Fi=Arg[Q^2], i.e., Q^2 = Q2 Exp[I*Fi].";

  (*------------------------------------------------------------*)
  (* FAPT global: *)
AFAPT1lglob::usage = "AFAPT1lglob[nu, k, Q2, L2nf3, fi] gives the 1-loop global FAPT coupling with the inclusion of logarithm of the coupling in the form [a1lglob^nu Log[a1lglob]^k]_{an.FAPT},  where [...]_{an.FAPT} is the FAPT-analytization of the term in [...]; L2nf3 is the (1-loop) squared Lambda scale (in GeV^2) at Nf=3 active quark flavors; other notations as in AFAPT1l.";
AFAPT2lglob::usage = "AFAPT2lglob[nu, k, Q2, L2nf3, fi] gives the 2-loop global FAPT coupling with the inclusion of logarithm of the coupling in the form [a2lglob^nu Log[a2lglob]^k]_{an.FAPT},  where [...]_{an.FAPT} is the FAPT-analytization of the term in [...]; L2nf3 is the (2-loop) squared Lambda scale (in GeV^2) at Nf=3 active quark flavors; other notations as in AFAPT1l.";
AFAPT3lglob::usage = "AFAPT3lglob[nu, k, Q2, L2nf3, fi] gives the 3-loop global FAPT coupling with the inclusion of logarithm of the coupling in the form [a3lglob^nu Log[a3lglob]^k]_{an.FAPT},  where [...]_{an.FAPT} is the FAPT-analytization of the term in [...]; L2nf3 is the (3-loop) squared Lambda scale (in GeV^2) at Nf=3 active quark flavors; other notations as in AFAPT1l.";
AFAPT4lglob::usage =  "AFAPT4lglob[nu, k, Q2, L2nf3, fi] gives the 4-loop global FAPT coupling with the inclusion of logarithm of the coupling in the form [a4lglob^nu Log[a4lglob]^k]_{an.FAPT},  where [...]_{an.FAPT} is the FAPT-analytization of the term in [...]; L2nf3 is the (4-loop) squared Lambda scale (in GeV^2) at Nf=3 active quark flavors; other notations as in AFAPT1l.";

  (*------------------------------------------------------------*)
  (* 2-delta analytic QCD (2d): *)
tA2d::usage = "tA2d[Nf, nu, Q2, Fi] gives the 2danQCD tilde coupling {tilde A}_{nu}(Q^2;Nf) in the c2-Lambert scheme [the underlying coupling is a2d of Eq.(14) with c2=c22din] at Q^2=|Q^2|*Exp[I Fi], with Q2=|Q^2| (in GeV^2), and Nf=3,4,5,6 is the number of active quark flavors.";

A2d1l::usage = "A2d1l[Nf, 0, nu, Q2, Fi] gives the N=1 ('1-loop') 2danQCD coupling A_{nu}(Q^2;Nf), i.e.,  A_{nu}={tilde A}_{nu} is the 1-loop 2danQCD analog of the pQCD power a(Q^2;Nf)^{nu};  Q2=|Q^2| in GeV^2; Q^2=Q2*Exp[I Fi]; nu is any real index (-1 < nu)."; 
A2d2l::usage = "A2d2l[Nf, n, nu, Q2, Fi] gives the N=2 ('2-loop') 2danQCD coupling A_{nu+n}(Q^2;Nf), the 2danQCD analog of the pQCD power a(Q^2;Nf)^{nu+n};  Q2=|Q^2| in GeV^2; Q^2=Q2*Exp[I Fi]; nu is any real index (-1 < nu), and n is an integer: n=0,1; A_{nu+n} is constructed from {tilde A}_{nu+n+m}'s using Eq.(28) of the manuscript with N=2 (and nu=nu_0): 
A2d2l[Nf,n,nu,Q2,Fi]=tA2d[Nf,nu+n,Q2,Fi]+ Sum(tk_m(nu+n)*tA2d[Nf,nu+n+m,Q2,Fi]) where the Sum runs from m=1 to m=N-1-n (where:N=2).";
A2d3l::usage = "A2d3l[Nf, n, nu, Q2, Fi] gives the N=3 ('3-loop') 2danQCD coupling A_{nu+n}(Q^2;Nf), the 2danQCD analog of the pQCD power a(Q^2;Nf)^{nu+n};  Q2=|Q^2| in GeV^2; Q^2=Q2*Exp[I Fi]; nu is any real index (-1 < nu), and n is an integer: n=0,1,2; A_{nu+n} is constructed from {tilde A}_{nu+n+m}'s using Eq.(28) of the manuscript with N=3 (and nu=nu_0): 
A2d3l[Nf,n,nu,Q2,Fi]=tA2d[Nf,nu+n,Q2,Fi]+ Sum(tk_m(nu+n)*tA2d[Nf,nu+n+m,Q2,Fi]) where the Sum runs from m=1 to m=N-1-n (where:N=3)."; 
A2d4l::usage = "A2d4l[Nf, n, nu, Q2, Fi] gives the N=4 ('4-loop') 2danQCD coupling A_{nu+n}(Q^2;Nf), the 2danQCD analog of the pQCD power a(Q^2;Nf)^{nu+n};  Q2=|Q^2| in GeV^2; Q^2=Q2*Exp[I Fi]; nu is any real index (-1 < nu), and n is an integer: n=0,1,2,3; A_{nu+n} is constructed from {tilde A}_{nu+n+m}'s using Eq.(28) of the manuscript with N=4 (and nu=nu_0): 
A2d3l[Nf,n,nu,Q2,Fi]=tA2d[Nf,nu+n,Q2,Fi]+ Sum(tk_m(nu+n)*tA2d[Nf,nu+n+m,Q2,Fi]) where the Sum runs from m=1 to m=N-1-n (where:N=4)." ;
A2d5l::usage = "A2d5l[Nf, n, nu, Q2, Fi] gives the N=5 ('5-loop') 2danQCD coupling A_{nu+n}(Q^2;Nf), the 2danQCD analog of the pQCD power a(Q^2;Nf)^{nu+n};  Q2=|Q^2| in GeV^2; Q^2=Q2*Exp[I Fi]; nu is any real index (-1 < nu), and n is an integer: n=0,1,2,3,4; A_{nu+n} is constructed from {tilde A}_{nu+n+m}'s using Eq.(28) of the manuscript with N=5 (and nu=nu_0): 
A2d3l[Nf,n,nu,Q2,Fi]=tA2d[Nf,nu+n,Q2,Fi]+ Sum(tk_m(nu+n)*tA2d[Nf,nu+n+m,Q2,Fi]) where the Sum runs from m=1 to m=N-1-n (where: N=5)."; 

  (*------------------------------------------------------------*)
  (* Massive perturbation theory (MPT): *)
tAMPT1l::usage = "tAMPT1l[Nf, n, del, Q2, M2, L2]= {tilde A}_{nu=n+del}(Q^2=Q2) in 1-loop MPT with m_{gl}^2=M2 (in GeV^2), and MSbar squared Lambda_{Nf}^2=L2; where Q2 is COMPLEX Q2=Q^2=|Q^2| Exp(I fi) (in GeV^2); nu=n+del is the index of generalized logarithmic derivative, where n=0,1,3,4;  0=<del<1.";
tAMPT2l::usage = "tAMPT2l[Nf, n, del, Q2, M2, L2]= {tilde A}_{nu=n+del}(Q^2=Q2) in 2-loop MPT with m_{gl}^2=M2 (in GeV^2), and MSbar squared Lambda_{Nf}^2=L2; where Q2 is COMPLEX Q2=Q^2=|Q^2| Exp(I fi) (in GeV^2); nu=n+del is the index of generalized logarithmic derivative, where n=0,1,3,4;  0=<del<1.";
tAMPT3l::usage = "tAMPT3l[Nf, n, del, Q2, M2, L2]= {tilde A}_{nu=n+del}(Q^2=Q2) in 3-loop MPT with m_{gl}^2=M2 (in GeV^2), and MSbar squared Lambda_{Nf}^2=L2; where Q2 is COMPLEX Q2=Q^2=|Q^2| Exp(I fi) (in GeV^2); nu=n+del is the index of generalized logarithmic derivative, where n=0,1,3,4;  0=<del<1.";
tAMPT4l::usage = "tAMPT4l[Nf, n, del, Q2, M2, L2]= {tilde A}_{nu=n+del}(Q^2=Q2) in 4-loop MPT with m_{gl}^2=M2 (in GeV^2), and MSbar squared Lambda_{Nf}^2=L2; where Q2 is COMPLEX Q2=Q^2=|Q^2| Exp(I fi) (in GeV^2); nu=n+del is the index of generalized logarithmic derivative, where n=0,1,3,4;  0=<del<1.";

AMPT1l::usage = "AMPT1l[Nf,Nu,Q2,M2,L2] gives the N=1 (1-loop) MPT coupling A_{nu}(Q^2;Nf), the MPT analog of the pQCD power a(Q^2;Nf)^{nu};  Q2=Q^2 in GeV^2 (complex number in general); M2 is the MPT scale (in GeV^2); 0<Nu<5; A_{nu} is constructed from MPT 1-loop {tilde A}_{nu+m}'s (tAA1l) using Eq.(28) of the manuscript with truncations as explained in the text at the end of Subsection 3.4 of the manuscript.";
AMPT2l::usage = "AMPT2l[Nf,Nu,Q2,M2,L2] gives the N=2 (2-loop) MPT coupling A_{nu}(Q^2;Nf), the MPT analog of the pQCD power a(Q^2;Nf)^{nu};  Q2=Q^2 in GeV^2 (complex number in general); M2 is the MPT scale (in GeV^2); 0<Nu<5; A_{nu} is constructed from {tilde A}_{nu+m}'s using Eq.(28) with truncations as explained in the text at the end of Subsection 3.4 of the manuscript.";
AMPT3l::usage = "AMPT3l[Nf,Nu,Q2,M2,L2] gives the N=3 (3-loop) MPT coupling A_{nu}(Q^2;Nf), the MPT analog of the pQCD power a(Q^2;Nf)^{nu};  Q2=Q^2 in GeV^2 (complex number in general); M2 is the MPT scale (in GeV^2); 0<Nu<5; A_{nu} is constructed from {tilde A}_{nu+m}'s using Eq.(28) with truncations as explained in the text at the end of Subsection 3.4 of the manuscript.";
AMPT4l::usage = "AMPT4l[Nf,Nu,Q2,M2,L2] gives the N=4 (4-loop) MPT coupling A_{nu}(Q^2;Nf), the MPT analog of the pQCD power a(Q^2;Nf)^{nu};  Q2=Q^2 in GeV^2 (complex number in general); M2 is the MPT scale (in GeV^2); 0<Nu<5; A_{nu} is constructed from {tilde A}_{nu+m}'s using Eq.(28) with truncations as explained in the text at the end of Subsection 3.4  of the manuscript.";
  (*------------------------------------------------------------*)


(* ::Section:: *)
(*Numerical constants: Inputs*)


(* Parameters specifically for 2danQCD model *)
(* You can choose for c22din any number between -6 and -2; c2dl=-4.9 is the preferred (default) case *) 

<<s0r.m; 
(* Li__nu.m gives a stable version of polylog[nu,u] *)
<<Li__nu.m;
c22din=-4.9;
$s0[3]=s0r[c22din];
$s0[4]=$s0[3]; 
$s0[5]=$s0[3]; $s0[6]=$s0[3];
(* s0=23.06; f1ss0[3]=0.271312; f2ss0[3]=0.5409; s1s0[3]=16.837; s2s0[3]=0.8077; *)
(* f1ss0[4]=0.34812; f2ss0[4]=0.761416; s1s0[4]=16.6753; s2s0[4]=0.130055; *)

(* parameters such as limits in numerical calculus *)

$lnsm=Log[10.^40];
(* $wwmax=10.^12; This is too high in Mathematica 11.1 for FAPT couploing calculations *)
$wwmax=10^12;
$wmax=10^8;
$wmax2d=10.^3;
$wmaxJ=10.^8.5;
$small=10.^-10;
$eepss=10.^-7;
$Nloo=8;


(* First pQCD: barMS, at 5-loop level, to be used at nf=3 and nf=4: *)

kappa=2;
(* The value of kappa can be changed, 1 =< kappa =< 3 *)
kappa2d=2;
(* The value of kappa2=2 should not be changed *)

NumDefanQCD = {
MZ -> 91.1876,
mc -> 1.27, 
mb -> 4.20, 
mt -> 163.0
};

NumDefanQCD1 = {
Q6 -> kappa mt/. NumDefanQCD,(*Q6=mu_thr(ng=5\[Rule]6)*)
Q5 -> kappa mb /. NumDefanQCD,(*Q5=mu_thr(nf=5\[Rule]4)*)
Q4 -> kappa mc /. NumDefanQCD,(*Q4=mu_thr(nf=4\[Rule]3)*)
aaMZ -> 0.1184/Pi, 
Q2thr[5] -> (kappa*mt)^2 /. NumDefanQCD, 
Q2thr[4] -> (kappa*mb)^2 /. NumDefanQCD, 
Q2thr[3] -> (kappa*mc)^2 /. NumDefanQCD
};


(* ::Section:: *)
(*\beta function coefficient in MSbar scheme, with effective *)
(*three-loop Gardi et al. beta function *)

Needs["NumericalCalculus`"]

$be0[nq_]:=(11 -2 nq/3)/4;
$be1[nq_]:=(102-38 nq/3)/16;
$be2MS[nq_]:=(2857/2 - 5033 nq/18 + 325 nq^2/54)/64;
$be3MS[nq_]:=(1/256)*( (149753/6 + 3564 Zeta[3]) - ((1078361/162) + (6508/27) Zeta[3])*nq + ((50065/162) + (6472/81) Zeta[3])*nq^2 + (1093/729)*nq^3 );
$c1[nq_]:=(102-38 nq/3)/(11 -2 nq/3)/4;
$c2MS[nq_]:=$be2MS[nq]/$be0[nq];
$c3MS[nq_]:=$be3MS[nq]/$be0[nq];

(* RSch we use is: c-Lambert RSch: *)
$c22d[nq_]:=c22din; $c32d[nq_]:=$c22d[nq]^2/$c1[nq]; $c42d[nq_]:=$c22d[nq]^3/$c1[nq]^2;
$c52d[nq_]:=$c22d[nq]^4/$c1[nq]^3; $c62d[nq_]:=$c22d[nq]^5/$c1[nq]^4;

(* betaP11 is c2-Lambert RSch beta function: *)
$betaP11[x_,nq_] :=  - $be0[nq] x^2 (1 + ($c1[nq] - $c22d[nq]/$c1[nq]) x)/(1 - ($c22d[nq]/$c1[nq]) x);
$btP11[x_,nq_]:=$betaP11[x,nq];
$betaMS4l[x_,nq_]:=- $be0[nq] x^2 (1 + $c1[nq] x + $c2MS[nq] x^2 + $c3MS[nq] x^3);

(* ::Subsection:: *)
(* 4-loop MSbar coefficient *)


K[1,0,nq_]:=1/$be0[nq]; 

K[2,0,nq_]:=0; K[2,1,nq_]:=-$c1[nq]/$be0[nq]^2; 

K[3,0,nq_]:=(-$c1[nq]^2 + $c2)/$be0[nq]^3; 
K[3,1,nq_]:= - $c1[nq]^2/$be0[nq]^3; 
K[3,2,nq_]:=-K[3,1,nq];

K[4,0,nq_]:=-(1/2) ($c1[nq]^3/$be0[nq]^4) (1 - $c3/$c1[nq]^3);
K[4,1,nq_]:=-($c1[nq]^3/$be0[nq]^4)(-2 + 3 $c2/$c1[nq]^2);
K[4,2,nq_]:=(5/2)($c1[nq]^3/$be0[nq]^4); 
K[4,3,nq_]:=-($c1[nq]^3/$be0[nq]^4);

K[5,0,nq_]:=(1/6/$be0[nq]^5)(7 $c1[nq]^4 - 18 $c1[nq]^2 $c2 + 10 $c2^2 - $c1[nq] $c3 + 2 $c4);
K[5,1,nq_]:=($c1[nq]/$be0[nq]^5)(4 $c1[nq]^3 - 3 $c1[nq] $c2 - 2 $c3);
K[5,2,nq_]:=-(3/2/$be0[nq]^5)($c1[nq]^4 - 4 $c1[nq]^2 $c2);
K[5,3,nq_]:=-(13/3) ($c1[nq]^4/$be0[nq]^5);  
K[5,4,nq_]:=$c1[nq]^4/$be0[nq]^5;

K[6, 0,nq_] := (17*$c1[nq]^5 - 18*$c1[nq]^3*$c2 - $c1[nq]*$c2^2 - 23*$c1[nq]^2*$c3 + 24*$c2*$c3 - 2*$c1[nq]*$c4 + 3*$c5)/(12*$be0[nq]^6);
K[6, 1,nq_] := (-11*$c1[nq]^5 + 72*$c1[nq]^3*$c2 - 50*$c1[nq]*$c2^2 - 7*$c1[nq]^2*$c3 - 10*$c1[nq]*$c4)/(6*$be0[nq]^6);
K[6, 2,nq_] := (-23*$c1[nq]^5 + 27*$c1[nq]^3*$c2 + 10*$c1[nq]^2*$c3)/(2*$be0[nq]^6);
K[6, 3,nq_] := (-11*$c1[nq]^5 - 60*$c1[nq]^3*$c2)/(6*$be0[nq]^6);
K[6, 4,nq_] := (77*$c1[nq]^5)/(12*$be0[nq]^6);
K[6, 5,nq_] := -($c1[nq]^5/$be0[nq]^6);

K[7, 0,nq_] := (-62*$c1[nq]^6 + 390*$c1[nq]^4*$c2 - 497*$c1[nq]^2*$c2^2 + 192*$c2^3 - 13*$c1[nq]^3*$c3 - 48*$c1[nq]*$c2*$c3 + 33*$c3^2 - 94*$c1[nq]^2*$c4 + 96*$c2*$c4 - 9*$c1[nq]*$c5 + 12*$c6)/(60*$be0[nq]^7);
K[7, 1,nq_] := (-62*$c1[nq]^6 + 126*$c1[nq]^4*$c2 - 47*$c1[nq]^2*$c2^2 + 62*$c1[nq]^3*$c3 - 72*$c1[nq]*$c2*$c3 - 4*$c1[nq]^2*$c4 - 9*$c1[nq]*$c5)/(6*$be0[nq]^7);
K[7, 2,nq_] := (-12*$c1[nq]^6 - 45*$c1[nq]^4*$c2 + 50*$c1[nq]^2*$c2^2 + 17*$c1[nq]^3*$c3 + 10*$c1[nq]^2*$c4)/(2*$be0[nq]^7);
K[7, 3,nq_] := ($c1[nq]^3*(127*$c1[nq]^3 - 222*$c1[nq]*$c2 - 60*$c3))/(6*$be0[nq]^7);
K[7, 4,nq_] := (5*(11*$c1[nq]^6 + 18*$c1[nq]^4*$c2))/(6*$be0[nq]^7);
K[7, 5,nq_] := (-87*$c1[nq]^6)/(10*$be0[nq]^7);
K[7, 6,nq_] := $c1[nq]^6/$be0[nq]^7;

K[8, 0,nq_]:= (-349*$c1[nq]^7 + 870*$c1[nq]^5*$c2 - 444*$c1[nq]^3*$c2^2 - 46*$c1[nq]*$c2^3 + 634*$c1[nq]^4*$c3 - 1386*$c1[nq]^2*$c2*$c3 + 760*$c2^2*$c3 - 49*$c1[nq]*$c3^2 + 12*$c1[nq]^3*$c4 - 108*$c1[nq]*$c2*$c4 + 100*$c3*$c4 - 168*$c1[nq]^2*$c5 + 170*$c2*$c5 - 16*$c1[nq]*$c6 + 20*$c7)/(120*$be0[nq]^8)
K[8, 1,nq_]:= -($c1[nq]*(62*$c1[nq]^6 + 490*$c1[nq]^4*$c2 - 1003*$c1[nq]^2*$c2^2 + 448*$c2^3 - 237*$c1[nq]^3*$c3 + 128*$c1[nq]*$c2*$c3 + 77*$c3^2 - 206*$c1[nq]^2*$c4 + 224*$c2*$c4 + 9*$c1[nq]*$c5 + 28*$c6))/(20*$be0[nq]^8);
K[8, 2,nq_]:= ($c1[nq]^2*(362*$c1[nq]^5 - 1152*$c1[nq]^3*$c2 + 629*$c1[nq]*$c2^2 - 332*$c1[nq]^2*$c3 + 504*$c2*$c3 + 88*$c1[nq]*$c4 + 63*$c5))/(12*$be0[nq]^8);
K[8, 3,nq_]:= ($c1[nq]^3*(211*$c1[nq]^4 + 93*$c1[nq]^2*$c2 - 350*$c2^2 - 179*$c1[nq]*$c3 - 70*$c4))/(6*$be0[nq]^8);
K[8, 4,nq_]:= (-223*$c1[nq]^7 + 638*$c1[nq]^5*$c2 + 140*$c1[nq]^4*$c3)/(8*$be0[nq]^8);
K[8, 5,nq_]:= (-323*$c1[nq]^7 - 315*$c1[nq]^5*$c2)/(15*$be0[nq]^8);
K[8, 6,nq_]:= (223*$c1[nq]^7)/(20*$be0[nq]^8);
K[8, 7,nq_]:= -($c1[nq]^7/$be0[nq]^8);

(* ::Section:: *)
(* QCD perturbative running couplings: a1l, a2l, a3l and a4l,
   for one, two, three and four-loop respectively *)

(* with u = |Q^2|/L2>0;  Q^2 = |Q^2| Exp[I*fi]
note that L2=LLamb^2 in Lambert scheme *) 

beta[nq_,x_] :=  - $be0[nq] x^2 (1 + ($c1[nq] - $c2/$c1[nq]) x)/(1 - ($c2/$c1[nq]) x);
 betared[nq_,x_]:= x^2 (1 + ($c1[nq] - $c2/$c1[nq]) x)/(1 - ($c2/$c1[nq]) x);

 zu[nq_,u_, fi_] := - (1/$c1[nq]/E) u^(-$be0[nq]/$c1[nq]) Exp[- I*($be0[nq]/$c1[nq]) fi];

 a[nq_,u_, fi_] := If[fi >= 0, -(1/$c1[nq]) 1/(1 -$c2/$c1[nq]^2 + ProductLog[-1, zu[nq,u, fi]]), -(1/$c1[nq]) 1/(1 - $c2/$c1[nq]^2 + ProductLog[1, zu[nq,u, fi]])];

a1l[nq_,u_,fi_]:= 1/$be0[nq]/Log[uu]/.uu->u Exp[I fi];
a1l[nq_, Q2_, L2_, fi_] := a1l[nq,Q2/L2,fi]

(* a1l[nq_,u_, fi_]:=a[nq,u, fi]/.{$c1[nq]->0,$c2->0} *)
(* We relate the two- and three-loop Lambda^2 Lambert L2L2l and L2L3l  
   with Lambda^2 MSbar L2 (our input) *)

L2L2l[L2_, nq_] := L2Lx /. FindRoot[a2l[nq, 10^10/L2Lx, 0] == a2lLog[nq, 10^10/L2], {L2Lx, 0.1, 0.6}];
L2L3l[L2_, nq_] := L3Lx /. FindRoot[a3l[nq, 10^10/L3Lx, 0] == a3lLog[nq, 10^10/L2], {L3Lx, 0.1, 0.6}];

a2lLog[nq_, u_] := N[Sum[K[i, j, nq] Log[Log[u]]^j/Log[u]^i, {i, 1, $Nloo}, {j, 0, i - 1}]] /. {$c2 -> 0, $c3 -> 0, $c4 -> 0, $c5 -> 0, $c6 -> 0, $c7 -> 0}
a2l[nq_,u_, fi_]:=a[nq,u, fi]/.{$c2->0}
a2l[nq_, Q2_, L2_, fi_] := a2l[nq, Q2/L2L2l[L2, nq], fi]

  a3lLog[nq_, u_] := N[Sum[K[i, j, nq] Log[Log[u]]^j/Log[u]^i, {i, 1, $Nloo}, {j, 0, i - 1}]] /. {$c2 -> $c2MS[nq], $c3 -> $c2MS[nq]^2/$c1[nq], $c4 -> $c2MS[nq]^3/$c1[nq]^2, $c5 -> $c2MS[nq]^4/$c1[nq]^3, $c6 -> $c2MS[nq]^5/$c1[nq]^4, $c7 -> $c2MS[nq]^6/$c1[nq]^5}
a3l[nq_,u_, fi_]:=a[nq,u, fi]/.{$c2->$c2MS[nq]}
a3l[nq_, Q2_, L2_, fi_] := a3l[nq, Q2/L2L3l[L2, nq], fi]
a2d[nq_,u_, fi_]:=a[nq,u, fi]/.{$c2->c22din}
a2d[nq_, Q2_, L2_, fi_]:=a2d[nq,Q2/L2L3l[L2, nq], fi]

  a4l[nq_,u_, fi_]:=N[Sum[K[i,j,nq]Log[Log[uu]]^j/Log[uu]^i,{i,1,$Nloo},{j,0,i-1}]]/.{uu->u Exp[I fi], $c2 -> $c2MS[nq], $c3->$c3MS[nq], $c4 -> 0, $c5 -> 0, $c6 -> 0, $c7 -> 0}
a4l[nq_, Q2_, L2_, fi_]:= a4l[nq,Q2/L2, fi]


(* ::Section:: *)
(* Global perturbative coupling *)

(* ::Subsection:: *)
(* we find, for input nf value, the (approximate) value of L2=Lambda(nf)^2 that appears in a[u=Q2/L2] for this nf:  *)
(* The 3-loop matching UPWARDS *)

lh=2*Log[kappa];

xi3=1.202057;
ct2=11/72;
ct3[nff_]:=- 82043*xi3/27648 + 564731/124416 - 2633*(nff)/31104;
k1= -lh/6;
k2= ( lh^2/36 - 19*lh/24 + ct2);
k3[nff_]:=-lh^3/216 - 131*lh^2/576 + lh*(-6793 + 281*(nff))/1728 +ct3[nff];

x1=-k1;
x2=-k2+2*k1^2;
x3[nff_]:=-k3[nff]+5*k2*k1-5*k1^3;

(* input: L2klnf3=Lambda[nf=3]^2 to k-loop order *)
On[General::stop];

L21l[4,L21lnf3_]:=L24/.FindRoot[a1l[4,Q2thr[3]/.NumDefanQCD1,L24,0]==a1l[3,Q2thr[3]/.NumDefanQCD1,L21lnf3,0],{L24,0.2,0.4}];
L22l[4,L22lnf3_]:=L24/.FindRoot[a2lLog[4,(Q2thr[3]/.NumDefanQCD1)/L24]==a2lLog[3,(Q2thr[3]/.NumDefanQCD1)/L22lnf3]*(1 + x1*a2lLog[3,(Q2thr[3]/.NumDefanQCD1)/L22lnf3] ),{L24,0.2,0.4}];
L23l[4,L23lnf3_]:=L24/.FindRoot[a3lLog[4,(Q2thr[3]/.NumDefanQCD1)/L24]==a3lLog[3,(Q2thr[3]/.NumDefanQCD1)/L23lnf3]*(1 + x1*a3lLog[3,(Q2thr[3]/.NumDefanQCD1)/L23lnf3] + x2*a3lLog[3,(Q2thr[3]/.NumDefanQCD1)/L23lnf3]^2),{L24,0.2,0.4}];
L24l[4,L24lnf3_]:=L24/.FindRoot[a4l[4,Q2thr[3]/.NumDefanQCD1,L24,0]==a4l[3,Q2thr[3]/.NumDefanQCD1,L24lnf3,0]*(1 + x1*a4l[3,Q2thr[3]/.NumDefanQCD1,L24lnf3,0] + x2*a4l[3,Q2thr[3]/.NumDefanQCD1,L24lnf3,0]^2 + x3[3]*a4l[3,Q2thr[3]/.NumDefanQCD1,L24lnf3,0]^3),{L24,0.2,0.4}];

L21l[5,L21lnf3_]:=L25/.FindRoot[a1l[5,Q2thr[4]/.NumDefanQCD1,L25,0]==a1l[4,Q2thr[4]/.NumDefanQCD1,L21l[4,L21lnf3],0],{L25,0.1,0.3}];
L22l[5,L22lnf3_]:=L25/.FindRoot[a2lLog[5,(Q2thr[4]/.NumDefanQCD1)/L25]==a2lLog[4,(Q2thr[4]/.NumDefanQCD1)/L22l[4,L22lnf3]]*(1 + x1*a2lLog[4,(Q2thr[4]/.NumDefanQCD1)/L22l[4,L22lnf3]]) ,{L25,0.1,0.3}];
L23l[5,L23lnf3_]:=L25/.FindRoot[a3lLog[5,(Q2thr[4]/.NumDefanQCD1)/L25]==a3lLog[4,(Q2thr[4]/.NumDefanQCD1)/L23l[4,L23lnf3]]*(1 + x1*a3lLog[4,(Q2thr[4]/.NumDefanQCD1)/L23l[4,L23lnf3]] + x2*a3lLog[4,(Q2thr[4]/.NumDefanQCD1)/L23l[4,L23lnf3]]^2),{L25,0.1,0.3}];
L24l[5,L24lnf3_]:=L25/.FindRoot[a4l[5,Q2thr[4]/.NumDefanQCD1,L25,0]==a4l[4,Q2thr[4]/.NumDefanQCD1,L24l[4,L24lnf3],0]*(1 + x1*a4l[4,Q2thr[4]/.NumDefanQCD1,L24l[4,L24lnf3],0] + x2*a4l[4,Q2thr[4]/.NumDefanQCD1,L24l[4,L24lnf3],0]^2 + x3[4]*a4l[4,Q2thr[4]/.NumDefanQCD1,L24l[4,L24lnf3],0]^3),{L25,0.1,0.3}];

L21l[6,L21lnf3_]:=L26/.FindRoot[a1l[6,Q2thr[5]/.NumDefanQCD1,L26,0]==a1l[5,Q2thr[5]/.NumDefanQCD1,L21l[5,L21lnf3],0],{L26,0.05,0.2}];
L22l[6,L22lnf3_]:=L26/.FindRoot[a2lLog[6,(Q2thr[5]/.NumDefanQCD1)/L26]==a2lLog[5,(Q2thr[5]/.NumDefanQCD1)/L22l[5,L22lnf3]]*(1 + x1*a2lLog[5,(Q2thr[5]/.NumDefanQCD1)/L22l[5,L22lnf3]]),{L26,0.05,0.2}];
L23l[6,L23lnf3_]:=L26/.FindRoot[a3lLog[6,(Q2thr[5]/.NumDefanQCD1)/L26]==a3lLog[5,(Q2thr[5]/.NumDefanQCD1)/L23l[5,L23lnf3]]*(1 + x1*a3lLog[5,(Q2thr[5]/.NumDefanQCD1)/L23l[5,L23lnf3]] + x2*a3lLog[5,(Q2thr[5]/.NumDefanQCD1)/L23l[5,L23lnf3]]^2),{L26,0.05,0.2}];
L24l[6,L24lnf3_]:=L26/.FindRoot[a4l[6,Q2thr[5]/.NumDefanQCD1,L26,0]==a4l[5,Q2thr[5]/.NumDefanQCD1,L24l[5,L24lnf3],0]*(1 + x1*a4l[5,Q2thr[5]/.NumDefanQCD1,L24l[5,L24lnf3],0] + x2*a4l[5,Q2thr[5]/.NumDefanQCD1,L24l[5,L24lnf3],0]^2 + x3[5]*a4l[5,Q2thr[5]/.NumDefanQCD1,L24l[5,L24lnf3],0]^3),{L26,0.05,0.2}];



(* ::Subsection:: *)
(* Global QCD perturbative running couplings: a1lglob, a2lglob, a3lglob and a4lglob,
   for one, two, three and four-loop respectively *)

a1lglob[Q2_,L2nf3_,fi_]:=If[Abs[Q2]<(Q2thr[3]/.NumDefanQCD1),a1l[3,Q2,L2nf3,fi],If[Abs[Q2]<(Q2thr[4]/.NumDefanQCD1),a1l[4,Q2,L21l[4,L2nf3],fi],If[Abs[Q2]<(Q2thr[5]/.NumDefanQCD1),a1l[5,Q2,L21l[5,L2nf3],fi],a1l[6,Q2,L21l[6,L2nf3],fi]]]];

a2lglob[Q2_,L2nf3_,fi_]:=If[Abs[Q2]<(Q2thr[3]/.NumDefanQCD1),a2l[3,Q2,L2nf3,fi],If[Abs[Q2]<(Q2thr[4]/.NumDefanQCD1),a2l[4,Q2,L21l[4,L2nf3],fi],If[Abs[Q2]<(Q2thr[5]/.NumDefanQCD1),a2l[5,Q2,L22l[5,L2nf3],fi],a2l[6,Q2,L22l[6,L2nf3],fi]]]];

a3lglob[Q2_,L2nf3_,fi_]:=If[Abs[Q2]<(Q2thr[3]/.NumDefanQCD1),a3l[3,Q2,L2nf3,fi],If[Abs[Q2]<(Q2thr[4]/.NumDefanQCD1),a3l[4,Q2,L23l[4,L2nf3],fi],If[Abs[Q2]<(Q2thr[5]/.NumDefanQCD1),a3l[5,Q2,L23l[5,L2nf3],fi],a3l[6,Q2,L23l[6,L2nf3],fi]]]];

a4lglob[Q2_,L2nf3_,fi_]:=If[Abs[Q2]<(Q2thr[3]/.NumDefanQCD1),a4l[3,Q2,L2nf3,fi],If[Abs[Q2]<(Q2thr[4]/.NumDefanQCD1),a4l[4,Q2,L24l[4,L2nf3],fi],If[Abs[Q2]<(Q2thr[5]/.NumDefanQCD1),a4l[5,Q2,L24l[5,L2nf3],fi],a4l[6,Q2,L24l[6,L2nf3],fi]]]];


(* ::Section:: *)
(* The spectral function for perturbative coupling *)
(* tr[nq,nu,k,s] = Im[a^nu Log[a]^k at [u,fi] = [s,-Pi]; with nf=nq *)


sL[nq_]:=$c1[nq]^(-$c1[nq]/$be0[nq]);
tr[nq_,nu_,k_,s_]:=If[s>0,Im[a[nq,s,-Pi]^nu Log[a[nq,s,-Pi]]^k],0];

 (* tr1l[nq_,nu_,k_,sig_,L2_]:=If[sig>0,Im[a1l[nq,sig,L2,-Pi+$eepss]^nu Log[a1l[nq,sig,L2,-Pi+$eepss]]^k],0]; *)
 tr1l[nq_,nu_,k_,sig_,L2_]:=If[sig>0,(1/$be0[nq])^nu Im[(Log[sig/L2] + I Pi)^nu]/(Log[sig/L2]^2+Pi^2)^nu,0];

tr2l[nq_,nu_,k_,sig_,L2_]:=If[sig>0,Im[a2l[nq,sig,L2,-Pi]^nu Log[a2l[nq,sig,L2,-Pi]]^k],0];

tr3l[nq_,nu_,k_,sig_,L2_]:=If[sig>0,Im[a3l[nq,sig,L2,-Pi]^nu Log[a3l[nq,sig,L2,-Pi]]^k],0];
tr12d[nq_,s_]:= If[s!=0,If[s > 0, Im[a2d[nq,s, -Pi]], Im[(1/$c1[nq])  1/(1 - $c22d[nq]/$c1[nq]^2 + ProductLog[-1, zu[nq,-s, 0]])]],0];
tr1x[nq_,xp_]:=tr12d[nq,Exp[xp]];

tr4l[nq_,nu_,k_,sig_,L2_]:=If[sig>0,Im[a4l[nq,sig,L2,-Pi+$eepss]^nu Log[a4l[nq,sig,L2,-Pi+$eepss]]^k],0];



(* ::Subsection:: *)
(* The global spectral function for perturbative coupling *)


tr1lglob[nu_,k_,sig_,L2nf3_]:=If[sig>0,Im[a1lglob[sig,L2nf3,-Pi+$eepss]^nu Log[a1lglob[sig,L2nf3,-Pi+$eepss]]^k],0];

tr2lglob[nu_,k_,sig_,L2nf3_]:=If[sig>0,Im[a2lglob[sig,L2nf3,-Pi]^nu Log[a2lglob[sig,L2nf3,-Pi]]^k],0];

tr3lglob[nu_,k_,sig_,L2nf3_]:=If[sig>0,Im[a3lglob[sig,L2nf3,-Pi]^nu Log[a3lglob[sig,L2nf3,-Pi]]^k],0];

tr4lglob[nu_,k_,sig_,L2nf3_]:=If[sig>0,Im[a4lglob[sig,L2nf3,-Pi+$eepss]^nu Log[a4lglob[sig,L2nf3,-Pi+$eepss]]^k],0];



(* ::Section:: *)
(* The general formula for FAPT coupling: *)


(* ::Subsection:: *)
(* A[nq,nu,k,u] = (a^{nu} (ln a)^k (Q2))^{(MA)}_an = AMA[nq,nu,k,u,fi]  *)
(*  (Q2= |Q2| Exp[I fi] = aQ2 Exp[I fi]: (nu is real; k is integer: k=0,1,2,3,... *)

     (* Note: below AFAPT-type of couplings, in the integrals at large w=Log[sig]>100, have exact forms of Q2 (->  IntegerPart[Q2 10^20] 10^-20) anf fi (-> IntegerPart[fi 10^2] 10^-20), in order to avoid instabilities in Mathematica 10.0.1 *)
 AFAPT1l[nq_,nu_,k_,Q2_,L2_,fi_]:=(1/Pi)NIntegrate[tr1l[nq,nu,k,sig,L2]/(sig + Q2 Exp[I fi]),{sig,0,1.}]+(1/Pi)NIntegrate[tr1l[nq,nu,k,Exp[w],L2]/(1 + Q2 Exp[I fi] Exp[-w]),{w,0,100.}]+(1/Pi)NIntegrate[tr1l[nq,nu,k,Exp[w],L2]/(1 + IntegerPart[Q2 10^20] 10^-20 Exp[I IntegerPart[fi 10^2] 10^-20] Exp[-w]),{w,100.,$wwmax}, MaxRecursion->12];
     (* AFAPT1l[nq_,nu_,k_,Q2_,L2_,fi_]:=(1/Pi)NIntegrate[tr1l[nq,nu,k,sig,L2]/(sig + Q2 Exp[I fi]),{sig,0,1.}]+(1/Pi)NIntegrate[tr1l[nq,nu,k,Exp[w],L2]/(1 + Q2 Exp[I fi] Exp[-w]),{w,0,100.}]+(1/Pi)NIntegrate[tr1l[nq,nu,k,Exp[w],L2]/(1 + Q2 Exp[I fi] Exp[-w]),{w,100.,$wwmax}, MaxRecursion->12]; *)

AFAPT2l[nq_,nu_,k_,Q2_,L2_,fi_]:=(1/Pi)NIntegrate[tr2l[nq,nu,k,sig,L2]/(sig + Q2 Exp[I fi]),{sig,0,1.}]+(1/Pi)NIntegrate[tr2l[nq,nu,k,Exp[w],L2]/(1 + Q2 Exp[I fi] Exp[-w]),{w,0,100.}]+(1/Pi)NIntegrate[tr2l[nq,nu,k,Exp[w],L2]/(1 + IntegerPart[Q2 10^20] 10^-20 Exp[I IntegerPart[fi 10^2] 10^-20] Exp[-w]),{w,100.,$wwmax}, MaxRecursion->12];
     (* AFAPT2l[nq_,nu_,k_,Q2_,L2_,fi_]:=(1/Pi)NIntegrate[tr2l[nq,nu,k,sig,L2]/(sig + Q2 Exp[I fi]),{sig,0,1.}]+(1/Pi)NIntegrate[tr2l[nq,nu,k,Exp[w],L2]/(1 + Q2 Exp[I fi] Exp[-w]),{w,0,100.}]+(1/Pi)NIntegrate[tr2l[nq,nu,k,Exp[w],L2]/(1 + Q2 Exp[I fi] Exp[-w]),{w,100.,$wwmax}, MaxRecursion->12]; *)

AFAPT3l[nq_, nu_, k_, Q2_, L2_, fi_] := (1/Pi) NIntegrate[tr3l[nq, nu, k, sig, L2]/(sig + Q2 Exp[I fi]), {sig, 0, 1.}] + (1/Pi) NIntegrate[tr3l[nq, nu, k, Exp[w], L2]/(1 + Q2 Exp[I fi] Exp[-w]), {w, 0, 100.}] + (1/Pi) NIntegrate[tr3l[nq, nu, k, Exp[w], L2]/(1 + IntegerPart[Q2 10^20] 10^-20 Exp[I IntegerPart[fi 10^2] 10^-20] Exp[-w]), {w, 100., $wwmax}, MaxRecursion -> 12];
     (* AFAPT3l[nq_,nu_,k_,Q2_,L2_,fi_]:=(1/Pi)NIntegrate[tr3l[nq,nu,k,sig,L2]/(sig + Q2 Exp[I fi]),{sig,0,1.}]+(1/Pi)NIntegrate[tr3l[nq,nu,k,Exp[w],L2]/(1 + Q2 Exp[I fi] Exp[-w]),{w,0,100.}]+(1/Pi)NIntegrate[tr3l[nq,nu,k,Exp[w],L2]/(1 + Q2 Exp[I fi] Exp[-w]),{w,100.,$wwmax}, MaxRecursion->12]; *)

AFAPT4l[nq_,nu_,k_,Q2_,L2_,fi_]:=(1/Pi)NIntegrate[tr4l[nq,nu,k,sig,L2]/(sig + Q2 Exp[I fi]),{sig,0,1.}]+(1/Pi)NIntegrate[tr4l[nq,nu,k,Exp[w],L2]/(1 + Q2 Exp[I fi] Exp[-w]),{w,0,100.}]+(1/Pi)NIntegrate[tr4l[nq,nu,k,Exp[w],L2]/(1 + IntegerPart[Q2 10^20] 10^-20 Exp[I IntegerPart[fi 10^2] 10^-20] Exp[-w]),{w,100.,$wwmax}, MaxRecursion->12];
     (* AFAPT4l[nq_,nu_,k_,Q2_,L2_,fi_]:=(1/Pi)NIntegrate[tr4l[nq,nu,k,sig,L2]/(sig + Q2 Exp[I fi]),{sig,0,1.}]+(1/Pi)NIntegrate[tr4l[nq,nu,k,Exp[w],L2]/(1 + Q2 Exp[I fi] Exp[-w]),{w,0,100.}]+(1/Pi)NIntegrate[tr4l[nq,nu,k,Exp[w],L2]/(1 + Q2 Exp[I fi] Exp[-w]),{w,100.,$wwmax}, MaxRecursion->12]; *)



(* ::Subsection:: *)
(* This gives A_{nu}(Q^2) in global FAPT*)

AFAPT1lglob[nu_, k_, Q2_, L2nf3_, fi_] := ((1/Pi) NIntegrate[tr1l[3, nu, k, sig, L2nf3]/(sig + Q2 Exp[I fi]), {sig, 0, (Q2thr[3] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr1l[4, nu, k, sig, L21l[4, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[3] /. NumDefanQCD1), (Q2thr[4] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr1l[5, nu, k, sig, L21l[5, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[4] /. NumDefanQCD1), (Q2thr[5] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr1l[6, nu, k, sig, L21l[6, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[5] /. NumDefanQCD1), $wmax}] + (1/Pi) NIntegrate[tr1l[6, nu, k, Exp[w], L21l[6, L2nf3]]/(1 + Exp[-w] Q2 Exp[I fi]), {w, Log[$wmax], 10.^3}]+ (1/Pi) NIntegrate[tr1l[6, nu, k, Exp[w], L21l[6, L2nf3]]/(1 + Exp[-w] IntegerPart[Q2 10^20] 10^-20 Exp[I  IntegerPart[fi 10^2] 10^-20]), {w, 10.^3, 10^2*$wmax}, MaxRecursion->12]);
     (* AFAPT1lglob[nu_, k_, Q2_, L2nf3_, fi_] := ((1/Pi) NIntegrate[tr1l[3, nu, k, sig, L2nf3]/(sig + Q2 Exp[I fi]), {sig, 0, (Q2thr[3] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr1l[4, nu, k, sig, L21l[4, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[3] /. NumDefanQCD1), (Q2thr[4] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr1l[5, nu, k, sig, L21l[5, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[4] /. NumDefanQCD1), (Q2thr[5] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr1l[6, nu, k, sig, L21l[6, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[5] /. NumDefanQCD1), $wmax}] + (1/Pi) NIntegrate[tr1l[6, nu, k, Exp[w], L21l[6, L2nf3]]/(1 + Exp[-w] Q2 Exp[I fi]), {w, Log[$wmax], 10.^3}]+ (1/Pi) NIntegrate[tr1l[6, nu, k, Exp[w], L21l[6, L2nf3]]/(1 + Exp[-w] Q2 Exp[I fi]), {w, 10.^3, 10^2*$wmax}, MaxRecursion->12]); *)


AFAPT2lglob[nu_, k_, Q2_, L2nf3_, fi_] := ((1/Pi) NIntegrate[tr2l[3, nu, k, sig, L2nf3]/(sig + Q2 Exp[I fi]), {sig, 0, (Q2thr[3] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr2l[4, nu, k, sig, L22l[4, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[3] /. NumDefanQCD1), (Q2thr[4] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr2l[5, nu, k, sig, L22l[5, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[4] /. NumDefanQCD1), (Q2thr[5] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr2l[6, nu, k, sig, L22l[6, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[5] /. NumDefanQCD1), $wmax}] + (1/Pi) NIntegrate[tr2l[6, nu, k, Exp[w], L22l[6, L2nf3]]/(1 + Exp[-w] Q2 Exp[I fi]), {w, Log[$wmax], 10.^3}]+ (1/Pi) NIntegrate[tr2l[6, nu, k, Exp[w], L22l[6, L2nf3]]/(1 + Exp[-w] IntegerPart[Q2 10^20] 10^-20 Exp[I  IntegerPart[fi 10^2] 10^-20]), {w, 10.^3, 10^2*$wmax}, MaxRecursion->12]);
     (* AFAPT2lglob[nu_, k_, Q2_, L2nf3_, fi_] := ((1/Pi) NIntegrate[tr2l[3, nu, k, sig, L2nf3]/(sig + Q2 Exp[I fi]), {sig, 0, (Q2thr[3] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr2l[4, nu, k, sig, L22l[4, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[3] /. NumDefanQCD1), (Q2thr[4] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr2l[5, nu, k, sig, L22l[5, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[4] /. NumDefanQCD1), (Q2thr[5] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr2l[6, nu, k, sig, L22l[6, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[5] /. NumDefanQCD1), $wmax}] + (1/Pi) NIntegrate[tr2l[6, nu, k, Exp[w], L22l[6, L2nf3]]/(1 + Exp[-w] Q2 Exp[I fi]), {w, Log[$wmax], 10.^3}]+ (1/Pi) NIntegrate[tr2l[6, nu, k, Exp[w], L22l[6, L2nf3]]/(1 + Exp[-w] Q2 Exp[I fi]), {w, 10.^3, 10^2*$wmax}, MaxRecursion->12]); *)

AFAPT3lglob[nu_, k_, Q2_, L2nf3_, fi_] := ((1/Pi) NIntegrate[tr3l[3, nu, k, sig, L2nf3]/(sig + Q2 Exp[I fi]), {sig, 0, (Q2thr[3] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr3l[4, nu, k, sig, L23l[4, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[3] /. NumDefanQCD1), (Q2thr[4] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr3l[5, nu, k, sig, L23l[5, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[4] /. NumDefanQCD1), (Q2thr[5] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr3l[6, nu, k, sig, L23l[6, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[5] /. NumDefanQCD1), $wmax}] + (1/Pi) NIntegrate[tr3l[6, nu, k, Exp[w], L23l[6, L2nf3]]/(1 + Exp[-w] Q2 Exp[I fi]), {w, Log[$wmax], 10.^3}]+ (1/Pi) NIntegrate[tr3l[6, nu, k, Exp[w], L23l[6, L2nf3]]/(1 + Exp[-w] IntegerPart[Q2 10^20] 10^-20 Exp[I  IntegerPart[fi 10^2] 10^-20]), {w, 10.^3, 10^2*$wmax}, MaxRecursion->12]);
     (* AFAPT3lglob[nu_, k_, Q2_, L2nf3_, fi_] := ((1/Pi) NIntegrate[tr3l[3, nu, k, sig, L2nf3]/(sig + Q2 Exp[I fi]), {sig, 0, (Q2thr[3] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr3l[4, nu, k, sig, L23l[4, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[3] /. NumDefanQCD1), (Q2thr[4] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr3l[5, nu, k, sig, L23l[5, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[4] /. NumDefanQCD1), (Q2thr[5] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr3l[6, nu, k, sig, L23l[6, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[5] /. NumDefanQCD1), $wmax}] + (1/Pi) NIntegrate[tr3l[6, nu, k, Exp[w], L23l[6, L2nf3]]/(1 + Exp[-w] Q2 Exp[I fi]), {w, Log[$wmax], 10.^3}]+ (1/Pi) NIntegrate[tr3l[6, nu, k, Exp[w], L23l[6, L2nf3]]/(1 + Exp[-w] Q2 Exp[I fi]), {w, 10.^3, 10^2*$wmax}, MaxRecursion->12]); *)

AFAPT4lglob[nu_, k_, Q2_, L2nf3_, fi_] := ((1/Pi) NIntegrate[tr4l[3, nu, k, sig, L2nf3]/(sig + Q2 Exp[I fi]), {sig, 0, (Q2thr[3] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr4l[4, nu, k, sig, L24l[4, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[3] /. NumDefanQCD1), (Q2thr[4] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr4l[5, nu, k, sig, L24l[5, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[4] /. NumDefanQCD1), (Q2thr[5] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr4l[6, nu, k, sig, L24l[6, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[5] /. NumDefanQCD1), $wmax}] + (1/Pi) NIntegrate[tr4l[6, nu, k, Exp[w], L24l[6, L2nf3]]/(1 + Exp[-w] Q2 Exp[I fi]), {w, Log[$wmax], 10.^3}]+ (1/Pi) NIntegrate[tr4l[6, nu, k, Exp[w], L24l[6, L2nf3]]/(1 + Exp[-w] IntegerPart[Q2 10^20] 10^-20 Exp[I  IntegerPart[fi 10^2] 10^-20]), {w, 10.^3, $wwmax}, MaxRecursion->12]);
     (* AFAPT4lglob[nu_, k_, Q2_, L2nf3_, fi_] := ((1/Pi) NIntegrate[tr4l[3, nu, k, sig, L2nf3]/(sig + Q2 Exp[I fi]), {sig, 0, (Q2thr[3] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr4l[4, nu, k, sig, L24l[4, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[3] /. NumDefanQCD1), (Q2thr[4] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr4l[5, nu, k, sig, L24l[5, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[4] /. NumDefanQCD1), (Q2thr[5] /. NumDefanQCD1)}] + (1/Pi) NIntegrate[tr4l[6, nu, k, sig, L24l[6, L2nf3]]/(sig + Q2 Exp[I fi]), {sig, (Q2thr[5] /. NumDefanQCD1), $wmax}] + (1/Pi) NIntegrate[tr4l[6, nu, k, Exp[w], L24l[6, L2nf3]]/(1 + Exp[-w] Q2 Exp[I fi]), {w, Log[$wmax], 10.^3}]+ (1/Pi) NIntegrate[tr4l[6, nu, k, Exp[w], L24l[6, L2nf3]]/(1 + Exp[-w] Q2 Exp[I fi]), {w, 10.^3, $wwmax}, MaxRecursion->12]); *)



(* ::Section:: *)
(* Non-integer powers analogous *)
(* ES[m,n]=S_m(n) = Sum_{k=1^n}(1/k^m) *)


(* ::Subsection:: *)
(*Coefficients k1[m],k2[m],k3[m],k4[m] *)
(*which allow us to relate {tilde A_m} with (A_ 1)^m, for real m*)


ES[m_,n_]:= ((-1)^(m-1)/(m-1)!)*( PolyGamma[m-1,n+1] - PolyGamma[m-1,1] );
(* Now functions Z[m,n]=Zm(n) = (1/n!) (d^m/dx^m) (Gamma[n+1+x]/Gamma[1+x])|_{x=0} : *)
Z1[n_]:=ES[1,n];
Z2[n_]:= ES[1,n]^2 - ES[2,n];
Z3[n_]:=ES[1,n]^3 - 3 ES[1,n] ES[2,n] + 2 ES[3,n];
Z4[n_]:=ES[1,n]^4 - 6 ES[1,n]^2 ES[2,n] + 3 ES[2,n]^2 + 8 ES[1,n] ES[3,n] - 6 ES[4,n];


(* Below we have: kNf[s] = k_N(s); tkNf[s]={tilde k}_N(s); for N=1,2,3,4; s > 0 real positive *)
k1f[s_,cc1_]:=cc1 s (-1 + HarmonicNumber[s]);

k2f[s_,cc1_,cc2_]:= ( ((s-1) s/2) cc2 + cc1^2( (s (s+1)/2)( (PolyGamma[s+1]-PolyGamma[1])^2 + (PolyGamma[1,s+1]-PolyGamma[1,1]) + 1) - s^2(PolyGamma[s+1]-PolyGamma[1])) );

k3f[s_,cc1_,cc2_,cc3_]:= ( (s(s-1)(s+4)/12) cc3 + (s(s+1)(s+2)/2) cc1 cc2 ( (s/(s+2)) ES[1,s+1] + (2/(s+1)) + (1/(s+2)) - (11/6) ) + (s(s+1)(s+2)/6) cc1^3 ( Z3[s+2] - 3 Z2[s+1] + 3 ES[1,s] - 1 ) );

k4f[s_,cc1_,cc2_,cc3_,cc4_]:=((s(s+1)(s+2)(s+3))/6) ( cc4 (1/12 - 2/((s+1)(s+2)(s+3))) + cc2^2 (5/12 - (3 s^2 + 6 s +1)/((s+1)(s+2)(s+3))) + (cc1 cc3/2) ( (1 - 6/((s+2)(s+3))) Z1[s+1] + 4/(s+1) - 1/(s+2) - 1/(s+3) - 13/6 ) + (cc1^2 cc2/2) ( (3(s+1)/(s+3)) Z2[s+2] + (12/(s+2)+6/(s+3)-11) Z1[s+1] - 6/(s+1) - 5/(s+2) - 11/(s+3) + 38/3 ) + (cc1^4/4) ( Z4[s+3] - 4 Z3[s+2] + 6 Z2[s+1] - 4 Z1[s] + 1 ) );

tk1f[s_,cc1_]:=-k1f[s,cc1];
tk2f[s_,cc1_,cc2_]:=k1f[s,cc1] k1f[s+1,cc1]-k2f[s,cc1,cc2];
tk3f[s_,cc1_,cc2_,cc3_]:=-k1f[s,cc1] k1f[s+1,cc1] k1f[s+2,cc1] + k1f[s,cc1] k2f[s+1,cc1,cc2] + k1f[s+2,cc1] k2f[s,cc1,cc2] - k3f[s,cc1,cc2,cc3];
tk4f[s_,cc1_,cc2_,cc3_,cc4_]:=k1f[s,cc1] k1f[s+1,cc1] k1f[s+2,cc1] k1f[s+3,cc1] - k1f[s,cc1] k1f[s+1,cc1] k2f[s+2,cc1,cc2] - k1f[s,cc1] k1f[s+3,cc1] k2f[s+1,cc1,cc2] - k1f[s+2,cc1] k1f[s+3,cc1] k2f[s,cc1,cc2] + k2f[s,cc1,cc2] k2f[s+2,cc1,cc2] + k1f[s,cc1] k3f[s+1,cc1,cc2,cc3] + k1f[s+3,cc1] k3f[s,cc1,cc2,cc3] - k4f[s,cc1,cc2,cc3,cc4];


(* ::Section:: *)
(* 2danQCD model: *)
(* Analytic running coupling in our model; Q^2=u LLamb^2 Exp[I fi]; u>0 *)

(* ::Subsection::Closed:: *)
(* The constants of our model for different flavors, we need to take the threshold into account *)
(* The running from Q=MZ down to Q=kappa2d mb *)
nff=5;
Clear[Bta];
Bta[z_]:=- $be0[nff] z^2 (1 + $c1[nff] z + $c2MS[nff] z^2 + $c3MS[nff] z^3);
xx5=2*Log[(MZ/.NumDefanQCD)/(Q5/.NumDefanQCD1)];
del=0.01;
(* We find aa that is a(Q5+0)=a(MSbar,kappa2d*mb+0;nff=5) *)
     sol1MS = NDSolve[{y'[x]==Bta[y[x]], y[xx5]==aaMZ/.NumDefanQCD1}, y, {x,-del,+del}]; (*'*)
aa=First[y[0.]/.sol1MS];
aaQ5plus=aa;

(* Now transition from  a(Q5+0) to  a(Q5-0)=a(kappa2d*mb-0; nf=5-1): *)
lhint2d=2*Log[kappa2d];
xxi3=1.2020569031595942;
cct2=11/72;
cct3=- 82043*xxi3/27648 + 564731/124416 - 2633*(nff-1)/31104;
kk1= -lhint2d/6;
kk2= ( lhint2d^2/36 - 19*lhint2d/24 + cct2);
kk3=-lhint2d^3/216 - 131*lhint2d^2/576 + lhint2d*(-6793 + 281*(nff-1))/1728 +cct3;
aanew=aa*(1 + kk1*aa + kk2*aa^2 + kk3*aa^3);  
aa=aanew;
aaQ5low=aa;

(* The running in the nf=4 regime from Q5=kappa2d*mb to Q4=kappa2d*mc: *)
nff=4;
Clear[Bta];
Bta[z_]:=- $be0[nff] z^2 (1 + $c1[nff] z + $c2MS[nff] z^2 + $c3MS[nff] z^3);
xx4=2*Log[(Q5/.NumDefanQCD1)/(Q4/.NumDefanQCD1)];
del=0.01; (*'*)
sol1MS = NDSolve[{y'[x]==Bta[y[x]], y[xx4]==aaQ5low}, y, {x,-del,+del}];
aa=First[y[0.]/.sol1MS];
(* aa is now a(Q4+0)=a(MSbar,kappa2d*mc+0;nf=4) *)
aaQ4plus=aa;

(* Now transition from  a(Q4+0) to  a(Q4-0)=a(kappa2d*mc-0; nf=5-1): *)
lhint2d=2*Log[kappa2d];
xxi3=1.2020569031595942;
cct2=11/72;
cct3=- 82043*xxi3/27648 + 564731/124416 - 2633*(nff-1)/31104;
kk1= -lhint2d/6;
kk2= ( lhint2d^2/36 - 19*lhint2d/24 + cct2);
kk3=-lhint2d^3/216 - 131*lhint2d^2/576 + lhint2d*(-6793 + 281*(nff-1))/1728 +cct3;
aanew=aa*(1 + kk1*aa + kk2*aa^2 + kk3*aa^3);  
aa=aanew;
(* aa is now a(Q4+0)=a(MSbar,kappa2d*mc+0;nf=4) *)
aaQ4low=aa;

(* ainMS[3]=a(MSbar,Q=kappa2d*mc-0, nf=3) *)
ainMS[3]=aaQ4low;
(* ainMS[4]=a(MSbar,Q=kappa2d*mc+0, nf=4) *)
ainMS[4]=aaQ4plus;

(* The (4-loop MSbar) RGE-running UPWARDS from Q=MZ to kappa2d*mt, and beyond: *)
nff=5;
Clear[Bta];
Bta[z_]:=- $be0[nff] z^2 (1 + $c1[nff] z + $c2MS[nff] z^2 + $c3MS[nff] z^3);
xx5=2*Log[(MZ/.NumDefanQCD)/(Q6/.NumDefanQCD1)];
     del=0.01; (*'*)
sol1MS = NDSolve[{y'[x]==Bta[y[x]], y[xx5]==aaMZ/.NumDefanQCD1}, y, {x,-del,+del}];
(* aa is now a(Q6-0)=a(MSbar,kappa2d*mt-0; nff=5) *)
aa=First[y[0.]/.sol1MS];
aaQ6low=aa;

(* Now (3-loop) transition from  a(Q6-0) to  a(Q6+0)=a(kappa2d*mt+0; nf=5+1): *)
lh2d=2*Log[kappa2d];
xi3=1.202057;
ct2=11/72;
(* ct3=- 82043*xi3/27648 + 564731/124416 - 2633*(nff)/31104; *)
kk1= -lh2d/6;
kk2= ( lh2d^2/36 - 19*lh2d/24 + ct2);
kk3=-lh2d^3/216 - 131*lh2d^2/576 + lh2d*(-6793 + 281*(nff))/1728 +ct3[nff];
xx1=-kk1;
xx2=-kk2+2*kk1^2;
xx3=-kk3+5*kk2*kk1-5*kk1^3;
aanew=aa*(1 + xx1*aa + xx2*aa^2 + xx3*aa^3);  
aa=aanew;
aaQ6plus=aa;

(* ainMS[5]=a(MSbar,Q=kappa2d*mt-0, nf=5) *)
ainMS[5]=aaQ6low;
(* ainMS[6]=a(MSbar,Q=kappa2d*mt+0, nf=6) *)
ainMS[6]=aaQ6plus; 

(* Now we calculate ain[3]=al(ourRSch;Q=kappa2d*mc-0;nf=3)/Pi, ain[4]=al(ourRSch;Q=kappa2d*mc+0;nf=4)/Pi, ain[5]=al(ourRSch;Q=kappa2d*mt-0;nf=5)/Pi, and ain[6]=al(ourRSch;Q=kappa2d*mt+0;nf=6)/Pi *)
(* 'ourRSch' means: the renormalization scheme of the 2danQCD; i.e., with c2=c22din, c3=c2^2/c1, etc. *)

(* ain0[nf]: this is the first (pert.) approx.to a(kappa2d*mc;nf;c22d[nf],c3,...) when nf=3,4; and the first (pert.) approx.to a(kappa2d*mt;nf;c22d[nf],c3,...) when nf=5,6;*) 
ain0=ainMS;
fr[x_,nq_]:=($betaMS4l[x,nq] + $be0[nq] x^2 (1 + $c1[nq] x))/(x^2 (1+$c1[nq] x) $betaMS4l[x,nq]); 
rhs1[nq_]:=1/ainMS[nq]+ $c1[nq] Log[$c1[nq] ainMS[nq]/(1 + $c1[nq] ainMS[nq])] + NIntegrate[fr[x,nq],{x,0,ainMS[nq]}];

(* fl1[x_]:=approx ($c22d[nq] + c3 x)/((1+$c1[nq] x) (1 + $c1[nq] x + $c22d[nq] x^2 + c3 x^3)); *)
  fl1[x_,nq_]:=($btP11[x,nq] + $be0[nq] x^2 (1 + $c1[nq] x))/(x^2 (1+$c1[nq] x) $btP11[x,nq]);
(* aux1[ain] + NIntegrate[fl1[xx], {xx, $eepss, ain}] - rhs1 = 0 has to be *)
  aux1[yyy_,nq_]:=1/yyy + $c1[nq] Log[$c1[nq] yyy/(1 + $c1[nq] yyy)];

Do[Module[{}, LHS[yyy_? NumberQ]:=aux1[yyy,nqu] + NIntegrate[fl1[xxx,nqu], {xxx, $eepss, yyy}]; sol1=FindRoot[LHS[yyy]-rhs1[nqu],{yyy,ain0[nqu]}]; ain[nqu]=yyy/.sol1],{nqu,3,6}];

(* We now obtain the Lambert-Lambda LLamb[3] and LLamb[4], from the pQCD expressions via the condition a2d[nq,u=|Q2|/L2,fi=0]=ain[nq] for nq=3, 4 and |Q2|=Q^2=(kappa2d*mc)^2 = mu2in^2: *)

  Do[Module[{}, ub0=(kappa2d*mc/.NumDefanQCD)^2/0.26^2; sol0=FindRoot[a2d[nqu,uub,0]==ain[nqu],{uub,ub0}]; ub[nqu]=uub/.sol0; LL2[nqu]=(kappa2d*mc/.NumDefanQCD)^2/ub[nqu]; LLamb[nqu]=Sqrt[LL2[nqu]]],{nqu,3,4}];
  Module[{}, ub0=(kappa2d*mt/.NumDefanQCD)^2/0.097^2; sol0=FindRoot[a2d[5,uub,0]==ain[5],{uub,ub0}]; ub[5]=uub/.sol0; LL2[5]=(kappa2d*mt/.NumDefanQCD)^2/ub[5]; LLamb[5]=Sqrt[LL2[5]]];
 Module[{}, ub0=(kappa2d*mt/.NumDefanQCD)^2/0.03^2; sol0=FindRoot[a2d[6,uub,0]==ain[6],{uub,ub0}]; ub[6]=uub/.sol0; LL2[6]=(kappa2d*mt/.NumDefanQCD)^2/ub[6]; LLamb[6]=Sqrt[LL2[6]]];

  (* Now the evaluation of the 'pQCD-onset-scale' M0[Nf]=Sqrt[s0[Nf] LL2[Nf]], (in GeV) , for Nf=3,4,5,6: *)
  Do[M0[nqu]=Sqrt[$s0[nqu]LL2[nqu]],{nqu,3,6}];

(* We solve the equations for A1 - a ~ (L2/Q^2)^5, considering s0 as input and we use notation: int_{-sL[nq]}^s0  ds s^j tr12d[nq,s] = Pj (j=0,1,2,3) *)
soll = Solve[f1s + f2s == P0 && f1s s1 + f2s s2 == P1 && f1s s1^2 + f2s s2^2 == P2 && f1s s1^3 + f2s s2^3 == P3, {f1s, f2s, s1, s2}];
sollS = FullSimplify[soll, Assumptions -> P0 > 0 && P1 > 0 && P2 > 0 && P3 > 0];
  
(* fjs = fj^2 (j=1,2) *)
f1so1=f1s /. sollS[[1]];
f2so1=f2s /. sollS[[1]];
s1o1=s1 /. sollS[[1]];
s2o1=s2 /. sollS[[1]];

f1so2=f1s /. sollS[[2]];
f2so2=f2s /. sollS[[2]];
s1o2=s1 /. sollS[[2]];
s2o2=s2 /. sollS[[2]];


 Do[Module[{},If[$c22d[nqu]<0,solzL[nqu]= FindRoot[-1 + $c22d[nqu]/$c1[nqu]^2==ProductLog[-1,zu[nqu,uuL,0]],{uuL,2^(nqu-2) sL[nqu]}]]; If[$c22d[nqu]<0,uL[nqu]=uuL/.solzL[nqu]]; If[$c22d[nqu]>=0,uL[nq]=0]],{nqu,3,6}];
  dsL[nq_]:=Min[Abs[(uL[nq]-sL[nq])]/10.,10.^-6];

(* The numerical precision of calculation of Pjs0 is ~10^-12; *)
  Do[P0s0[nq]=(1/Pi) NIntegrate[tr12d[nq,s],{s,-sL[nq]+dsL[nq],0,$s0[nq]}]+ (dsL[nq]/2/Pi) NIntegrate[a2d[nq,Sqrt[sL[nq]^2 + (dsL[nq])^2 + 2 sL[nq] dsL[nq] Cos[Phi]],Arg[sL[nq]+dsL[nq] Exp[I Phi]]] Exp[I Phi],{Phi,-Pi,Pi}]+  HeavisideTheta[-$c22d[nq]-$small] (dsL[nq]/2/Pi) NIntegrate[a2d[nq,Sqrt[uL[nq]^2 + (dsL[nq])^2 + 2 uL[nq] dsL[nq] Cos[Phi]],Arg[uL[nq]+dsL[nq] Exp[I Phi]]] Exp[I Phi],{Phi,-Pi,Pi}],{nq,3,6}];
  Do[P1s0[nq]=(1/Pi)NIntegrate[s tr12d[nq,s],{s,-sL[nq]+dsL[nq],0,$s0[nq]}]+ (dsL[nq]/2/Pi) NIntegrate[a2d[nq,Sqrt[sL[nq]^2 + (dsL[nq])^2 + 2 sL[nq] dsL[nq] Cos[Phi]],Arg[sL[nq]+dsL[nq] Exp[I Phi]]] Exp[I Phi] (-sL[nq]-dsL[nq] Exp[I Phi]),{Phi,-Pi,Pi}]+  HeavisideTheta[-$c22d[nq]-$small] (dsL[nq]/2/Pi) NIntegrate[a2d[nq,Sqrt[uL[nq]^2 + (dsL[nq])^2 + 2 uL[nq] dsL[nq] Cos[Phi]],Arg[uL[nq]+dsL[nq] Exp[I Phi]]] Exp[I Phi] (-uL[nq]-dsL[nq] Exp[I Phi]),{Phi,-Pi,Pi}],{nq,3,6}]; 
  Do[P2s0[nq]=(1/Pi)NIntegrate[s^2 tr12d[nq,s],{s,-sL[nq]+dsL[nq],0,$s0[nq]}]+ (dsL[nq]/2/Pi) NIntegrate[a2d[nq,Sqrt[sL[nq]^2 + (dsL[nq])^2 + 2 sL[nq] dsL[nq] Cos[Phi]],Arg[sL[nq]+dsL[nq] Exp[I Phi]]] Exp[I Phi] (-sL[nq]-dsL[nq] Exp[I Phi])^2,{Phi,-Pi,Pi}]  +  HeavisideTheta[-$c22d[nq]-$small] (dsL[nq]/2/Pi) NIntegrate[a2d[nq,Sqrt[uL[nq]^2 + (dsL[nq])^2 + 2 uL[nq] dsL[nq] Cos[Phi]],Arg[uL[nq]+dsL[nq] Exp[I Phi]]] Exp[I Phi] (-uL[nq]-dsL[nq] Exp[I Phi])^2,{Phi,-Pi,Pi}],{nq,3,6}]; 
  Do[P3s0[nq]=(1/Pi)NIntegrate[s^3 tr12d[nq,s],{s,-sL[nq]+dsL[nq],0,$s0[nq]}]+ (dsL[nq]/2/Pi) NIntegrate[a2d[nq,Sqrt[sL[nq]^2 + (dsL[nq])^2 + 2 sL[nq] dsL[nq] Cos[Phi]],Arg[sL[nq]+dsL[nq] Exp[I Phi]]] Exp[I Phi] (-sL[nq]-dsL[nq] Exp[I Phi])^3,{Phi,-Pi,Pi}]+  HeavisideTheta[-$c22d[nq]-$small] (dsL[nq]/2/Pi) NIntegrate[a2d[nq,Sqrt[uL[nq]^2 + (dsL[nq])^2 + 2 uL[nq] dsL[nq] Cos[Phi]],Arg[uL[nq]+dsL[nq] Exp[I Phi]]] Exp[I Phi] (-uL[nq]-dsL[nq] Exp[I Phi])^3,{Phi,-Pi,Pi}],{nq,3,6}]; 

Do[Module[{},P0s0[nqu]=Re[P0s0[nqu]]; P1s0[nqu]=Re[P1s0[nqu]]; P2s0[nqu]=Re[P2s0[nqu]]; P3s0[nqu]=Re[P3s0[nqu]]],{nqu,3,6}];
Do[Module[{},f1ss0[nqu]= f1so1/.{P0->P0s0[nqu],P1->P1s0[nqu],P2->P2s0[nqu],P3->P3s0[nqu]}; f2ss0[nqu]= f2so1/.{P0->P0s0[nqu],P1->P1s0[nqu],P2->P2s0[nqu],P3->P3s0[nqu]}; s1s0[nqu]=s1o1/.{P0->P0s0[nqu],P1->P1s0[nqu],P2->P2s0[nqu],P3->P3s0[nqu]}; s2s0[nqu]=s2o1/.{P0->P0s0[nqu],P1->P1s0[nqu],P2->P2s0[nqu],P3->P3s0[nqu]}],{nqu,3,6}];
Do[lns0[nqu]=Log[$s0[nqu]],{nqu,3,6}];

(* Now calculation of the scales M1[Nf]=Sqrt[s1s0[Nf]LL2[Nf]] and M2[Nf]=Sqrt[s2s0[Nf]LL2[Nf]] (in GeV) where the two delta functions delta(sigma-Mj[Nf]) appear: *)
Do[M1[nqu]=Sqrt[s1s0[nqu]LL2[nqu]],{nqu,3,6}];
Do[M2[nqu]=Sqrt[s2s0[nqu]LL2[nqu]],{nqu,3,6}];


(* ::Subsection::Closed:: *)
(* First, an (efficient) calculation of A1(Q^2;nq) =A1an2d[nq,Q2,fi]: *)

Do[Module[{},J0[nq]=(1/Pi) NIntegrate[tr12d[nq,Exp[w]],{w,$lnsm,$wmaxJ}]; J1[nq]=(1/Pi) NIntegrate[tr12d[nq,Exp[w]] Exp[-w],{w,$lnsm,$wmax2d}]; J2[nq]=(1/Pi) NIntegrate[tr12d[nq,Exp[w]] Exp[-2 w],{w,$lnsm,$wmax2d}]; J3[nq]=(1/Pi) NIntegrate[tr12d[nq,Exp[w]] Exp[-3 w],{w,$lnsm,$wmax2d}]],{nq,3,6}];

(* A1[Q^2=Q2 Exp[I fi]; Nf)=A1an2d[Nf,Q2,fi]; note: Q2=|Q^2|; fi=Arg[Q^2] *)
A1an02d[nq_,Q2_,fi_]:=(1/Pi) NIntegrate[tr12d[nq,Exp[w]]/(1 + (Q2/LL2[nq]) Exp[I fi] Exp[-w]),{w,lns0[nq],$wmaxJ}] + f1ss0[nq]/( (Q2/LL2[nq]) Exp[I fi] +s1s0[nq]) + f2ss0[nq]/((Q2/LL2[nq]) Exp[I fi] +s2s0[nq]);
    (* A more efficient form of calculating A1: *)
A1an2d[nq_,Q2_,fi_]:= ( (1/Pi) NIntegrate[tr12d[nq,Exp[w]]/(1 + (Q2/LL2[nq]) Exp[I fi] Exp[-w]),{w,lns0[nq],$lnsm}] + f1ss0[nq]/( (Q2/LL2[nq]) Exp[I fi] +s1s0[nq])  + f2ss0[nq]/( (Q2/LL2[nq]) Exp[I fi] +s2s0[nq]) + J0[nq] - (Q2/LL2[nq]) Exp[I fi] J1[nq] + (Q2/LL2[nq])^2 Exp[I 2 fi] J2[nq] );

(* diffan(Q^2) = A1(Q^2) - a(Q^2), is very small for large Q^2 (> 1 GeV^2): *)
diffan[nq_,Q2_,fi_]:=A1an2d[nq,Q2,fi]-a2d[nq,(Q2/LL2[nq]),fi];
reldiffan[nq_,Q2_,fi_]:=A1an2d[nq,Q2,fi]/a2d[nq,Q2/LL2[nq],fi]-1;


(* ::Subsection::Closed:: *)
(* Calculation of tilde higher power analogs tA2d= {tilde A}_\nu(Q^2) for any real \nu in 2danQCD *)
 tr1MA1l[nq_,s_]:=(Pi/$be0[nq])(1/(Log[s]^2+Pi^2));
  
Delt[nu_]:=10.^-2 + HeavisideTheta[nu-3] (nu-3) 2 10.^-2; 
 Limt1l[u_,nu_,fi_]:=-(Zeta[1 - nu]/Gamma[nu]) + ((u Exp[I fi]-1) Zeta[-nu])/Gamma[nu] - ((u Exp[I fi]-1)^2 (Zeta[-1 - nu] + Zeta[-nu]))/(2 Gamma[nu]);

(* tAMA1l[nq,nu,u,fi]={tilde A}_{nu}^{MA,1-l.}(Q^2= u LL2[nq] Exp[I fi], nq) *)
  tAMA1l[nq_,nu_,u_,fi_]:=If[Abs[u Exp[I fi]-1]>Delt[nu],(1/$be0[nq]^nu)(1/(Log[u] + I fi)^nu - polylog[-nu+1,(1/u) Exp[-I fi]]/Gamma[nu]),(1/$be0[nq]^nu) Limt1l[u,nu,fi]];
 
(* tA2d[nq,nu,Q2,fi]= {tilde A}_{nu}(Q^2=Q2 Exp[I fi];nq); Q2=|Q^2| *)
(* This is Eq.(22) of JPG39,065005(2012) by G.Cvetic and A.V.Kotikov *)
(* tA2d[nq_,nu_,Q2_,fi_]:= tAMA1l[nq,nu,Q2/LL2[nq],fi] + ((-1)/(Pi*($be0[nq]^(nu-1))*Gamma[nu]))*((Pi f1ss0[nq]/s1s0[nq]) polylog[-nu+1,-s1s0[nq] Exp[-I fi] LL2[nq]/Q2] + (Pi f2ss0[nq]/s2s0[nq]) polylog[-nu+1,-s2s0[nq] Exp[-I fi] LL2[nq]/Q2] - NIntegrate[(tr1MA1l[nq,ssig]/ssig) polylog[-nu+1,-ssig Exp[-I fi] LL2[nq]/Q2],{ssig,0,$s0[nq]}] + NIntegrate[((tr12d[nq,ssig]-tr1MA1l[nq,ssig])/ssig) polylog[-nu+1,-ssig Exp[-I fi] LL2[nq]/Q2],{ssig,$s0[nq],10 $s0[nq]}] + NIntegrate[(tr12d[nq,Exp[w]]-tr1MA1l[nq,Exp[w]]) polylog[-nu+1,-Exp[w] Exp[-I fi] LL2[nq]/Q2],{w,Log[10 $s0[nq]],10^3}]);   *)

 tA2dpr[nq_,nupp_,Q2_,fi_]:= tAMA1l[nq,nupp,Q2/LL2[nq],fi] + ((-1)/(Pi*($be0[nq]^(nupp-1))*Gamma[nupp]))*((Pi f1ss0[nq]/s1s0[nq]) polylog[-nupp+1,-s1s0[nq] Exp[-I fi] LL2[nq]/Q2] + (Pi f2ss0[nq]/s2s0[nq]) polylog[-nupp+1,-s2s0[nq] Exp[-I fi] LL2[nq]/Q2] - NIntegrate[(tr1MA1l[nq,ssig]/ssig) polylog[-nupp+1,-ssig Exp[-I fi] LL2[nq]/Q2],{ssig,0,$s0[nq]}] + NIntegrate[((tr12d[nq,ssig]-tr1MA1l[nq,ssig])/ssig) polylog[-nupp+1,-ssig Exp[-I fi] LL2[nq]/Q2],{ssig,$s0[nq],10 $s0[nq]}] + NIntegrate[(tr12d[nq,Exp[w]]-tr1MA1l[nq,Exp[w]]) polylog[-nupp+1,-Exp[w] Exp[-I fi] LL2[nq]/Q2],{w,Log[10 $s0[nq]],$wmax2d}, MaxRecursion->12]);
 tA2d[nq_,nupp_,Q2_,fi_]:=If[fi==0,Re[tA2dpr[nq,nupp,Q2,0]],tA2dpr[nq,nupp,Q2,fi]];

(* ::Subsection::Closed:: *)
(* Calculation of higher power analogs A2dNl[Nf,M+nu,Q2,fi] = A_{M+nu}(Q^2; Nf) in N-loop approach (<-> N^{N-1}LO truncation approach); where Q^2 = Q2 Exp[i fi]; Q2=|Q^2|(= u LL2[Nf]) in GeV^2; \nu > -1 and real; M=0,1,...,(N-1) *)
A2d5l[nq_,4,nu_,Q2_,fi_]:=tA2d[nq,4+nu,Q2,fi];
A2d5l[nq_,3,nu_,Q2_,fi_]:=tA2d[nq,3+nu,Q2,fi] + tk1f[3+nu,$c1[nq]] tA2d[nq,4+nu,Q2,fi];
A2d5l[nq_,2,nu_,Q2_,fi_]:=tA2d[nq,2+nu,Q2,fi] + tk1f[2+nu,$c1[nq]] tA2d[nq,3+nu,Q2,fi] + tk2f[2+nu,$c1[nq],$c22d[nq]] tA2d[nq,4+nu,Q2,fi];
A2d5l[nq_,1,nu_,Q2_,fi_]:=tA2d[nq,1+nu,Q2,fi] + tk1f[1+nu,$c1[nq]] tA2d[nq,2+nu,Q2,fi] + tk2f[1+nu,$c1[nq],$c22d[nq]] tA2d[nq,3+nu,Q2,fi] + tk3f[1+nu,$c1[nq],$c22d[nq],$c32d[nq]] tA2d[nq,4+nu,Q2,fi];
A2d5l[nq_,0,nu_,Q2_,fi_]:=tA2d[nq,nu,Q2,fi] + tk1f[nu,$c1[nq]] tA2d[nq,1+nu,Q2,fi] + tk2f[nu,$c1[nq],$c22d[nq]]  tA2d[nq,2+nu,Q2,fi]+ tk3f[nu,$c1[nq],$c22d[nq],$c32d[nq]] tA2d[nq,3+nu,Q2,fi] + tk4f[nu,$c1[nq],$c22d[nq],$c32d[nq],$c42d[nq]] tA2d[nq,4+nu,Q2,fi];

A2d4l[nq_,3,nu_,Q2_,fi_]:=tA2d[nq,3+nu,Q2,fi];
A2d4l[nq_,2,nu_,Q2_,fi_]:=tA2d[nq,2+nu,Q2,fi] + tk1f[2+nu,$c1[nq]] tA2d[nq,3+nu,Q2,fi];
A2d4l[nq_,1,nu_,Q2_,fi_]:=tA2d[nq,1+nu,Q2,fi] + tk1f[1+nu,$c1[nq]] tA2d[nq,2+nu,Q2,fi] + tk2f[1+nu,$c1[nq],$c22d[nq]] tA2d[nq,3+nu,Q2,fi];
A2d4l[nq_,0,nu_,Q2_,fi_]:=tA2d[nq,nu,Q2,fi] + tk1f[nu,$c1[nq]] tA2d[nq,1+nu,Q2,fi] + tk2f[nu,$c1[nq],$c22d[nq]]  tA2d[nq,2+nu,Q2,fi]+ tk3f[nu,$c1[nq],$c22d[nq],$c32d[nq]] tA2d[nq,3+nu,Q2,fi];

A2d3l[nq_,2,nu_,Q2_,fi_]:=tA2d[nq,2+nu,Q2,fi];
A2d3l[nq_,1,nu_,Q2_,fi_]:=tA2d[nq,1+nu,Q2,fi] + tk1f[1+nu,$c1[nq]] tA2d[nq,2+nu,Q2,fi];
A2d3l[nq_,0,nu_,Q2_,fi_]:=tA2d[nq,nu,Q2,fi] + tk1f[nu,$c1[nq]] tA2d[nq,1+nu,Q2,fi] + tk2f[nu,$c1[nq],$c22d[nq]]  tA2d[nq,2+nu,Q2,fi];

A2d2l[nq_,1,nu_,Q2_,fi_]:=tA2d[nq,1+nu,Q2,fi];
A2d2l[nq_,0,nu_,Q2_,fi_]:=tA2d[nq,nu,Q2,fi] + tk1f[nu,$c1[nq]] tA2d[nq,1+nu,Q2,fi];

A2d1l[nq_,0,nu_,Q2_,fi_]:=tA2d[nq,nu,Q2,fi];


(* ::Section:: *)
(* Analytic massive running coupling MPT by Shirkov *)
(* Introduce xx=Log[(Q2+M2)/L2]=Log[uu] *)
(* And derivatives of massive coupling: (d/d Log[Q2])^N, for N=0,1,2,3: *)
     (* NOte: d/d Log[Q2] = (1 - all Exp[xx]) d/dxx, with xx=Log[(Q2/L2+all]=Log[(Q2+M2)/L2]: *)
(* NOte: ((-1)^N/(bet0^N N!)) ader[N,xx,nf] = {tilde A}_{N+1}(Q2=(Exp[xx]-all) bL2[nf],nf) *)
   (* ader[N,xx,nf] = (d/d Log[Q2])^N a(Q2+M2) *)

ader1l[0,xx_,M2_,L2_,nff_]:=1/$be0[nff]/xx;
ader1l[1,xx_,M2_,L2_,nff_]:=(1 - (M2/L2) Exp[-xxf]) D[1/$be0[nff]/xxf,xxf]/.{xxf->xx};
ader1l[2,xx_,M2_,L2_,nff_]:=(1 - (M2/L2) Exp[-xxf]) D[ader1l[1,xxf,M2,L2,nff],xxf]/.{xxf->xx};
ader1l[3,xx_,M2_,L2_,nff_]:=(1 - (M2/L2) Exp[-xxf]) D[ader1l[2,xxf,M2,L2,nff],xxf]/.{xxf->xx};
ader1l[4,xx_,M2_,L2_,nff_]:=(1 - (M2/L2) Exp[-xxf]) D[ader1l[3,xxf,M2,L2,nff],xxf]/.{xxf->xx};

ader2l[0,xx_,M2_,L2_,nff_]:=a2l[nff,Exp[xx]*L2/L2L2l[L2,nff],0];
ader2l[1,xx_,M2_,L2_,nff_]:=(1 - (M2/L2) Exp[-xxf]) D[a2l[nff,Exp[xxf]*L2/L2L2l[L2,nff],0],xxf]/.{xxf->xx};
ader2l[2,xx_,M2_,L2_,nff_]:=(1 - (M2/L2) Exp[-xxf]) D[ader2l[1,xxf,M2,L2,nff],xxf]/.{xxf->xx};
ader2l[3,xx_,M2_,L2_,nff_]:=(1 - (M2/L2) Exp[-xxf]) D[ader2l[2,xxf,M2,L2,nff],xxf]/.{xxf->xx};
ader2l[4,xx_,M2_,L2_,nff_]:=(1 - (M2/L2) Exp[-xxf]) D[ader2l[3,xxf,M2,L2,nff],xxf]/.{xxf->xx};

ader3l[0,xx_,M2_,L2_,nff_]:=a3l[nff,Exp[xx]*L2/L2L3l[L2,nff],0];
ader3l[1,xx_,M2_,L2_,nff_]:=(1 - (M2/L2) Exp[-xxf]) D[a3l[nff,Exp[xxf]*L2/L2L3l[L2,nff],0],xxf]/.{xxf->xx};
ader3l[2,xx_,M2_,L2_,nff_]:=(1 - (M2/L2) Exp[-xxf]) D[ader3l[1,xxf,M2,L2,nff],xxf]/.{xxf->xx};
ader3l[3,xx_,M2_,L2_,nff_]:=(1 - (M2/L2) Exp[-xxf]) D[ader3l[2,xxf,M2,L2,nff],xxf]/.{xxf->xx};
ader3l[4,xx_,M2_,L2_,nff_]:=(1 - (M2/L2) Exp[-xxf]) D[ader3l[3,xxf,M2,L2,nff],xxf]/.{xxf->xx};

ader4l[0,xx_,M2_,L2_,nff_]:=a4l[nff,Exp[xx],0];
ader4l[1,xx_,M2_,L2_,nff_]:=(1 - (M2/L2) Exp[-xxf]) D[a4l[nff,Exp[xxf],0],xxf]/.{xxf->xx};
ader4l[2,xx_,M2_,L2_,nff_]:=(1 - (M2/L2) Exp[-xxf]) D[ader4l[1,xxf,M2,L2,nff],xxf]/.{xxf->xx};
ader4l[3,xx_,M2_,L2_,nff_]:=(1 - (M2/L2) Exp[-xxf]) D[ader4l[2,xxf,M2,L2,nff],xxf]/.{xxf->xx};
ader4l[4,xx_,M2_,L2_,nff_]:=(1 - (M2/L2) Exp[-xxf]) D[ader4l[3,xxf,M2,L2,nff],xxf]/.{xxf->xx};


(* ::Subsection::Closed:: *)
(* For the use now, and for later cross-checks, we construct first: {tilde A}_N(Q2,nf)=tAAint[N,Q2,nf] for N=1,2,3,4,5 (integer index): *)
   (* Note: xx=Log[(Q2/L2[nf]+all[nf])]; Q2=L2[nf](Exp[xx]-all[nf]) *)
tAAint1l[1,Q2_,M2_,L2_,nff_]:= a1l[nff,Q2+M2,L2,0];
tAAint1l[2,Q2_,M2_,L2_,nff_]:=  (-1/$be0[nff]) QQ2 D[tAAint1l[1,QQ2,M2,L2,nff],QQ2]/.{QQ2->Q2};
tAAint1l[3,Q2_,M2_,L2_,nff_]:=  (-1/$be0[nff]/2) QQ2 D[tAAint1l[2,QQ2,M2,L2,nff],QQ2]/.{QQ2->Q2};
tAAint1l[4,Q2_,M2_,L2_,nff_]:=  (-1/$be0[nff]/3) QQ2 D[tAAint1l[3,QQ2,M2,L2,nff],QQ2]/.{QQ2->Q2};
tAAint1l[5,Q2_,M2_,L2_,nff_]:=  (-1/$be0[nff]/4) QQ2 D[tAAint1l[4,QQ2,M2,L2,nff],QQ2]/.{QQ2->Q2};

tAAint2l[1,Q2_,M2_,L2_,nff_]:= a2l[nff,Q2+M2,L2,0];
tAAint2l[2,Q2_,M2_,L2_,nff_]:=  (-1/$be0[nff]) QQ2 D[tAAint2l[1,QQ2,M2,L2,nff],QQ2]/.{QQ2->Q2};
tAAint2l[3,Q2_,M2_,L2_,nff_]:=  (-1/$be0[nff]/2) QQ2 D[tAAint2l[2,QQ2,M2,L2,nff],QQ2]/.{QQ2->Q2};
tAAint2l[4,Q2_,M2_,L2_,nff_]:=  (-1/$be0[nff]/3) QQ2 D[tAAint2l[3,QQ2,M2,L2,nff],QQ2]/.{QQ2->Q2};
tAAint2l[5,Q2_,M2_,L2_,nff_]:=  (-1/$be0[nff]/4) QQ2 D[tAAint2l[4,QQ2,M2,L2,nff],QQ2]/.{QQ2->Q2};

tAAint3l[1,Q2_,M2_,L2_,nff_]:= a3l[nff,Q2+M2,L2,0];
tAAint3l[2,Q2_,M2_,L2_,nff_]:=  (-1/$be0[nff]) QQ2 D[tAAint3l[1,QQ2,M2,L2,nff],QQ2]/.{QQ2->Q2};
tAAint3l[3,Q2_,M2_,L2_,nff_]:=  (-1/$be0[nff]/2) QQ2 D[tAAint3l[2,QQ2,M2,L2,nff],QQ2]/.{QQ2->Q2};
tAAint3l[4,Q2_,M2_,L2_,nff_]:=  (-1/$be0[nff]/3) QQ2 D[tAAint3l[3,QQ2,M2,L2,nff],QQ2]/.{QQ2->Q2};
tAAint3l[5,Q2_,M2_,L2_,nff_]:=  (-1/$be0[nff]/4) QQ2 D[tAAint3l[4,QQ2,M2,L2,nff],QQ2]/.{QQ2->Q2};

tAAint4l[1,Q2_,M2_,L2_,nff_]:= a4l[nff,Q2+M2,L2,0];
tAAint4l[2,Q2_,M2_,L2_,nff_]:=  (-1/$be0[nff]) QQ2 D[tAAint4l[1,QQ2,M2,L2,nff],QQ2]/.{QQ2->Q2};
tAAint4l[3,Q2_,M2_,L2_,nff_]:=  (-1/$be0[nff]/2) QQ2 D[tAAint4l[2,QQ2,M2,L2,nff],QQ2]/.{QQ2->Q2};
tAAint4l[4,Q2_,M2_,L2_,nff_]:=  (-1/$be0[nff]/3) QQ2 D[tAAint4l[3,QQ2,M2,L2,nff],QQ2]/.{QQ2->Q2};
tAAint4l[5,Q2_,M2_,L2_,nff_]:=  (-1/$be0[nff]/4) QQ2 D[tAAint4l[4,QQ2,M2,L2,nff],QQ2]/.{QQ2->Q2};



      (* Now we use Eq.(15) of GCAK (JPG 2012), for nu+1 -> nu=del+N, N=0,1,2,3,4; 0<del<1 *)
(* To get {tilde A}_{nu=N+del}(Q2,nf)=tAA[N,del,Q2,nf] in mQCD: *)
   (* NOte: xx=Log[Q2/bL2[3]+all[3]] = Log[(Q2+M2)/L2[3]]; Q2=Q^2 in GeV^2: *)
   (*--------------*)
      (* Use ziMIN=0.04 for N=0, ziMIN=0.04*(N+1) for N=1,2,3,4; and ziMAX=1-10.^-4, in tAAlong *)

 tAA1l[N_,del_,Q2_,M2_,L2_,nff_]:=(1/$be0[nff]^(N-1+del))*(Gamma[1+del]/Gamma[N+del])*(Sin[Pi del]/(Pi del))*( (-1)^N NIntegrate[(Exp[1/zi]/(E*zi^2))*(1/(Exp[1/zi]/E - 1)^del)*ader1l[N,Log[(Q2/L2/E) Exp[Exp[1/zi-1]]+M2/L2],M2,L2,nff],{zi,0.25,1.-10.^-4}, AccuracyGoal -> 8] + (-1)^N NIntegrate[(Exp[1/zi]/(E*zi^2))*(1/(Exp[1/zi]/E - 1)^del)*ader1l[N,Log[(Q2/L2/E)]+Exp[1/zi-1],0,L2,nff],{zi,0.04*(N+1),0.25}, AccuracyGoal -> 8] + $be0[nff]^N N! ((10.^-4)^(1-del)/(1-del)) tAAint1l[N+1,Q2,M2,L2,nff]*(1 + (3/2)(10^-4) (1 -del)) + (N!/$be0[nff])(Exp[N+del]/(N+del))Exp[-(N+del)/(N+1)/(0.04)]   );

tAA2l[N_,del_,Q2_,M2_,L2_,nff_]:=(1/$be0[nff]^(N-1+del))*(Gamma[1+del]/Gamma[N+del])*(Sin[Pi del]/(Pi del))*( (-1)^N NIntegrate[(Exp[1/zi]/(E*zi^2))*(1/(Exp[1/zi]/E - 1)^del)*ader2l[N,Log[(Q2/L2/E) Exp[Exp[1/zi-1]]+M2/L2],M2,L2,nff],{zi,0.25,1.-10.^-4}, AccuracyGoal -> 8] + (-1)^N NIntegrate[(Exp[1/zi]/(E*zi^2))*(1/(Exp[1/zi]/E - 1)^del)*ader2l[N,Log[(Q2/L2/E)]+Exp[1/zi-1],0,L2,nff],{zi,0.04*(N+1),0.25}, AccuracyGoal -> 8] + $be0[nff]^N N! ((10.^-4)^(1-del)/(1-del)) tAAint2l[N+1,Q2,M2,L2,nff]*(1 + (3/2)(10^-4) (1 -del)) + (N!/$be0[nff])(Exp[N+del]/(N+del))Exp[-(N+del)/(N+1)/(0.04)]   );

tAA3l[N_,del_,Q2_,M2_,L2_,nff_]:=(1/$be0[nff]^(N-1+del))*(Gamma[1+del]/Gamma[N+del])*(Sin[Pi del]/(Pi del))*( (-1)^N NIntegrate[(Exp[1/zi]/(E*zi^2))*(1/(Exp[1/zi]/E - 1)^del)*ader3l[N,Log[(Q2/L2/E) Exp[Exp[1/zi-1]]+M2/L2],M2,L2,nff],{zi,0.25,1.-10.^-4}, AccuracyGoal -> 8] + (-1)^N NIntegrate[(Exp[1/zi]/(E*zi^2))*(1/(Exp[1/zi]/E - 1)^del)*ader3l[N,Log[(Q2/L2/E)]+Exp[1/zi-1],0,L2,nff],{zi,0.04*(N+1),0.25}, AccuracyGoal -> 8] + $be0[nff]^N N! ((10.^-4)^(1-del)/(1-del)) tAAint3l[N+1,Q2,M2,L2,nff]*(1 + (3/2)(10^-4) (1 -del)) + (N!/$be0[nff])(Exp[N+del]/(N+del))Exp[-(N+del)/(N+1)/(0.04)]   );

tAA4l[N_,del_,Q2_,M2_,L2_,nff_]:=(1/$be0[nff]^(N-1+del))*(Gamma[1+del]/Gamma[N+del])*(Sin[Pi del]/(Pi del))*( (-1)^N NIntegrate[(Exp[1/zi]/(E*zi^2))*(1/(Exp[1/zi]/E - 1)^del)*ader4l[N,Log[(Q2/L2/E) Exp[Exp[1/zi-1]]+M2/L2],M2,L2,nff],{zi,0.25,1.-10.^-4}, AccuracyGoal -> 8] + (-1)^N NIntegrate[(Exp[1/zi]/(E*zi^2))*(1/(Exp[1/zi]/E - 1)^del)*ader4l[N,Log[(Q2/L2/E)]+Exp[1/zi-1],0,L2,nff],{zi,0.04*(N+1),0.25}, AccuracyGoal -> 8] + $be0[nff]^N N! ((10.^-4)^(1-del)/(1-del)) tAAint4l[N+1,Q2,M2,L2,nff]*(1 + (3/2)(10^-4) (1 -del)) + (N!/$be0[nff])(Exp[N+del]/(N+del))Exp[-(N+del)/(N+1)/(0.04)]   );

      (*--------------------------------_*)
tAMPT1l[Nq_,n_,del_,Q2_,M2_,L2_]:=tAA1l[n,del+$small,Q2,M2,L2,Nq];
tAMPT2l[Nq_,n_,del_,Q2_,M2_,L2_]:=tAA2l[n,del+$small,Q2,M2,L2,Nq];
tAMPT3l[Nq_,n_,del_,Q2_,M2_,L2_]:=tAA3l[n,del+$small,Q2,M2,L2,Nq];
tAMPT4l[Nq_,n_,del_,Q2_,M2_,L2_]:=tAA4l[n,del+$small,Q2,M2,L2,Nq];

      (* tAMPTXl[Nf,n,del,Q2,M2,L2]= {tilde A}_{nu=n+del}(Q^2=Q2) in MPT with m_{gl}^2=M2 (in GeV^2), and MSbar Lambda_{Nf}^2=L2, at X-loop (X=1,2,3,4); where Q2 is COMPLEX Q^2=|Q^2| Exp(I fi) (in GeV^2); nu = n+del, where n=0,1,3,4;  0=<del<1; *)

      (*--------------------------------_*)
(* ::Subsection::Closed:: *)
 (* Now construction of A_{nu}(Q2,nf) for general (noninteger) nu: *)
 
 (* AN[del,Q2,nf] = A_{del+N}(Q2), (nu = del+N), where Q2 is Q^2 in GeV, can be complex *)
 (* where 0=<del<1, N=1,2,3,4; nf=3,4,5,6 *)


AMPT1l[nff_,nu_,Q2_,M2_,L2_]:=   If[0<=nu<1,tAA1l[0,nu+$small,Q2,M2,L2,nff],If[1<=nu<2,tAA1l[1,nu-1+$small,Q2,M2,L2,nff],If[2<=nu<3,tAA1l[2,nu-2+$small,Q2,M2,L2,nff],If[3<=nu<4,tAA1l[3,nu-3+$small,Q2,M2,L2,nff],If[4<=nu<5,tAA1l[4,nu-4+$small,Q2,M2,L2,nff]]]]]];

 AMPT2l[nff_,nu_,Q2_,M2_,L2_]:=   
If[0<=nu<1,
tAA2l[0,nu+$small,Q2,M2,L2,nff]+ tk1f[nu+$small,cc1] tAA2l[1,nu+$small,Q2,M2,L2,nff],
If[1<=nu<2,
tAA2l[1,nu-1+$small,Q2,M2,L2,nff]+tk1f[nu+$small,cc1] tAA2l[2,nu-1+$small,Q2,M2,L2,nff],
If[2<=nu<3,
tAA2l[2,nu-2+$small,Q2,M2,L2,nff]+tk1f[nu+$small,cc1] tAA2l[3,nu-2+$small,Q2,M2,L2,nff],
If[3<=nu<4,
tAA2l[3,nu-3+$small,Q2,M2,L2,nff]+tk1f[nu+$small,cc1] tAA2l[4,nu-3+$small,Q2,M2,L2,nff],
If[4<=nu<5,
tAA2l[4,nu-4+$small,Q2,M2,L2,nff]
]]]]]/.{cc1->$c1[nff]};

 AMPT3l[nff_,nu_,Q2_,M2_,L2_]:=   
If[0<=nu<1,
tAA3l[0,nu+$small,Q2,M2,L2,nff]+ tk1f[nu+$small,cc1] tAA3l[1,nu+$small,Q2,M2,L2,nff]+ tk2f[nu+$small,cc1,cc2]  tAA3l[2,nu+$small,Q2,M2,L2,nff],
If[1<=nu<2,
tAA3l[1,nu-1+$small,Q2,M2,L2,nff]+tk1f[nu+$small,cc1] tAA3l[2,nu-1+$small,Q2,M2,L2,nff]+ tk2f[nu+$small,cc1,cc2]  tAA3l[3,nu-1+$small,Q2,M2,L2,nff],
If[2<=nu<3,
tAA3l[2,nu-2+$small,Q2,M2,L2,nff]+tk1f[nu+$small,cc1] tAA3l[3,nu-2+$small,Q2,M2,L2,nff]+ tk2f[nu+$small,cc1,cc2]  tAA3l[4,nu-2+$small,Q2,M2,L2,nff],
If[3<=nu<4,
tAA3l[3,nu-3+$small,Q2,M2,L2,nff]+tk1f[nu+$small,cc1] tAA3l[4,nu-3+$small,Q2,M2,L2,nff],
If[4<=nu<5,
tAA3l[4,nu-4+$small,Q2,M2,L2,nff]
   ]]]]]/.{cc1->$c1[nff],cc2->$c2MS[nff]};


AMPT4l[nff_,nu_,Q2_,M2_,L2_]:=   
If[0<=nu<1,
tAA4l[0,nu+$small,Q2,M2,L2,nff]+ tk1f[nu+$small,cc1] tAA4l[1,nu+$small,Q2,M2,L2,nff]+ tk2f[nu+$small,cc1,cc2]  tAA4l[2,nu+$small,Q2,M2,L2,nff]+tk3f[nu+$small,cc1,cc2,cc3] tAA4l[3,nu+$small,Q2,M2,L2,nff],
If[1<=nu<2,
tAA4l[1,nu-1+$small,Q2,M2,L2,nff]+tk1f[nu+$small,cc1] tAA4l[2,nu-1+$small,Q2,M2,L2,nff]+ tk2f[nu+$small,cc1,cc2]  tAA4l[3,nu-1+$small,Q2,M2,L2,nff]+tk3f[nu+$small,cc1,cc2,cc3] tAA4l[4,nu-1+$small,Q2,M2,L2,nff],
If[2<=nu<3,
tAA4l[2,nu-2+$small,Q2,M2,L2,nff]+tk1f[nu+$small,cc1] tAA4l[3,nu-2+$small,Q2,M2,L2,nff]+ tk2f[nu+$small,cc1,cc2]  tAA4l[4,nu-2+$small,Q2,M2,L2,nff],
If[3<=nu<4,
tAA4l[3,nu-3+$small,Q2,M2,L2,nff]+tk1f[nu+$small,cc1] tAA4l[4,nu-3+$small,Q2,M2,L2,nff],
If[4<=nu<5,
tAA4l[4,nu-4+$small,Q2,M2,L2,nff]
   ]]]]]/.{cc1->$c1[nff],cc2->$c2MS[nff],cc3->$c3MS[nff]};



EndPackage[];
