(* 2danQCDrtalMath11.m  This is the version which works for Mathematica 11 *)
(* It uses the default PolyLog, and not any more the polylog of Li__nu.m *)
(*--------------------------------------------------------------*)
(* Program which calculates the couplings of the 2delta anQCD model (2danQCD), in the Lambert scheme with a given (chosen) c2(=beta_2/beta_0) < 0 (so that there is uL pole of a2d pQCD coupling)  and for chosen s0[3] and alpha_s(MZ^2;MSbar); it calculates rtau(D=0) with dBG22, TPS and LB+bLB methods *)
(*----------------------------------------------------------------*)
(* Couplings, functions: *)
LL2::usage= "LL2[Nf]=Lambda_L^2, this is square the Lambert Lambda scale Lambda_L (with Nf flavors), in GeV^2, appearing in the pQCD running coupling a(Q^2) (see: a2d) underlying for the considered Nf-flavor 2danQCD";
alphaMSMZ::usage= "alphaMSMZ=alpha_s(Q^2=MZ^2,MSbar). In this case, alphaMSMZ=0.1185, or 0.1181, or 0.1189"
c22din::usage= "c22din=c2=beta_2/beta_0 (chosen to be c2ddin=-4.9) is the value of the c2 coefficient (for Nf=3,4,5,6) in the Lambert renormalization scheme used in the considered 2danQCD [c2=-4.9; c3=c2^2/c1. c4=c2^3/c1^2, etc., and cj=beta_j/beta_0; normalization convention: beta_0=(11 -2*Nf/3)/4; beta_1=(102-38*Nf/3)/16; etc.]"
s0::usage= "s0[Nf]=s_0 is a dimensionless parameter for Nf flavor 2danQCD (Nf=3,4,5,6); it is equal for Nf=3,4,5,6; the pQCD-onset scale M_0^2 in Nf=3 2danQCD is then: M_0^2(Nf=3)=s0[3]*LL2[3] (in GeV^2), and similarly for other cases of Nf"
s1s0::usage= "s1s0[Nf]=s_1 is a dimensionless parameter for Nf flavor 2danQCD (Nf=3,4,5,6); the location of the first delta function of the discontinuity function rho_1(sigma)=Im A(Q^2=-sigma - I epsilon) of the A-coupling of the Nf-flavor 2danQCD is then at sigma=M_1^2=s1s0[Nf]*LL2[Nf] (in GeV^2)"
s2s0::usage= "s2s0[Nf]=s_2 is a dimensionless parameter for Nf flavor 2danQCD (Nf=3,4,5,6); the location of the second delta function of the discontinuity function rho_1(sigma)=Im A(Q^2=-sigma - I epsilon) of the A-coupling of the Nf-flavor 2danQCD is then at sigma=M_2^2=s2s0[Nf]*LL2[Nf] (in GeV^2)"
f1ss0::usage= "f1ss0[Nf]=f_1^2 is a dimensionless strength parameter at the first delta in rho_1(sigma) discontinuity function of A(Q^2), in the Nf-flavor 2danQCD; the first delta function in rho_1(sigma) is Pi*f1ss0[Nf]*LL2[Nf]*delta(sigma - M1^2[Nf]), where M1^2[Nf]= s1s0[Mf]*LL2[Nf] (in GeV^2)"
f2ss0::usage= "f2ss0[Nf]=f_2^2 is a dimensionless strength parameter at the second delta in rho_1(sigma) discontinuity function of A(Q^2), in the Nf-flavor 2danQCD; the second delta function in rho_1(sigma) is Pi*f2ss0[Nf]*LL2[Nf]*delta(sigma - M2^2[Nf]), where M2^2[Nf]= s2s0[Mf]*LL2[Nf] (in GeV^2)"  
a2d::usage= "a2d[Nf,u,fi]= alpha_s(Q^2=u*LL2[Nf]*Exp[I fi)/Pi =a(Q^2), where Q^2=|Q^2|*Exp[I fi] and |Q^2|=u*LL2[Nf] (see: LL2); a2d is the underlying pQCD coupling in the Lambert renormalization scheme [c2=c22din=-4.9; c3=c2^2/c1. c4=c2^3/c1^2, etc., and cj=beta_j/beta_0; normalization conventions: beta_0=(11 -2*Nf/3)/4; beta_1=(102-38*Nf/3)/16; etc.], the pQCD coupling underlying the Nf-flavor 2danQCD coupling A (see: A1).";
tr12d::usage= "tr12d[Nf,s] is the discontinuity function rho_1^{(pt)}(sigma)=Im a(Q^2=-sigma - I epsilon) of the underlying pQCD coupling a(Q^2) (see: a2d) in the considered Nf-flavor 2danQCD, where s=sigma/LL2[Nf] (see: LL2)"
A1::usage= "A1[Nf,Q2,fi] = A(Q^2) = A_1(Q^2), where Q^2=Q2 Exp[I fi] (Q2=|Q^2|, in GeV^2) and A(Q^2) is the 2danQCD analog of the underlying perturbative coupling a(Q^2) = a2d[Nf,Q2/LL2[Nf],fi] (see: a2d); at high |Q^2|> 1 GeV^2 the two couplings practically coincide";
tA2::usage= "tA2[Nf,Q2,fi]={tilde A}_2(Q^2) for Q^2=Q2 Exp[I fi] (Q2=|Q^2|, in GeV^2); this is the first logarithmic derivative of A_1 (see: A1): {tilde A}_2(Q^2) = (-1)/(beta_0 1!) (d/d ln Q^2) A_1(Q^2)";
tA3::usage= "tA3[Nf,Q2,fi]={tilde A}_3(Q^2) for Q^2=Q2 Exp[I fi]; this is the second logarithmic derivative of A_1 (see: A1): {tilde A}_3(Q^2) = (-1)^2/(beta_0^2 2!) (d/d ln Q^2)^2 A_1(Q^2)";
tA4::usage= "tA4[Nf,Q2,fi]={tilde A}_4(Q^2) for Q^2=Q2 Exp[I fi]; this is the third logarithmic derivative of A_1 (see: A1): {tilde A}_4(Q^2) = (-1)^3/(beta_0^3 3!) (d/d ln Q^2)^3 A_1(Q^2)";
tA2d::usage= "tA2d[Nf,nu,Q2,fi]= {tilde A}_{nu}(Q^2), with Q^2=Q2 Exp[I fi]  (Q2=|Q^2|); {tilde A}_{nu} can have nu noninteger; {tilde A}_{nu} is the generalization (to noninteger nu) of the logarithmic derivative (-1)^{nu-1}/(beta_0^{nu-1} (nu-1)! (d/d ln Q^2)^(nu-1) A_1(Q^2)";
A2dNX::usage= "A2dNX[Nf,M,nu,Q2,fi] = A_{nu+M}(Q^2) is the 2danQCD analog of the power of the pQCD coupling a2d[Nf,Q2/LL2[Nf],fi]^(nu+M)=a(Q^2)^(nu+M) (see: a2d), for NX=N1,N2,N3,N4; here, -1 =< nu =< 1 (in general noninteger), and M is integer, M=0,1,..,X; where Q^2=Q2 Exp[I fi] and Q2=|Q^2|=u*LL2[nq] (in GeV^2). Note that for X=1 (N1), A2dN1[...] are the couplings (A_{nu}, A_{nu+1}) in the NLO sense, when in the TPS of a considered observable we only have two terms known (~a^nu and ~a^(nu+1)); for X=2 (N2), A2dN2[...] are the couplings (A_{nu}, A_{nu+1}, A_{nu+2}) in the N^2LO sense, i.e.,  when in the truncated perturbation series of a considered observable we have three terms known (~a^nu, ~a^(nu+1) and ~a^(nu+2)); etc.; Nf indicates the chosen number of active flavors of the 2danQCD (Nf=3,4,5,6), usually Nf=3."
AYNX::usage= "AYNX[Nf,Q2,fi] = A_{Y}(Q^2) is the 2danQCD analog of the power of the pQCD coupling a2d[Nf,Q2/LL2[Nf],fi]^Y=a(Q^2)^Y, for NX=N1,N2,N3; here, Y is integer, Y=1,..,X+1; where Q^2=Q2*Exp[I fi] and Q2=|Q^2|= u*LL2[nq] (in GeV^2). Note that X=1 (N1) is used when we have in the truncated perturbation series of a considered observable only two terms (~a and ~a^2) known; X=2 (N2) is used when three terms are known (~a, ~a^2 and ~a^3); X=3 (N3) is used when four terms are known (~a, ~a^2, ~a^3 and ~a^4). Nf indicates the chosen number of active flavors of the 2danQCD (Nf=3,4,5,6), usually Nf=3."
tdj::usage= "tdj={tilde d}_j (j=1,2,3) are the coefficients for the Nf=3 canonical (and massless) Adler function in the Lambert renormalization scheme used in the considered 2danQCD and in the underlying pQCD [c2=-4.9; c3=c2^2/c1. c4=c2^3/c1^2, etc., and cj=beta_j/beta_0; normalizations: beta_0=(11 -2 Nf/3)/4; beta_1=(102-38 Nf/3)/16; etc.]"
tdjbLB::usage= "tdjbLB={tilde d}_j - {tilde d}_j^(LB) are the 'beyond-the-leading-beta0' (bLB) part of {tilde d}_j Nf=3 coefficients of the canonical (and massless) Adler function in the Lambert renormalization scheme used in the considered 2danQCD and in the underlying pQCD (see also: tdj)."
rtauLBbLB::usage= "rtauLBbLB=r_{tau}(D=0) calculated (in the considered Lambert renormalization scheme) in the considered 2danQCD with the 'leading-beta_0 (LB) + beyond-the leading-beta_0 (bLB)' approach with Nf=3. Note that here rtauLBbLB=0.2010, or 0.2030, or 0.1990"
rtaudBG22::usage= "rtaudBG22=r_{tau}(D=0) calculated in the considered 2danQCD with the diagonal [2/2] Baker-Gammel approach."
rtauTPS::usage= "rtauTPS=r_{tau}(D=0) calculated in the considered 2danQCD with the truncated perturbation series approach."
aMSnf3::usage= "aMSnf3[Q2] is the pQCD coupling a(Q^2=Q2)=alpha_s(Q2)/Pi for positive Q2 (Q2 in GeV^2) in MSbar renormalization scheme, for the number of active flavors Nf=3; it corresponds to alpha_s(MZ^2;MSbar)=0.1185, or 0.1181, or 0.1189,  and is obtained by 4-loop RGE evolution in MSbar scheme from Q2=M_Z^2 downwards, with 3-loop quark threshold conditions imposed at threshold momenta Q2=(2*mqbar)^2, where mqbar=mqbar(Q2=mqbar^2) are the MSbar quark masses of b (for Nf=5->4) and of c (for Nf=4->3); we used mbbar=4.2 GeV and mcbar=1.27 GeV."
  (*---------------------------------------------------------------------------*)  

(* 2danQCDrt201al01185.m: as 2danQCDforMathematicaal01185_ext.m, but now no explicit calculation of rtau, to save time *)
(* This program gives, for nf=3: rtau(DS=0,mq=0,D=0)=0.201; it has as the reference value alphaMSMZ=0.1185; and uses the scheme parameter c2=-4.9, for Lambert scheme *)
(*-------*)
(* After running this program (by the command: << 2danQCDforMathematicaal01185_extSH.m;), we can evaluate the following couplings: *)
(* A1[nf,Q2,fi] = A(Q^2=Q2 Exp[I fi]; nf), where Q2=|Q^2|, and nf=number of active quarks (nf=3,4,5,6) *)
(* tAXan[nf,Q2,fi] = {tilde A}_{X}(Q^2=Q2 Exp[I fi]; nf), where Q2=|Q^2|, and nf=number of active quarks; X=2,3,4; *)
(* The general {tilde A}_{nu} coup0.1ling: *)
(* tA2d[nf,nu,Q2,fi]= {tilde A}_{nu}(Q^2=Q2 Exp[I fi];nf); Q2=|Q^2| and nf=number of active quarks *)
(* --------------------------*)
(* --------------------------*)
(* --------------------------*)
(* --------------------------*)
(* 2danQCDforMathematicaal01185_ext.m: as 2danQCDforMathematica_ext.m, but for alphaMSMZ=0.1185 (and not 0.1184), and for fixed (not variable) c22din=-4.9 *)
(* It is 'ext' in the sense that it generates 2danQCD model not just for nq=3, but also for nq=4,5,6; it imposes the condition rtau(DS=0,mq=0,D=0)=0.201 (NOT any more: 0.203) *)
(* --------------------------*)
(* 2danQCDforMathematica_ext.m: as 2danQCDforMathematica.m, but not calling subroutines RGErundown.m (+ RGE3matchdown.m) and RGErunup.m (+RGEmatch.m) -- these parts are now explicitly incorporated into the program; also notation is such that '2d' (meaning: 2-delta) is added to most of the relevant functions: *)
(* c2in -> c22din; cN[nq]->cN2d[nq] (N=2,3,4,5,6); zu[...]->zu2d[...]; *)
(* L2[nq] -> LL2[nq]  (=LambertLambda[nq]^2, in GeV^2) *)
(* kappa->kappa2d;  *)
(* a[...]->a2d[...]; tr1[...]->tr12d[...]; adis[...] -> adis2d[...]; *)
(* adis2[...]->adis2d2[...]; adis2simpl[...]->adis2d2simpl[...]; *)
(* tA[nq,nupp,u,fi] -> tA2d[nq,nupp,Q2,fi]; tAwm[nq,nupp,u,fi,wm]->tA2dwm[nq,nupp,Q2,fi,wm];  where  Q2=u LL2[nq] = |Q^2| in GeV^2 *)
(* AMNX[nq,nupp,u,fi] -> A2dNX[nq,M,nupp,Q2,fi] for NX=N2,N3,N4; and M=0,1,..,X; where Q2=u LL2[nq] = |Q^2| in GeV^2 *)
(* Where: A2dNX[nq,M,nu,Q2,fi]=A_{M+nu}(Q^2=Q2 Exp[I fi};Nf=nq)^{(NXLO)} in "NXLO" order (Q2=|Q^2|; fi=Arg[Q^2] where Q^2 in GeV^2), i.e., "NXLO" means that these couplings are used for the evaluation of truncated sums of the type: *)
(* --------------------------*)
(* D[Q2]_{pt}^{NXLO} = a^{nu} + d1 a^{nu+1} + d2 a^{nu+2} + ... + dX a^{nu+X} [where: a=a(mu^2)] *)
(* Note: D[Q2]_{mpt}^{NXLO} = {tilde a}_{nu} + {tilde d1} {tilde a}_{nu+1} + {tilde d2} {tilde a}_{nu+2} + ... {tilde dX} {tilde a}_{nu+X}; *)
(* {tilde d1}= {tilde k}_1(nu) + d1; *)
(* {tilde d2}={tilde k}_2(nu)+ d1 {tilde k}_1(nu+1) + d2;  ... ; *)
(* {tilde dX}={tilde k}_X(nu) + Sum_{NN=1}^{X-1} dNN {tilde k}_{X-NN}(nu+NN) + dX; *)
(* Then:  D[Q2]_{mpt}^{NXLO}-D[Q2]_{pt} = O(a^{nu+X+1}) *)

(*  D[Q2]_{an}^{NXLO} = A_{nu}^{(NXLO)} + d1 A_{nu+1}^{(NXLO)} + d2 A_{nu+2}^{(NXLO)} + ... + dX A_{nu+X}^{(NXLO)} *)
(*  D[Q2]_{man}^{NXLO} = tA_{nu} + {tilde d1} tA_{nu+1} + {tilde d2} tA_{nu+2} + ... + {tilde dX} tA_{nu+X} *)
(* where A_{nu+m}^{(NXLO)}=A_{nu+m}(mu^2)^{(NXLO)}; tA_{nu+m)=tA_{nu+m}(mu^2) *)

(* Then it is straightforward to see that: D[Q2]_{an}^{NXLO} = D[Q2]_{man}^{NXLO} *)
(*-------------------------------------------*)
(* 2danQCDforMathematica.m: as 2danQCDnf3genRSch, but now not just for a general c2-Lambert renormalization scheme (-6 < c2 < -2, cj=c2^{j-1}/c1^{j-2} (j=3,4,...)), but also for ALL 4 CHOICES of Nf: Nf=3,4,5,6 *)
(* Namely, for Nf->Nf+1, s0[Nf+1] is fixed by demanding A1[kappa2d*mq+0,Nf+1] - A1[kappa2d*mb-0,Nf] = a2d[kappa2d*mq+0,Nf+1]- a2d[kappa2d*mq-0,Nf], where the latter difference of corresponding underlying pQCD couplings a2d[Nf=4] and a2d[Nf=3] is obtained by using the CKS 3-loop matching in MSbar [CKS, PLB79,2184(1997) Eqs.(6)-(7)] translated to the c2-Lambert RSch with the same scales kappa2d2d*mq  (kappa2d=2 used; mq=mqbar(mqbar)); for Nf=3 -> Nf+1=4 we have mq=mc=1.27 (in GeV); for Nf=4 ->Nf+1=5 we have mb=4.20; for Nf=5 -> Nf+1=6 we have mt=163; note, these are MSbar masses, not pole masses *)
(*-------------------------------------------*)
(* 2danQCDnf3genRSch.m: as 2delanQCDAnu_downw_ext.m, but only for nf=3 and for general RSch parameter c2: -2 < c2 < -6; i.e., s0 becomes a function of c2 *)
(*-------------------------------------------*)
(* It calculates, in 2delanQCDA, the values of {tilde A}_{nu}(Q2) = tA2d[nu,u=|Q2|/L2, fi], where L2=LambdaLambert^2, and Q2=Q^2=|Q2| Exp[I fi]; u=|Q2|/L2 *)
(* And it calculates the (4-loop truncated) values of A_{N+nu}(Q2)=AN[nu,u=|Q2|/L2,fi] for N=0,1,2,3,4 and 0<nu<1 *)
(*----------------------------------------------------------------------*)
(* Our 2delanQCD model, which uses Gardi et al. beta function [use of Lambert functions] in c2-Lambert Rsch *)
(*-------------------------------------------------------------------*)
(* Search for strings 'Input' to put in the values of the Input parameters *)
(*-------------------------------------------------------------------*)
(*  Input c2=c22din (to define our c2-Lambert renormalization scheme; cj=c2^{j-1}/c1^{j-2}, j=3,4,5,... ;  note that c2in has the same value for all Nf=3,4,5,6, in our approach; we have the allowed integral of variation of c2in: -6 < c2in < -2 (so that A1[0]< 1.00 and M0< 1.5 GeV): *)
 (* c22din=-4.9; *)
 c22din=-4.9;
 (*---------------------------_*)
Needs["NumericalCalculus`"]


 (* Now the corresponding s0[Nf]=M0^2/LL2[Nf], where M0 is the pQCD-onset scale in the discontinuity function; and L2 is LambertLambda^2; Nf=3 at first; also extended to Nf=4,5,6 in order to simulate at high energies the pQCD, but now more suitable for the application of the logarithmic derivatives {tilde A}_{nu}(Q^2; Nf) for nu general (not necessarily integer) *)

 (* Input: s0[3], alphaMSMZ: *)
 (* either s0[3]=28.3; alphaMSMZ=0.1189; *)
 (* or s0[3]=25.6; alphaMSMZ=0.1185; *)
 (* or s0[3]=22.7; alphaMSMZ=0.1181; *)
 (* in all these cases we get: rtau(LB+bLB)=0.201 *)
 (* s0[3]=23.0757; *)
 s0[3]=22.7;
 (* s0[3]=s0[Nf=3], will be adjusted such that such that rtau=0.201 (=rtau(LB+bLB)), in NF=3 2danQCD; note that we take usually alphaMSMZ=0.1185 *)
 alphaMSMZ=0.1181;   (* alphaMSMZ=alpha_s(MZ^2;MSbar) *)
 (*-----------------------------------------------_*)
 (* Input: s0[4], s0[5], s0[6]: *)
 (* s0[4]=2*s0[3]; s0[5]=4*s0[3];  s0[6]=20*s0[3]; *)
s0[4]=s0[3]; 
 (* s0[5]=8*s0[3];  s0[6]=30*s0[3];  *)
 s0[5]=s0[3]; s0[6]=s0[3];

 (* Note: for s0[5]=s0[3] we get s0s2[5]<0; s0[6]=s0[3] -> s0s2[6]<0 (very negative) -- should increase very mucg s0[5], and even much more inclrease s0[6] *) 

 (* If s0[5]=80.; s0[6]=500; --> M0[5]=Sqrt[s0[5]] LLamb[5]=0.857593; M0[6]=Sqrt[s0[6] ] LLamb[6]=0.520049; and we get both s2s0 positive: s2s0[5]= s2s0[5]; s2s0[6]=41.6262 (-> M2[6]=Sqrt[s2s0[6]] LLamb[6] = 0.150); all f1ss0 and f2ss0 are positive *)

$lnsm=Log[10.^40];
$wwmax=10.^12;
$wmax=10.^8;
$wmax2d=10.^8;
$wmaxJ=10.^10.;
(*----------------------------------------------------------------------*)
(* nq= number of active quark flavors *)
be0[nq_]:=(11 -2 nq/3)/4;
 be1[nq_]:=(102-38 nq/3)/16;
 c1[nq_]:=(102-38 nq/3)/(11 -2 nq/3)/4;

(* be0 = 9/4; c1 = 16/9; *)
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 );
(* be3MS[nq_]:=(29243.0 - 6946.30 nq + 405.089 nq^2 + 1.49931 nq^3)/256; *)
 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 our c2-Lambert RSch beta functiuon: *)
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 is 4-loop MSbar beta function: *) 
betaMS4l[x_,nq_]:=- be0[nq] x^2 (1 + c1[nq] x + c2MS[nq] x^2 + c3MS[nq] x^3);
(*--------------------------------------------------------------------*)
(*--------------------------------------------------------------------*)
(* Now calculation of tkNf[s] = {tilde k}_N(s): *)

	       (* For this we first need: the coefficients k1[m],k2[m],k3[m],k4[m] which allow us to relate {tilde A_m} with (A_1)^m, for m not necessarily integer: *)	       
(* Note: m=1,2,3,...; n can be real noninteger: n=s: *)
(* ES[m,n]=S_m(n) = Sum_{k=1^n}(1/k^m) *)
ES[m_,n_]:= ((-1)^(m-1)/(m-1)!)*( PolyGamma[m-1,n+1] - PolyGamma[m-1,1] );
(* m=1,2,3,... ; Note: PolyGamma[0,x]=PolyGamma[x] *)
 
(* 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];
(* The above Z-ES relations I obtained by my own calculations *) 
(*--------------------------------------------------------*)
 k1f[nq_,s_]:=c1[nq] s (-1 + HarmonicNumber[s]);

k2f[nq_,s_]:= ( ((s-1) s/2) c22d[nq] + c1[nq]^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[nq_,s_]:= ( (s(s-1)(s+4)/12) c32d[nq] + (s(s+1)(s+2)/2) c1[nq] c22d[nq] ( (s/(s+2)) ES[1,s+1] + (2/(s+1)) + (1/(s+2)) - (11/6) ) + (s(s+1)(s+2)/6) c1[nq]^3 ( Z3[s+2] - 3 Z2[s+1] + 3 ES[1,s] - 1 ) );

k4f[nq_,s_]:=((s(s+1)(s+2)(s+3))/6) ( c42d[nq] (1/12 - 2/((s+1)(s+2)(s+3))) + c22d[nq]^2 (5/12 - (3 s^2 + 6 s +1)/((s+1)(s+2)(s+3))) + (c1[nq] c32d[nq]/2) ( (1 - 6/((s+2)(s+3))) Z1[s+1] + 4/(s+1) - 1/(s+2) - 1/(s+3) - 13/6 ) + (c1[nq]^2 c22d[nq]/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 ) + (c1[nq]^4/4) ( Z4[s+3] - 4 Z3[s+2] + 6 Z2[s+1] - 4 Z1[s] + 1 ) );

tk1f[nq_,s_]:=-k1f[nq,s];
tk2f[nq_,s_]:=k1f[nq,s] k1f[nq,s+1]-k2f[nq,s];
tk3f[nq_,s_]:=-k1f[nq,s] k1f[nq,s+1] k1f[nq,s+2] + k1f[nq,s] k2f[nq,s+1] + k1f[nq,s+2] k2f[nq,s] - k3f[nq,s];
tk4f[nq_,s_]:=k1f[nq,s] k1f[nq,s+1] k1f[nq,s+2] k1f[nq,s+3] - k1f[nq,s] k1f[nq,s+1] k2f[nq,s+2] - k1f[nq,s] k1f[nq,s+3] k2f[nq,s+1] - k1f[nq,s+2] k1f[nq,s+3] k2f[nq,s] + k2f[nq,s] k2f[nq,s+2] + k1f[nq,s] k3f[nq,s+1] + k1f[nq,s+3] k3f[nq,s] - k4f[nq,s];
(*----------------------------------------------------------------------*)
(*----------------------------------------------------------------------*)
zu2d[nq_,u_, fi_] := - (1/c1[nq]/E) u^(-be0[nq]/c1[nq]) Exp[- I*(be0[nq]/c1[nq]) fi];
(* u=|Q2|/L22[nq] = |Q^2|/LL2[nq], where LL2[nq]=(LambertLambda[nq])^2 *)

 a2l[nq_,u_, fi_] := If[fi >= 0, -(1/c1[nq]) 1/(1 + ProductLog[-1, zu2d[nq,u, fi]]), -(1/c1[nq]) 1/(1 + ProductLog[1, zu2d[nq,u, fi]])];

a2d[nq_,u_, fi_] := If[fi >= 0, -(1/c1[nq]) 1/(1 -c22d[nq]/c1[nq]^2 + ProductLog[-1, zu2d[nq,u, fi]]), -(1/c1[nq]) 1/(1 - c22d[nq]/c1[nq]^2 + ProductLog[1, zu2d[nq,u, fi]])];
(* u = |Q^2|/LL2[nq] > 0; Q^2=Q2 = LL2[nq] Exp[I fi]; Q^2/L[nq]^2 = u Exp[I fi] *)
 (* a2d[nq,u,fi] =alpha_s[Q2;nq]/Pi is the solution with betaP11[x,nq]= - be0[nq] x^2 (1 + (c1[nq] - c22d[nq]/c1[nq]) x)/(1 - (c22d[nq]/c1[nq]) x), i.e., an effective "3-loop" solution, in the Lambert RSch with a chosen c22d[nq]=c22din (-> c3=c22d[nq]^2/c1[nq], c4=c22d[nq]^3/c1[nq]^2, c5=c22d[nq]^4/c1[nq]^3,  c6=c22d[nq]^5/c1[nq]^4; etc. *)

 (* z=u Exp[I fi] = Q^2=Q2 = ; zu2d[nq,u,fi]=zf[nq,z]=zf[nq,u Exp[I fi]]; *)
 (* a2d[nq,u,fi] = af[nq,z] = af[nq,u Exp[I fi]]: *)
  zf[nq_,z_]:= - (1/c1[nq]/E) z^(-be0[nq]/c1[nq]);  
 af[nq_,z_]:= If[Arg[z] >= 0, -(1/c1[nq]) 1/(1 -c22d[nq]/c1[nq]^2 + ProductLog[-1, zf[nq,z]]),-(1/c1[nq]) 1/(1 - c22d[nq]/c1[nq]^2 + ProductLog[1, zf[nq,z]])]; (* valid only for fi>= 0 *);
  (* z = u Exp[I fi]; u=|z| = absQ2/LL2[nq]  >0; fi = Arg[z]=Abs[Q2] *)

(* Beware: ProductLog[1,z] cannot be applied for z real, *)
  (* z has to be in the lower half plane; *)
  (* ProductLog[-1,z] can be applied for z in the upper half AND for z real *)
(*----------------------------------------------------------------------_*)
(*----------------------------------------------------------------------_*)
(*----------------------------------------------------------------------_*)
  (*-------------------------------------------------------*)
 (* We have here: c22d[nq] < 0. As a consequence, we have a (Landau) pole for Q^2>0, i.e., for zu2d[nq,u,0]<0 (at u=uL[nq]>0): -1/e < zL[nq]=zu2d[nq,uL[nq],0] < 0 *)
  (* Let us find uL[nq] = QL^2/L[nq]^2, i.e. zL[nq]: *)

 sL[nq_]:=c1[nq]^(-c1[nq]/be0[nq]); (* sL =approx 0.634696 for nq=3; sL such that: tr12d[s]=0 for s<-sL *)
 Do[Module[{},If[c22d[nqu]<0,solzL[nqu]= FindRoot[-1 + c22d[nqu]/c1[nqu]^2==ProductLog[-1,zu2d[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}];

 (* Note, for nqu=6: uL is quite large (7.742 >> sL[6]=1.04011) *)

   (* The dispersion relation for a(perturbative) has to include a [circular, of small radius sg0=dsL*L^2 in Q^2 plane] integration around the pole at Q^2=L0^2=sL L^2, because that pole is quite strong in the n-loop case: *)
  (*--------------------*)
 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, zu2d[nq,-s, 0]])]],0];
 (* tr12d[nq,s]=rho_1(sigma = s L[nq]^2;nq) = Im[a2d[nq,s, -Pi] = Im[a2d[Q^2=-s L[nq]^2 - I eps;nq]], when s=sigma/L[nq]^2 > 0; *)
    (* s = u Exp[- I Pi] if s > 0; *)
  (* tr12d[nq,s] is applicable even when s<0; s= -u Exp[I 0] if s < 0 *)
  (* Note that tr12d[nq,-sL[nq]]=0; tr12d[nq,0]=0 *)
  (* tr12d[nq,s]=0 for all s<-sL exactly *)
(*-----------------------------*)

  (* Now the dispersively calculated a2d[nq,u,fi] = a2d[Q^2=(u*L2)Exp[I fi];nq]  (u>0) *)
(* Beware: the dispersively calculated a has to include careful integration around s=-sL[nq] and s=-uL[nq] - integration along the two cicles of small radius sg0=dsL[nq]*L[nq]^2 in Q^2 plane *) 

   small=10.^-10;
   (* Note: HeavisideTheta2d[-c22d[nq]-small] = 0 for c22d[nq]>0 and c22d[nq]=0;  =1 for c22d[nq] < -small *)
  dsL[nq_]:=Min[Abs[(uL[nq]-sL[nq])]/10.,10.^-6];
    (*-----*)
  (* adis2d[nq,u,fi] is dispersively calculated a2d[nq,u,fi]: *)
  adis2d[nq_,u_,fi_]:=(1/Pi) ( NIntegrate[tr12d[nq,s]/(s + u Exp[I fi]),{s,-sL[nq]+dsL[nq],0}]+NIntegrate[tr12d[nq,s]/(s + u Exp[I fi]),{s,0,1.}] + NIntegrate[tr12d[nq,Exp[w]]/(1 + u Exp[I fi] Exp[-w]),{w,0.,Infinity}, MaxRecursion->30,WorkingPrecision->16]) + (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]/(u Exp[I fi] - 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]/(u Exp[I fi] - uL[nq] - dsL[nq] Exp[I Phi]),{Phi,-Pi,Pi}];
   (* adis2d[u,fi]=a2d[nq,u,fi]; except on the cut, i.e., for u < sL[nq] ; and at u=uL[nq] *)
    (*-----*)
  (* Now adis2d2[nq,u,fi] ia dispersively calculated by using the residue formula for integral around s=-uL[nq] (when c22d[nq] <-small): *)
   adis2d2[nq_,u_,fi_]:=(1/Pi) ( NIntegrate[tr12d[nq,s]/(s + u Exp[I fi]),{s,-sL[nq]+dsL[nq],0}]+NIntegrate[tr12d[nq,s]/(s + u Exp[I fi]),{s,0,1.}] + NIntegrate[tr12d[nq,Exp[w]]/(1 + u Exp[I fi] Exp[-w]),{w,0.,$wmaxJ}]) + (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]/(u Exp[I fi] - sL[nq] - dsL[nq] Exp[I Phi]),{Phi,-Pi,Pi}] + HeavisideTheta[-c22d[nq]-small] (1/(u Exp[I fi] - uL[nq])) NResidue[af[nq,z],{z,uL[nq]}];
  (* Note, adis2d has:  MaxRecursion->30,WorkingPrecision->16, which makes it better than the (below) adis2d2 and adis2simpl; importane at least in Mathematica 11.3 *)
  (*-----*)
   (* adis2simpl is as adis2, but integration around s=-sL[nq] is done in the simple manner, in the hope that the singularity at sL[nq] is integrable directly *)   
adis2d2simpl[nq_,u_,fi_]:=(1/Pi) ( NIntegrate[tr12d[nq,s]/(s + u Exp[I fi]),{s,-sL[nq],0}]+ NIntegrate[tr12d[nq,s]/(s + u Exp[I fi]),{s,0,1.}] + NIntegrate[tr12d[nq,Exp[w]]/(1 + u Exp[I fi] Exp[-w]),{w,0.,$wmaxJ}]) + HeavisideTheta[-c22d[nq]-small] (1/(u Exp[I fi] - uL[nq])) NResidue[af[nq,z],{z,uL[nq]}];

   (* Note, when uL[nq]-sL[nq] < 0.01, then the numerical calculation of residue (NResidue) in z=uL[nq] is not relable any more; therefore, better NOT use NResidue, but rather the contour integration in circle around Q2/L2=uL[nq] -- i.e., adis; not: adis2, or adissimpl *)
    (*----------------------*)
  A1APT[nq_,u_,fi_]:=(1/Pi) NIntegrate[tr12d[nq,s]/(s + u Exp[I fi]),{s,0,1.}] + (1/Pi) NIntegrate[tr12d[nq,Exp[w]]/(1 + u Exp[I fi] Exp[-w]),{w,0.,$wwmax}];
   (*    wmax=10.^3; *)
    (* Numerically, |adis-a| ~ 10^-10 - 10^-9 should be, IF a2d[nq,u,fi] has no Landau singularities outside the real Q^2-axis *)
   (* The numerical precision of evaluation of adis2d[u,0] is ~sim 10^-9 *)
(*------------------------------------------------*)
(*--------------------------------------------------------------------*)
   (* THE BASIC IDEA IS: TO HAVE c22din as input parameter; and also Nf as input; first the calculation is performed at Nf0=3, using <<s0r.m to obtain s0[Nf=3] for this c22din; then keeping the same c22din, adjust s0[Nf] (Nf=4,5,6) so that at the thresholds 2*mq+-0 we get foe A1 the same values of discontinuities as in pQCD (CKS) *)
(*--------------------------------------------------------------------*)
MZ=91.1876;
mt=163.;
mb=4.2;
mc=1.27;
  (* mb=mbbar(mbbar^2);  mc=mcbar(mcbar^2); mt=mtbat(mtbar^2); *)
(*--------------------------------------------------------------------*)
  (* kappa2d; fixed by default *)
(* kappa2d=2 will be now the only possibility here: *)
(* Q_{thr} = kappa2d mq (q=b,c) in RGErun.m *)
 kappa2d=2;
 (* kappa2d=2; *)
 (*------------------------------------------------*)
(* mb=mbbar(mbbar); mc=mcbar(mcbar); thresholds at kappa2d*mb and kappa2d*mc *)
  (* kappa2d=2 *)

kappa2d=2.;
  (*--------------------------------*)
  (* First we want to obtain ainMS[3]=a(MSbar,mu=kappa2d*mc;Nf=3): *)
  (* Now (4-loop MSbar beta) RGE-running down from mu=MZ to mu=kappa2d*mc(=Q4), making (3-loop) threshold matchings at mu=kappa2d*mb(=Q5) from Nf=5 to Nf=4; and at mu=kappa*mc(=Q4) from Nf=4 to Nf=3 *)

QMZ=MZ;
aaMZ=alphaMSMZ/Pi;
(* aaMZ=a(QMZ^2) *)
   Q5=kappa2d mb; (* Q4=mu_thr(nf=5->4) *)
   Q4=kappa2d mc; (* Q3=mu_thr(nf=4->3) *)
 (*--------------*)
  (* First 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[QMZ/Q5];
del=0.01;

sol1MS = NDSolve[{y'[x]==Bta[y[x]], y[xx5]==aaMZ}, y, {x,-del,+del}];
aa=First[y[0.]/.sol1MS];
aaQ5plus=aa;
(* aa ' is now a(Q5+0)=a(MSbar,kappa2d*mb+0;nff=5) *)
 (*---------------*)
 (* 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;
  (* aaQ5low = a(Q5-0;nf=5-1)=a(MSbar,kappa2d*mb-0,nf=4) *)
 (*---------------*)
 (* Now the running in the nf=4 regime from Q5=kappa*mb to Q4=kappa*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/Q4];
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,kappa*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;
aaQ4low=aa;
  (* aaQ4low = a(Q4-0;nf=4-1)=a(MSbar,kappa2d*mc,nf=3) *)
  (*------------------------------------*)
  ainMS[3]=aaQ4low;  (* ainMS[3]=0.0844636; *)
  (* This is: ainMS[3]= a(MSbar,Q=kappa2d*mc-0, with nf=3); obtained by RGE running down with 4-loop MSbar beta function, and 3-loop threshold matching at Qthr=kappa2d*mq (mq=mb, mc) *)

   (* Note that in RGErundown.m the thresholds (CSK) are implemented with running DOWNWARD, that gives somewhat different results than in the case when thresholds are implemented in the RGE-running upwards *)
  ainMS[4]=aaQ4plus;  (* ainMS[4]=0.0870619 *)
  (* ainMS[3]=a(MSbar,Q=kappa2d*mc-0, nf=3) *)
  (* ainMS[4]=a(MSbar,Q=kappa2d*mc+0, nf=4) *)
 (*-----------------------------------------------------------------*)
   (* Now the running in the nf=3 regime from Q4=kappa*mc to Q (Q =< kappa*mc): *)
   (* We can assume that 0 =< Q^2 =< (kappa*mc)^2, i.e., the nf=3 regime is approximately the regime of three massless quarks *)

nff=3;
Clear[Bta];
Bta[z_]:=- be0[nff] z^2 (1 + c1[nff] z + c2MS[nff] z^2 + c3MS[nff] z^3);


  xx3=0;  (* xx=2*Log[Q/Q4] *)
delext=20.;
sol1a = NDSolve[{y'[x]==Bta[y[x]], y[0]==aaQ4low}, y, {x,xx3-delext,+delext}];
aaMSQ=First[y[xx3]/.sol1a];

aMSnf3[Qf2_]:=First[y[Log[Qf2/Q4^2]]/.sol1a];
(* aaMSQ = a(Q^2;MSbar; nf=3); we ' assume that Q =< kappa*mc *)
  (* aMSnf3[Q2] is a[Q2] in MSbar scheme, nf=3, where Q2 is in GeV^2 *)
aMSnf3plot[Q2_? NumberQ]:=If[Q2>0.3706,aMSnf3[Q2],10.];
(* aMSnf3plot[Q2] = a(Q2) in nf=3 MSbar scheme, for Q2>0 *)
  (*----------------------------------------------------------_*)
  (* Now the (4-loop MSbar) RGE-running UPWARDS from Q=MZ to kappa2d*mt, and beyond: *)
  (* we will obtain in this way: ainMS[5]=a[(kappa2d*mt)^2-0;nq=5;MSbar];ainMS[6]=a[(kappa2d*mt)^2+0;nq=6;MSbar] *)  

QMZ=MZ;
aaMZ=alphaMSMZ/Pi;
Q6=kappa2d mt;
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[QMZ/Q6];
del=0.01;
sol1MS = NDSolve[{y'[x]==Bta[y[x]], y[xx5]==aaMZ}, y, {x,-del,+del}];
aa=First[y[0.]/.sol1MS];
aaQ6low=aa;
(* aa ' is now a(Q6-0)=a(MSbar,kappa2d*mt-0; nff=5) *)
 (*---------*)
  (* 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;
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;
  (* Note, the new aaQ6plus=(MSbar,Q6+0, nnf+1), the old aa was a(MSbar, Q6-0, nnf) *)
 (*------------------*)
 ainMS[5]=aaQ6low;
 ainMS[6]=aaQ6plus;
  (* ainMS[5]=a(MSbar,Q=kappa2d*mt-0, nf=5) *)
  (* ainMS[6]=a(MSbar,Q=kappa2d*mt+0, nf=6) *)
  (*---------------------*)
  (* Now we calculate ain[3]=al(ourRSch;Q=kappa2d*mc-0;nf=3)/Pi *)
  (* and ain[4]=al(ourRSch;Q=kappa2d*mc+0;nf=4)/Pi *)
  (* and we calculate ain[5]=al(ourRSch;Q=kappa2d*mt-0;nf=5)/Pi *)
  (* and ain[6]=al(ourRSch;Q=kappa2d*mt+0;nf=6)/Pi *)

  ain0=ainMS;
  (* 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;*) 
  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]); 
 eepss=10.^(-7);
  aux1[yyy_,nq_]:=1/yyy + c1[nq] Log[c1[nq] yyy/(1 + c1[nq] yyy)];
(* aux1[ain] + NIntegrate[fl1[xx], {xx, eepss, ain}] - rhs1 = 0 has to be *)
  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}];
  (* ain[3]=a(Q^2=(kappa2d mc)^2-0; nf=3; c22d[nf],c3,c4,...) pQCD, in RSch of our betaP11[a,3]; we will use the same LLamb[3]=Lambda(nf=3) that is obtained from this a(Q^2;c22d[nf],...;nf=3) also in our A1(Q^2;c22d[nf],...;nf) *)
  (* ain[4]=a(Q^2=(kappa2d mc)^2+0; nf=4; c22d[nf],c3,c4,...) pQCD, in RSch of our betaP11[a,4]; *)
  (* we note that ain[3] and ain[4] are related by the CKS 3-loop matching (PRL79,2184 Eqs.(6)-(7)); in RGErundown.m, the matchings were used by calling <<RGE3lmatchdown.m *)
  (*-------------------------------------------*)
  (* 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)^2/0.26^2; sol0=FindRoot[a2d[nqu,uub,0]==ain[nqu],{uub,ub0}]; ub[nqu]=uub/.sol0; LL2[nqu]=(kappa2d*mc)^2/ub[nqu]; LLamb[nqu]=Sqrt[LL2[nqu]]],{nqu,3,4}];

  Module[{}, ub0=(kappa2d*mt)^2/0.097^2; sol0=FindRoot[a2d[5,uub,0]==ain[5],{uub,ub0}]; ub[5]=uub/.sol0; LL2[5]=(kappa2d*mt)^2/ub[5]; LLamb[5]=Sqrt[LL2[5]]];

 Module[{}, ub0=(kappa2d*mt)^2/0.03^2; sol0=FindRoot[a2d[6,uub,0]==ain[6],{uub,ub0}]; ub[6]=uub/.sol0; LL2[6]=(kappa2d*mt)^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}];

  apert[nq_,Q2_,fi_]:=a2d[nq,Q2/LL2[nq],fi];  (* = a[Q^2=Q2 Exp[I fi]], the underlying pQCD coupling of 2d anQCD *)
 (*----------------------------------------------------------------*)
 (*----------------------------------------------------------------*)
  (* Now solving equations for A1 - a ~ (L2/Q^2)^5, considering s0 as input: *)
     (* We use notation: int_{-sL[nq]}^s0  ds s^j tr12d[nq,s] = Pj (j=0,1,2,3) *)

     (* 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];
     (* Notation: fjs = fj^2 (j=1,2) *)

     (* There are two solutions: *)

     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]];

     (* However, the solution 2 is equivalent to solution 1, because: f1so1=f2so2, f2so1=f1so2; s1o1=s2o2; s2o1=s1o2 *)
  (*---------------------------------*)
  (*----------------------------------*)
   (* Note: correction in the formulas below: sL[nq]^k -> (-sL[nq])^k; uL[nq]^k -> (-uL[nq])^k; for k=0,1,2,3; Correction made on March 16, 2012 *)

  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}]; 

   (* The numerical precision of calculation of Pjs0 is ~10^-12; *)
(* This is IMPORTANT to keep in mind *)

  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}];
     (* This is the solution (as numbers), for a given numerical value of s0 *)
     (* Note: f1ss0 = f1^2 for this s0; f2ss0 = f_2^2 for this s0 *)
     (* s1s0 = s1 for this s0; s2s0 = s2 for this s0 *)

  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}];
  (*---------------------------------------------------------------*)
    (* w = ln(s) = ln(sigma/LLamb[nq]^2); s = sigma/LLamb[nq]^2 = Exp[w] *)
    (* Now constants needed for an approximate expression for A1an[u,fi], for u < umax (~ 10^6 ?): *) 
    sm=10.^40; lnsm=Log[sm];
  Do[Module[{},J0[nq]=(1/Pi) NIntegrate[tr12d[nq,Exp[w]],{w,lnsm,Infinity},MaxRecursion->30,WorkingPrecision->16]; J1[nq]=(1/Pi) NIntegrate[tr12d[nq,Exp[w]] Exp[-w],{w,lnsm,Infinity},MaxRecursion->30,WorkingPrecision->16]; J2[nq]=(1/Pi) NIntegrate[tr12d[nq,Exp[w]] Exp[-2 w],{w,lnsm,Infinity},MaxRecursion->30,WorkingPrecision->16]; J3[nq]=(1/Pi) NIntegrate[tr12d[nq,Exp[w]] Exp[-3 w],{w,lnsm,Infinity},MaxRecursion->30,WorkingPrecision->16]],{nq,3,6}];
 (* For J0[nq],MaxRecursion->30,WorkingPrecision->16 is very IMPORTANT to get the true result (nonzero), at least this is important in Mathematica 11.3 *)
  (*-----*)

   A1[nq_,Q2_,fi_]:= ( (1/Pi) NIntegrate[tr12d[nq,Exp[w]]/(1 + (Q2/LL2[nq]) Exp[I fi] Exp[-w]),{w,lns0[nq],lnsm},MaxRecursion->30,WorkingPrecision->16] + 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] );
 (* Note: A1(Q^2)=A_{1}^{(an2d)}(Q^2) is the running coupling of 2d anQCD *)
 (* A1[Q^2=Q2 Exp[I fi]; Nf)=A1an[Nf,Q2,fi]; note: Q2=|Q^2|; fi=Arg[Q^2] *)
  (* Note; |Q^2|=Q2 is input variable for analytic couppligs; while for the (underlying) pQCD couplings a2d the input variable is instead: u = Q2/LL2[nq] = |Q^2|/LambertLambda[nq]^2 *)

    (* Apparently work for |u| < 10.^38 or so; and any fi *)

  diffan[nq_,Q2_,fi_]:=A1[nq,Q2,fi]-a2d[nq,(Q2/LL2[nq]),fi];
  reldiffan[nq_,Q2_,fi_]:=A1[nq,Q2,fi]/a2d[nq,Q2/LL2[nq],fi]-1;
diff[nq_,Q2_,fi_]:=adis2d[nq,Q2/LL2[nq],fi]-a2d[nq,Q2/LL2[nq],fi];
 reldiff[nq_,Q2_,fi_]:=adis2d[nq,Q2/LL2[nq],fi]/a2d[nq,Q2/LL2[nq],fi]-1; 
  (*  reldiff[nq_,Q2_,fi_]:=diff[nq,Q2,fi]/a2d[nq,Q2/LL2[nq],fi]; *)
  (* u = |Q2|/LL2[nq] *)
    (*----------------------------------------------------_*)
    (*----------------------------------------------------_*)
  
    (* Now calculation of higher power analogs: *)

  (* Calculation of our {tilde A2}^{analytic}[Q2,fi], in our model, *)
    (* where Q2=|Q^2| and fi=Arg[Q^2]: *)
  tA2an0[nq_,Q2_,fi_]:=(1/be0[nq]/Pi) ((Q2/LL2[nq]) Exp[I fi]) ( NIntegrate[ Exp[-w] tr12d[nq,Exp[w]]/(1 + (Q2/LL2[nq]) Exp[I fi] Exp[-w])^2,{w,lns0[nq],lnsm}] + (f1ss0[nq]*Pi)/((Q2/LL2[nq]) Exp[I fi]+s1s0[nq])^2 + (f2ss0[nq]*Pi)/((Q2/LL2[nq]) Exp[I fi]+s2s0[nq])^2 );
    (* Now approximate form: *)
 tA2[nq_,Q2_,fi_]:=(1/be0[nq]/Pi) ((Q2/LL2[nq]) Exp[I fi]) ( NIntegrate[ Exp[-w] tr12d[nq,Exp[w]]/(1 + (Q2/LL2[nq]) Exp[I fi] Exp[-w])^2,{w,lns0[nq],lnsm},MaxRecursion->30,WorkingPrecision->16]  + (f1ss0[nq]*Pi)/((Q2/LL2[nq]) Exp[I fi]+s1s0[nq])^2 + (f2ss0[nq]*Pi)/((Q2/LL2[nq]) Exp[I fi]+s2s0[nq])^2 + J1[nq] - 2  (Q2/LL2[nq]) Exp[I fi] J2[nq] + 3 (Q2/LL2[nq])^2 Exp[I 2 fi] J3[nq] );

  
  (* Calculation of our {tilde A3}^{analytic}[Q2,fi], in our model, *)
  (* where Q2=|Q^2| and fi=Arg[Q^2]: *)
  tA3[nq_,Q2_,fi_]:=(1/2/be0[nq]^2/Pi)((Q2/LL2[nq]) Exp[I fi]) ( NIntegrate[ Exp[-w] tr12d[nq,Exp[w]]/(1 + (Q2/LL2[nq]) Exp[I fi] Exp[-w])^2,{w,lns0[nq],lnsm},MaxRecursion->30,WorkingPrecision->16] -2  NIntegrate[ Exp[- w] tr12d[nq,Exp[w]]/(1 + (Q2/LL2[nq]) Exp[I fi] Exp[-w])^3,{w,lns0[nq],lnsm},MaxRecursion->30,WorkingPrecision->16] + (f1ss0[nq]*Pi) (-s1s0[nq] + (Q2/LL2[nq]) Exp[I fi])/(s1s0[nq] + (Q2/LL2[nq]) Exp[I fi])^3  + (f2ss0[nq]*Pi) (-s2s0[nq] + (Q2/LL2[nq]) Exp[I fi])/(s2s0[nq] + (Q2/LL2[nq]) Exp[I fi])^3 );

  (* Calculation of our {tilde A4}^{analytic}[nq,Q2,fi], in our model, *)
    (* where Q2=|Q^2| and fi=Arg[Q^2]: *)
  tA4[nq_,Q2_,fi_]:=(1/6/be0[nq]^3/Pi)((Q2/LL2[nq]) Exp[I fi]) ( NIntegrate[ Exp[-w] (1 - 4 Exp[-w] (Q2/LL2[nq]) Exp[I fi] + Exp[-2 w] (Q2/LL2[nq])^2 Exp[2 I fi]) tr12d[nq,Exp[w]]/(1 + (Q2/LL2[nq]) Exp[I fi] Exp[-w])^4,{w,lns0[nq],lnsm},MaxRecursion->30,WorkingPrecision->16] + (f1ss0[nq]*Pi) (s1s0[nq]^2 - 4 s1s0[nq] (Q2/LL2[nq]) Exp[I fi] + (Q2/LL2[nq])^2 Exp[2 I fi] )/(s1s0[nq] + (Q2/LL2[nq]) Exp[I fi])^4  + (f2ss0[nq]*Pi) (s2s0[nq]^2 - 4 s2s0[nq] (Q2/LL2[nq]) Exp[I fi] + (Q2/LL2[nq])^2 Exp[2 I fi] )/(s2s0[nq] + (Q2/LL2[nq]) Exp[I fi])^4  );

 (* tAjan[nq,Q2,fi] = {tilde A}_j(Q^2;nq) with Q^2 = Q2 Exp[I fi] *) 
   (*----------------------------------------------------_*)
(*--------------------------------------------------*)
  (*----------------------------------------------------_*)
(*--------------------------------------------------*)
  (*----------------------------------------------------_*)
(*--------------------------------------------------*)
  (* Now calculation of rtau: *)
mtau2=1.77686^2;   
 (*-----------------------*)
bet2MS= (2857/2 - 5033 nfq/18 + 325 nfq^2/54)/64;
bet2MS=bet2MS/.{nfq->- 6 bet0 + 33/2};  
(* This acts only if bet2 written in terms of nfq *)
     bb20 = bet2MS/.{bet0->0};
     bb21=Coefficient[bet2MS,bet0^1];
     bb22=Coefficient[bet2MS,bet0^2];
  
  CC=0.;   (* CC= C2 = C3 is RScl parameter = Log[mu^2/Q^2] *) 

b11=19/4; b10=-107/16; (* be1 = b11 bet0 + b10 *)
	       bet2 = c22d[3] be0[3];
	       bet3 = c32d[3] be0[3];
(* Now for Adler function dv(Q^2), in RScl mu^2=Q^2 and bMS:('b' scheme); *)
(* note: bcvk[i,j] = {\bar c}^{(k)}[i,j] for dv(Q^2), in 'b' scheme *)
	       bcv1[1,0]=1/12; bcv1[1,1]=(11/2 - 4 Zeta[3]);
	       bcv1[2,0]=-(12521/288) + 13 Zeta[3];
               bcv1[2,1]=1/12 (401 - 424 Zeta[3] + 200 Zeta[5]);
	       bcv1[2,2]=302/9 - (76 Zeta[3])/3;
  d3bMS=( nqf^3 (-6131/5832 + 203 Zeta[3]/324 + 5 Zeta[5]/18) + nqf^2 (1045381/15552 - 40655 Zeta[3]/864 + 5 (Zeta[3]^2)/6 - 260 Zeta[5]/27) + nqf (-13044007/10368 + 12205 Zeta[3]/12 - 55 (Zeta[3])^2 + 29675 Zeta[5]/432 + 665 Zeta[7]/72) + (144939499/20736 - 5693495 Zeta[3]/864 + 5445 (Zeta[3]^2)/8 + 65945 Zeta[5]/288 - 7315 Zeta[7]/48) ); (* see Eq.(4) of PRL101,012002 (2008) *)
d3bMS=d3bMS/.{nqf->- 6 bet0 + 33/2};
  	       bcv1[3,0]=d3bMS/.{bet0->0};
               bcv1[3,1]=Coefficient[d3bMS,bet0];
               bcv1[3,2]=Coefficient[d3bMS,bet0^2];
               bcv1[3,3]=Coefficient[d3bMS,bet0^3];
(* The exact values of all bcv1[3,j] incorporated on: Oct.21,2009 *)
(* bcv1[3,2]=49.5237; bcv1[3,3]= 2.18004; known earlier exactly; *)
(* bcv1[3,0]=46.1992; bcv1[3,1]=-131.04; were approximate values by PMS and ECH *)
(* Now the exact values are: bcv1[3,0]=32.727; bcv1[3,1]=-115.199 *)
	       bcv2[1,0]=(bcv1[2,0]-b10 bcv1[1,1])/bcv1[1,0];
               bcv2[1,1]=(bcv1[2,1]-b11 bcv1[1,1])/bcv1[1,0];
               bcv2[2,0]=(bcv1[3,0] - bb20 bcv1[1,1])/bcv1[1,0];
               bcv2[2,1]=(bcv1[3,1] - (5/2) b10 bcv1[2,2] - bb21 bcv1[1,1])/bcv1[1,0];
               bcv2[2,2]=(bcv1[3,2] - (5/2) b11 bcv1[2,2] - bb22 bcv1[1,1])/bcv1[1,0];
  (*------------------------*)
		   d1bar= bcv1[1,0]+be0[3] bcv1[1,1];
                   d2bar= bcv1[2,0]+be0[3] bcv1[2,1] + be0[3]^2 bcv1[2,2];
		   d3bar=d3bMS/.{bet0->be0[3]}; (* i.e., numerical value of d2bMS *)
	       (* Now we do the change to {tilde d}_{j} = tdj, in {bar MS} and mu^2=Q^2: *)
   td1bar=d1bar; td2bar=d2bar - c1[3] d1bar;
    td3bar=d3bar - (5/2) c1[3] d2bar + ((5/2)c1[3]^2 - c2MS[3]) d1bar;

  (* Now we perform the change of RSch from {bar MS} to our (c2,c3,...) and of RScl from Q^2 to mu^2 (CC=ln(mu^2/Q^2): tdjbar->tdj: *)
 
     td1=td1bar+be0[3] CC;
     td2=td2bar + 2 be0[3] td1bar CC + be0[3]^2 CC^2 - (c22d[3] - c2MS[3]);
     td3=( (td3bar + 3 be0[3] td2bar CC + 3 be0[3]^2 td1bar CC^2 + be0[3]^3 CC^3) + (-3 (td1bar + be0[3] CC) + (5/2) c1[3])(c22d[3] - c2MS[3]) - (1/2)(c32d[3] - c3MS[3]) );

		   (* Note: Adler function is now: D[Q^2]=A1 + td1 {tilde A}_2 + td2 {tilde A}_3 + td3 {tilde A}_4, where everything in our RSch (c2, c3, ...) and in mu^2 = Q^2 Exp[CC] *)

		(* Now d1, d2, d3 in our RSch (c2,c3,...) and RScl mu^2=Q^2 Exp[CC]: *)
		d1=td1;d2=td2+c1[3] td1; d3=td3 + (5/2) c1[3] td2 + c22d[3] td1;

   (* Now the cv1[n,n] in the RScl mu^2=Q^2 Exp[CC] (where Q^2=mtau^2 Exp[I fi] in the contour integrations for rtau): *)
		   cv1[1,1]=bcv1[1,1] + CC;
		   cv1[2,2]=bcv1[2,2] + 2 CC bcv1[1,1] + CC^2;
		   cv1[3,3]=bcv1[3,3] + 3 CC bcv1[2,2] + 3 CC^2 bcv1[1,1] + CC^3;
		   (* Now calculation of Tj s; we should have: Tj=ttv{j+1}: *)
		   T1=td1 - be0[3] cv1[1,1];
		   T2=td2 - be0[3]^2 cv1[2,2];
		   T3=td3 - be0[3]^3 cv1[3,3];

   td1bLB=T1; td2bLB=T2; td3bLB=T3;
(* We have (c22din=-4.9): td1=1.63982, T1=0.0833333; td2=12.8268, T2=-2.88436; td3=34.5787, T3=9.74666;  this indicates that LB is a relatively good approximation, at least up to ~a^4, for rtau in the Lambert RSch (with c2=-4.9) *)
(*------------------------------------------------------------*)
(*-----------------------------------------------------*)
		(* Now calculating rtau using [2/2] diagonal Baker-Gammel (dBG) for Adler: *)

  (* First the diagonal Baker-Gammel (dBG) of DAdlan in 4danQCD: *)
   dvmptser[xx_]:=xx + td1 xx^2 + td2 xx^3 + td3 xx^4;
		(* Adler is D[Q2]= A1[Q2]+ td1 tA2[Q2] + td2 tA3[Q2] + td3 tA4[Q2]; note, here and from now on we assume that eCC=1, i.e., mu^2=Q^2 *)

   (* This dvmptser[xx] is the AUXILIARY power series -> leads via dPA to dBGA: *)
 Needs["NumericalCalculus`"];

dPAdvmpt[nn_,xx_]:=PadeApproximant[SetPrecision[dvmptser[xxp],30],{xxp,0,nnp}]/.{nnp->nn,xxp->xx};


 (*-------------*)
   (* Now calculation of alt[j] and ut[j] (j=1,2), where: *)
   (* dPAdvmpt[n,x] = Sum_{j=1}^n alt[n,j] x/(1 + x ut[n,j] be0) = Sum_{j=1}^n alt[n,j] a(1-l.)(Q[n,j]^2), where Q[n,j]^2 = Q^2 Exp[ut[n,j]]; note, we have mu^2=Q^2 in our dv[Q^2]=dvtmpt[Q^2] *)

  (* We use: Sum_{j=1}^2 alt[n,j]/(1 + x ut[n,j] be0)= Apart[dPAdvmpt[n,x]/x] in order to identify alt[n,j] and ut[n,j]'s;  for n=1; and n=2; *)
(* ' *)

alt[1,1]=1.;
ut[1,1]=-1.6398212048969847647356898797/be0[3];  
(* this is for [1/1] dBG of Adler function *)

(* Apart[(A + B x)/((x - k1) (x - k2))]; use the formula obtained from: Apart[(A + B x)/((x - k1) (x - k2))] *)

AA=-Numerator[Simplify[dPAdvmpt[2, x]/x]]/.{x->0};
BB=-Coefficient[Numerator[Simplify[dPAdvmpt[2, x]/x]],x];
deno=NSolve[Denominator[Simplify[dPAdvmpt[2, x]/x]] == 0, x];		
kk1=x/.deno[[1]];
kk2=x/.deno[[2]];
alt[2,1]= (AA + BB kk1)/(kk1-kk2)/(-kk1);
alt[2,2]= (-AA - BB kk2)/(kk1-kk2)/(-kk2);
ut[2,1]=-1/(kk1)/be0[3];
ut[2,2]=-1/(kk2)/be0[3];

 (* this is for [2/2] dBG of Adler function *)

(* Cross-checking: *)
Appt[nn_,xx_] := Sum[alt[nn,jj]/(1 + xx ut[nn,jj] be0[3]), {jj, 1, nn}];
(* Then check: Together[Simplify[N[Appt[n,x] - Apart[dPAdvmpt[n,x]/x]]]] is practically zero, for n=1,2  *)
 (* Yes, I checked, high precision  ~ 10^-16 *)

  Do[Do[kapps[nn,jj]=Exp[ut[nn,jj]],{jj,1,nn}],{nn,1,2}];

(* We have, as a result (when CC(=Log[mu^2/Q^2]) = 0), i.e., eCC=1: *)
 (*  kapps[1,1]=0.4824830808572873476119276474; *)

(* when c22din=-4.9: *)
(*  alt[2,1]=0.35404279856468405851182235 *)
(*  alt[2,2]=0.6459572014353159414881777 *)
(*  kapps[2,1]=3.26299318382473852552078473 *)
(*  kapps[2,2]=0.169236389981400880666520744 *)

		(* D[Q2]=A1(Q2) + td1 {tilde A}_2(Q2) + td2 {tilde A}_3(Q2) + td3 {tilde A}_4(Q2) + O(A_5) = alt[2,1] A1(kapps[2,1]Q2) + alt[2,2] A1([kapps[2,2]Q2) + O(A5), where the last expression is [2/2] dBG *) 

(*---------*)
(* Beware: only correct if eCC=1: *)
(* theta_max=Pi-eeps taken *)
 eeps=10.^-4;

tauk22=(1/2/Pi) 2 Re [( (1/Pi) NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) tr12d[3,Exp[w]]/(1 + (kapps[2,2] mtau2/LL2[3]) Exp[I fi] Exp[-w]),{fi,0,+Pi-eeps},{w,lns0[3],lnsm},MaxRecursion->30,WorkingPrecision->16] +  NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) ( f1ss0[3]/( (kapps[2,2] mtau2/LL2[3]) Exp[I fi] +s1s0[3])  + f2ss0[3]/( (kapps[2,2] mtau2/LL2[3]) Exp[I fi] +s2s0[3]) + J0[3] - (kapps[2,2] mtau2/LL2[3]) Exp[I fi] J1[3] + (kapps[2,2] mtau2/LL2[3])^2 Exp[I 2 fi] J2[3] ),{fi,0,+Pi-eeps}] )];
(* tauk22= (1/2/Pi) 2 Re[NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) A1[3,kapps[2,2] mtau2, fi],{fi,0,+Pi-eeps}]]; *)
tauk21=(1/2/Pi) 2 Re [( (1/Pi) NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) tr12d[3,Exp[w]]/(1 + (kapps[2,1] mtau2/LL2[3]) Exp[I fi] Exp[-w]),{fi,0,+Pi-eeps},{w,lns0[3],lnsm},MaxRecursion->30,WorkingPrecision->16] +  NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) ( f1ss0[3]/( (kapps[2,1] mtau2/LL2[3]) Exp[I fi] +s1s0[3])  + f2ss0[3]/( (kapps[2,1] mtau2/LL2[3]) Exp[I fi] +s2s0[3]) + J0[3] - (kapps[2,1] mtau2/LL2[3]) Exp[I fi] J1[3] + (kapps[2,1] mtau2/LL2[3])^2 Exp[I 2 fi] J2[3] ),{fi,0,+Pi-eeps}] )];
(* tauk21= (1/2/Pi) 2 Re[NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) A1[3,kapps[2,1] mtau2, fi],{fi,0,+Pi-eeps}]]; *)
		 
rtaudBG22=alt[2,1] tauk21 +  alt[2,2] tauk22;
(* rtaudBG22 is rtau(D=0) calculated with dBGA[2/2] approach for the Adler function *)
(* ddBG22[fi_]:=alt[2,1] a4dk21complex[fi] + alt[2,2]  a4dk22complex[fi]; *)

tauk11=(1/2/Pi) 2 Re [( (1/Pi) NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) tr12d[3,Exp[w]]/(1 + (kapps[1,1] mtau2/LL2[3]) Exp[I fi] Exp[-w]),{fi,0,+Pi-eeps},{w,lns0[3],lnsm},MaxRecursion->30,WorkingPrecision->16] +  NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) ( f1ss0[3]/( (kapps[1,1] mtau2/LL2[3]) Exp[I fi] +s1s0[3])  + f2ss0[3]/( (kapps[1,1] mtau2/LL2[3]) Exp[I fi] +s2s0[3]) + J0[3] - (kapps[1,1] mtau2/LL2[3]) Exp[I fi] J1[3] + (kapps[1,1] mtau2/LL2[3])^2 Exp[I 2 fi] J2[3] ),{fi,0,+Pi-eeps}] )];
(* tauk11= (1/2/Pi) 2 Re[NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) a4dk11complex[fi],{fi,0,+Pi-eeps}]]; *)

rtaudBG11= tauk11;
 (*--------------------------------------------*)
 (*--------------------------------------------*)
(* Now the TPS evaluation (with {tilde A}_n) of rtau(D=0): *)		 

		 rtauA1=(1/2/Pi) 2 Re [( (1/Pi) NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) tr12d[3,Exp[w]]/(1 + ( mtau2/LL2[3]) Exp[I fi] Exp[-w]),{fi,0,+Pi-eeps},{w,lns0[3],lnsm},MaxRecursion->30,WorkingPrecision->16] +  NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) ( f1ss0[3]/( ( mtau2/LL2[3]) Exp[I fi] +s1s0[3])  + f2ss0[3]/( ( mtau2/LL2[3]) Exp[I fi] +s2s0[3]) + J0[3] - ( mtau2/LL2[3]) Exp[I fi] J1[3] + ( mtau2/LL2[3])^2 Exp[I 2 fi] J2[3] ),{fi,0,+Pi-eeps},MaxRecursion->30,WorkingPrecision->16] )];
(* rtauA1=(1/2/Pi) 2 Re[ NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) A1[3,mtau2,fi],{fi,0,+Pi-eeps}]]; *)

		 (* rtautA2alt= (1/2/Pi) 2  (1/be0[3]/Pi) Re [  ( NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi])  ((mtau2/LL2[3]) Exp[I fi]) Exp[-w] tr12d[3,Exp[w]]/(1 + (mtau2/LL2[3]) Exp[I fi] Exp[-w])^2,{fi,0,+Pi-eeps},{w,lns0[3],lnsm},MaxRecursion->30,WorkingPrecision->16]  +  NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) ((mtau2/LL2[3]) Exp[I fi]) ((f1ss0[3]*Pi)/((mtau2/LL2[3]) Exp[I fi]+s1s0[3])^2 + (f2ss0[3]*Pi)/((mtau2/LL2[3]) Exp[I fi]+s2s0[3])^2 + J1[3] - 2  (mtau2/LL2[3]) Exp[I fi] J2[3] + 3 (mtau2/LL2[3])^2 Exp[I 2 fi] J3[3] ),{fi,0,+Pi-eeps},MaxRecursion->30,WorkingPrecision->16] ) ]; *)
		 rtautA2= (2/be0[3]) (1/2/Pi) 2 Re [( (1/Pi) NIntegrate[Exp[I fi] (1 + Exp[I fi])^2 (1 - 2 Exp[I fi]) tr12d[3,Exp[w]]/(1 + ( mtau2/LL2[3]) Exp[I fi] Exp[-w]),{fi,0,+Pi-eeps},{w,lns0[3],lnsm},MaxRecursion->30,WorkingPrecision->16] +  NIntegrate[Exp[I fi] (1 + Exp[I fi])^2 (1 - 2 Exp[I fi]) ( f1ss0[3]/( ( mtau2/LL2[3]) Exp[I fi] +s1s0[3])  + f2ss0[3]/( ( mtau2/LL2[3]) Exp[I fi] +s2s0[3]) + J0[3] - ( mtau2/LL2[3]) Exp[I fi] J1[3] + ( mtau2/LL2[3])^2 Exp[I 2 fi] J2[3] ),{fi,0,+Pi-eeps},MaxRecursion->30,WorkingPrecision->16] )];
(* rtautA2= (1/2/Pi) 2 Re[NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) tA2[3,mtau2,fi],{fi,0,+Pi-eeps}]], with integration by parts (pp) applied once; *)

		 (* rtautA3alt= (1/2/Pi) (1/2/be0[3]^2/Pi) 2 Re[ ( NIntegrate[ (1 + Exp[I fi])^3 (1 - Exp[I fi]) ((mtau2/LL2[3]) Exp[I fi]) Exp[-w] tr12d[3,Exp[w]]/(1 + (mtau2/LL2[3]) Exp[I fi] Exp[-w])^2,{fi,0,+Pi-eeps},{w,lns0[3],lnsm},MaxRecursion->30,WorkingPrecision->16] -2  NIntegrate[ (1 + Exp[I fi])^3 (1 - Exp[I fi]) ((mtau2/LL2[3]) Exp[I fi]) Exp[- w] tr12d[3,Exp[w]]/(1 + (mtau2/LL2[3]) Exp[I fi] Exp[-w])^3,{fi,0,+Pi-eeps},{w,lns0[3],lnsm},MaxRecursion->30,WorkingPrecision->16] +   NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) ((mtau2/LL2[3]) Exp[I fi]) ( (f1ss0[3]*Pi) (-s1s0[3] + (mtau2/LL2[3]) Exp[I fi])/(s1s0[3] + (mtau2/LL2[3]) Exp[I fi])^3  + (f2ss0[3]*Pi) (-s2s0[3] + (mtau2/LL2[3]) Exp[I fi])/(s2s0[3] + (mtau2/LL2[3]) Exp[I fi])^3 ),{fi,0,+Pi-eeps},MaxRecursion->30,WorkingPrecision->16] )]; *)
		 rtautA3= (1/be0[3]^2) (1/2/Pi) 2 Re [( (1/Pi) NIntegrate[Exp[I fi] (1 + Exp[I fi]) (1 - Exp[I fi] - 8 Exp[2 I fi]) tr12d[3,Exp[w]]/(1 + ( mtau2/LL2[3]) Exp[I fi] Exp[-w]),{fi,0,+Pi-eeps},{w,lns0[3],lnsm},MaxRecursion->30,WorkingPrecision->16] +  NIntegrate[Exp[I fi] (1 + Exp[I fi]) (1 -  Exp[I fi] - 8 Exp[2 I fi]) ( f1ss0[3]/( ( mtau2/LL2[3]) Exp[I fi] +s1s0[3])  + f2ss0[3]/( ( mtau2/LL2[3]) Exp[I fi] +s2s0[3]) + J0[3] - ( mtau2/LL2[3]) Exp[I fi] J1[3] + ( mtau2/LL2[3])^2 Exp[I 2 fi] J2[3] ),{fi,0,+Pi-eeps},MaxRecursion->30,WorkingPrecision->16] )];	
(* rtautA3= (1/2/Pi) 2 Re[NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) tA3[3,mtau2,fi],{fi,0,+Pi-eeps}]], with integration by parts (pp) applied twice; *)

		 (* rtautA4alt=(1/2/Pi) 2 (1/6/be0[3]^3/Pi) Re[ ( NIntegrate[  (1 + Exp[I fi])^3 (1 - Exp[I fi]) ((mtau2/LL2[3]) Exp[I fi]) Exp[-w] (1 - 4 Exp[-w] (mtau2/LL2[3]) Exp[I fi] + Exp[-2 w] (mtau2/LL2[3])^2 Exp[2 I fi]) tr12d[3,Exp[w]]/(1 + (mtau2/LL2[3]) Exp[I fi] Exp[-w])^4,{fi,0,+Pi-eeps},{w,lns0[3],lnsm},MaxRecursion->30,WorkingPrecision->16] + NIntegrate[ (1 + Exp[I fi])^3 (1 - Exp[I fi]) ((mtau2/LL2[3]) Exp[I fi]) ( (f1ss0[3]*Pi) (s1s0[3]^2 - 4 s1s0[3] (mtau2/LL2[3]) Exp[I fi] + (mtau2/LL2[3])^2 Exp[2 I fi] )/(s1s0[3] + (mtau2/LL2[3]) Exp[I fi])^4  + (f2ss0[3]*Pi) (s2s0[3]^2 - 4 s2s0[3] (mtau2/LL2[3]) Exp[I fi] + (mtau2/LL2[3])^2 Exp[2 I fi] )/(s2s0[3] + (mtau2/LL2[3]) Exp[I fi])^4 ),{fi,0,+Pi-eeps}] )]; *)
		 rtautA4=(1/3/be0[3]^3) (1/2/Pi) 2 Re [( (1/Pi) NIntegrate[Exp[I fi] (1 - 27 Exp[2 I fi] - 32 Exp[3 I fi]) tr12d[3,Exp[w]]/(1 + ( mtau2/LL2[3]) Exp[I fi] Exp[-w]),{fi,0,+Pi-eeps},{w,lns0[3],lnsm},MaxRecursion->30,WorkingPrecision->16] +  NIntegrate[Exp[I fi] (1 - 27 Exp[2 I fi] - 32 Exp[3 I fi])  ( f1ss0[3]/( ( mtau2/LL2[3]) Exp[I fi] +s1s0[3])  + f2ss0[3]/( ( mtau2/LL2[3]) Exp[I fi] +s2s0[3]) + J0[3] - ( mtau2/LL2[3]) Exp[I fi] J1[3] + ( mtau2/LL2[3])^2 Exp[I 2 fi] J2[3] ),{fi,0,+Pi-eeps},MaxRecursion->30,WorkingPrecision->16] )];
		 (* rtautA4=(1/2/Pi) 2 Re[NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) tA4[3,mtau2,fi],{fi,0,+Pi-eeps},MaxRecursion->30,WorkingPrecision->16]]with integration by parts (pp) applied three times; *)

rtauTPS=rtauA1+td1 rtautA2+td2 rtautA3+td3 rtautA4;

(*------------------------------------------------------*)
(*------------------------------------------------------*)
(* Now follows calculation of the leading-beta(LB) and beyond LB (bLB) rtau; these calculations appear to be less predictive than rtaudBG22: *)
(* - using the contour integration of the corresp. dv-terms : *)
	       (* we will use RScl |mu^2| = Q^2 Exp[CC] *)

CCbMS = -5/3;
eCbMS= Exp[CCbMS];(*CCbMS=-5/3*)
<<tFtau.m;
  (*---------*)

mtau2=1.77686^2; 

 kap=LL2[3]/mtau2/eCbMS;
 tFtaux[xp_]:=tFtau[kap Exp[xp]];

  (* tr1x[nq_,xp_]:=tr12d[nq,Exp[xp]]; *)
  tr1x[nq_,xp_]:=tr12d[nq,Exp[xp]]; (* for xp>x0=Log[s0]; otherwise we will set equal tr1x[xp]=0 *)
 (* x0=Log[s0]; *)
   rtauLB=(1/Pi) NIntegrate[tFtaux[xp] tr1x[3,xp],{xp,lns0[3],10.^3},MaxRecursion->30,WorkingPrecision->25]+(1/Pi) NIntegrate[tFtaux[xp] tr1x[3,xp],{xp,10.^3,10.^10},MaxRecursion->30,WorkingPrecision->20]+ f1ss0[3] tFtau[kap s1s0[3]]/s1s0[3] +  f2ss0[3] tFtau[kap s2s0[3]]/s2s0[3];
rtauLBalt=(1/Pi) NIntegrate[tFtaux[xp] tr1x[3,xp],{xp,lns0[3],Infinity},MaxRecursion->30,WorkingPrecision->20]+ f1ss0[3] tFtau[kap s1s0[3]]/s1s0[3] +  f2ss0[3] tFtau[kap s2s0[3]]/s2s0[3]; 
		 
		 (* It was IMPORTANT to include here MaxRecursion->30,WorkingPrecision->25, rtauLB changes dramatically (from =appr 0.17 to =appr 0.20); this change made on Sept.5, 2018; this problem seems to occur only in Mathematica 11 (11.3); in earlier versions of Mathematica (9,10), it seems that rtauLB was calculated correctly without specifying (high) WorkingPrecision; this is furstrating *)
		 


		 xxt=mtau2 Exp[CC]; 
		    (*     eps=10^-10; *)
		    eps=10.^-4;
	       mxpts=1000;
     (* theta_max=Pi-eps taken *)
eCC=Exp[CC];


(*   tau1=2 (1/2/Pi) Re[ NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) A1[3, eCC mtau2,fi],{fi,0.,+Pi-eps}]]; *)

		    (* Now approximate contour integration, taking only NN points on the semicircle in the Q^2-plane and interpolating over them: *)
		    (* 'v3ap' *)
NN=30;
    ttabA1=Table[{kk Pi/NN, (1 + Exp[I kk Pi/NN])^3 (1 - Exp[I kk Pi/NN]) A1[3,eCC (mtau2),kk Pi/NN]},{kk,0,NN-1}];
    ttabtA2=Table[{kk Pi/NN, (1 + Exp[I kk Pi/NN])^3 (1 - Exp[I kk Pi/NN]) tA2[3,eCC (mtau2),kk Pi/NN]},{kk,0,NN-1}];
    ttabtA3=Table[{kk Pi/NN, (1 + Exp[I kk Pi/NN])^3 (1 - Exp[I kk Pi/NN]) tA3[3,eCC (mtau2),kk Pi/NN]},{kk,0,NN-1}];
    ttabtA4=Table[{kk Pi/NN, (1 + Exp[I kk Pi/NN])^3 (1 - Exp[I kk Pi/NN]) tA4[3,eCC (mtau2),kk Pi/NN]},{kk,0,NN-1}];

 (* Now the interpolated approximants on the semicircle: *)
		    A1con[ffi_]:=Interpolation[ttabA1][ffi];
		    tA2con[ffi_]:=Interpolation[ttabtA2][ffi];
		    tA3con[ffi_]:=Interpolation[ttabtA3][ffi];
		    tA4con[ffi_]:=Interpolation[ttabtA4][ffi];
tau1ap=2 (1/2/Pi) Re[ NIntegrate[A1con[fi],{fi,0.,+Pi-eps}]];
tau2ap=2 (1/2/Pi) Re[ NIntegrate[tA2con[fi],{fi,0.,+Pi-eps}]];
tau3ap=2 (1/2/Pi) Re[ NIntegrate[tA3con[fi],{fi,0.,+Pi-eps}]];
tau4ap=2 (1/2/Pi) Re[ NIntegrate[tA4con[fi],{fi,0.,+Pi-eps}]];
rtaubLB1ap=tau2ap T1;
rtaubLB2ap=tau3ap T2;
rtaubLB3ap=tau4ap T3;
rtauLBbLBap=rtauLB + T1 tau2ap+ T2 tau3ap + T3 tau4ap;
rtaubLBap=T1 tau2ap+ T2 tau3ap + T3 tau4ap;
rtauLB
rtaubLB1ap
rtaubLB2ap
rtaubLB3ap
rtauLBbLBap
		 
rtaubLB1=T1 rtautA2
rtaubLB2=T2 rtautA3		 
rtaubLB3=T3 rtautA4
rtauLBbLB=rtauLB + T1 rtautA2+ T2 rtautA3 + T3 rtautA4;		 
rtaubLB=T1 rtautA2+ T2 rtautA3 + T3 rtautA4;
rtauLB
rtaubLB
rtauLBbLB


(*----------------------------------------------------*)
     (* Further comments, March 16, 2012: *)
  (* Note: a2d[nq,Q2/LL2[nq],0] - adis2d[nq,Q2/LL2[nq],0] ~ 10^-9 for all Q2, including Q2 ~ 1-1000 GeV^2; i.e., numerical precision of our calculation is |delta a| ~ |delta A1| ~ 10^-9 *)
   (* On the other hand, diff[Q^2] < ~ sim 10^-9 for Q2=10^3; so, we can trust the value od  diff[Q2] only for Q2 <= 10^2 GeV^2 *) 
   (* However, when s0>>1, we have formally  diff[Q2]~ (L2/Q2)^5; practically this is not so, because the coefficients at (L2/Q2)^n for n=5,6,... are very large (and they grow out of control for n even larger) *) tau1
   (* March 16-20, 2012 *)
(*----------------------------------------------------*)
(*----------------------------------------------------*)
f[x_] := -Exp[x]/(1 - Exp[x]);
F[z_, n_] := -Derivative[n][f][x] /. {x -> Log[z]}
  (* --> Simplify[F[z, n] - PolyLog[-n, z]] gives zero for any n nonnegative integer *)
    (* However, Derivative[n][F][z, n] /. {n -> 2.}, and subsequent %/.z-> 5. or sth like that, does not give any numerical result (even when N[...] applied) *)
    (* Further, Simplify[F[z, 1.1] - PolyLog[-1.1, z]], and replacing z->5.  (even when N[...] applied) - does not give a number *)
    (* -- which would suggest that Mathematica is not able to extend  F[z, n] to noninteger n, as PolyLog[-n,z]; and cannot make derivatives with respect to n of F[z,n] *)
(*----------------------------------------------------*)
(*----------------------------------------------------*)
  (* Now calculation of tA_{nu}(Q^2,nq) for noninteger nu: *)

 (* s = sigma/L2; u = |Q^2|/L2; Arg[Q^2]=fi  *)
  tr1MA1l[nq_,s_]:=(Pi/be0[nq])(1/(Log[s]^2+Pi^2));

  (* In  tAMA1l[nu,u,fi] we take into account that: *)
  (* Limit[1/Log[u]^(nu + 1) - PolyLog[-nu, 1/u]/Gamma[nu + 1], u -> 1] = -Zeta[-nu]/Gamma[nu+1], which holds for all real nu, it seems *)

  (* <<Li__nu.m; *)
  (* Li__nu.m gives a stable version of polylog[nu,u] *)

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_]:=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]];
(* tAMA1l[nq,nu,u,fi]={tilde A}_{nu}^{MA,1-l.}(Q^2= u LL2[nq] Exp[I fi], nq) *) 

  (* NOte, at u =approx +1, we have singularity of the type infty-infty, so around there we replace tAMA1l[nu,u,fi] with Limit[tAMA1l[nu,u,fi], u Exp[I fi] ->1] *)

  (* tA2d[nq,nu,Q2,fi] = {tilde A}_{nu}(Q^2=Q2 Exp[I f];nq), with Q2=|Q^2|; Q^2 = Q2 Exp[I fi] (u=Q2/LL2[nq]) *)
 (* ssig Exp[-I fi]/u = sigma/Q^2 *)



 (* The above integration cannot be applied for nupp=1 (indeterminate expressions) --> there we have to choose, e.g., nupp=1.000000001 *)

  (* Below: alternative tA (tAalt), where nupp can be any integer without danger -- using polylog[nu,x] instead of GenLia[-nu,-x]; (however, tAalt is may not good for making derivatives of tAalt[nu] with respect to nu - check this) *)
  (* It seems to be OK, we will declare now tAal as: tA: *)

		 tA2d[nq_,nupp_,Q2_,fi_]:= If[Abs[nupp]<=2.,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]],10.^9},MaxRecursion->30,WorkingPrecision->16]),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]],10.^6},MaxRecursion->30,WorkingPrecision->16])];
  tA[nq_,nupp_,Q2_,fi_]:=  tA2d[nq,nupp,Q2,fi];
  
  (* tA2d[nq,nu,Q2,fi]= {tilde A}_{nu}(Q^2=Q2 Exp[I fi];nq); Q2=|Q^2| *)
  (* Valid for all -2 < nupp; even for integer nupp *)
 (* Note, in Mathematica 11.3 it is important to include MaxRecursion->30,WorkingPrecision->16 in this integration, otherwise large imprecisions or even error occur in the integration *) 

  
(* Now as the above tA, but with the upper integration limit (wm) adjustable: *)
  (* usually wm >= 10.^5 *)

 tA2dwm[nq_,nupp_, Q2_, fi_, wm_] := 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]], wm},MaxRecursion->30,WorkingPrecision->16]);
  (* tA2dwm[nq,nu,Q2,fi,wm]= {tilde A}_{nu}(Q^2=Q2 Exp[I fi];nq); Q2=|Q^2| *)
  (* Use: 10^5 < wm < 10^10 *)
  (* Note, in Mathematica 11.3 it is important to include MaxRecursion->30,WorkingPrecision->16 in this integration, otherwise large imprecisions or even error occur in the integration *)
  (*---------------------*)
  
 (* Now the spacelike couplings without tilde, with precision up to O(tA_{nupp+4}), beyond which we truncate: *)
(* I could assume that 0 =< nupp < 1, or 1 =< nupp <2, or ... *)

 (* A special case in the formulas below is: nupp=0, i.e., when the indices are integers: AN[0,u,fi] = A[N,u,fi]=A_N[u,fi] *)
 
 (* --------------------------------------------------_*)
 (* First the case of N^4LO truncation ('N4', i.e., up to FOUR correction terms; recommended application of these N4-coefficients: if we have available N^4LO series of an observable: D(Q2)_{pt} = a^{nu} + d1 a^{nu+1} + d2 a^{nu+2} + d3 a^{nu+3} + d4 a^{nu+4}): *)
 (* At present, no observable with available N^4LO perturbation expansion exists *)

  (* A2dN4[Nf,M,nupp,Q2,fi] = A2d[Nf,M+nupp,Q2,fi] = A_{nupp+M}(Q^2; Nf) in N^4LO truncation approach; where Q^2 = Q2 Exp[i fi]; Q2=|Q^2|(= u LL2[Nf]) in GeV^2 *)

  (* In practice, we take M=0,1,2,3,4; and usually 0<nupp<2 (then n=M+nupp noninteger) *)

A2dN4[nq_,4,nupp_,Q2_,fi_]:=tA2d[nq,4+nupp,Q2,fi];
A2dN4[nq_,3,nupp_,Q2_,fi_]:=tA2d[nq,3+nupp,Q2,fi] + tk1f[nq,3+nupp] tA2d[nq,4+nupp,Q2,fi];
A2dN4[nq_,2,nupp_,Q2_,fi_]:=tA2d[nq,2+nupp,Q2,fi] + tk1f[nq,2+nupp] tA2d[nq,3+nupp,Q2,fi] + tk2f[nq,2+nupp] tA2d[nq,4+nupp,Q2,fi];
A2dN4[nq_,1,nupp_,Q2_,fi_]:=tA2d[nq,1+nupp,Q2,fi] + tk1f[nq,1+nupp] tA2d[nq,2+nupp,Q2,fi] + tk2f[nq,1+nupp] tA2d[nq,3+nupp,Q2,fi] + tk3f[nq,1+nupp] tA2d[nq,4+nupp,Q2,fi];
A2dN4[nq_,0,nupp_,Q2_,fi_]:=tA2d[nq,nupp,Q2,fi] + tk1f[nq,nupp] tA2d[nq,1+nupp,Q2,fi] + tk2f[nq,nupp]  tA2d[nq,2+nupp,Q2,fi]+ tk3f[nq,nupp] tA2d[nq,3+nupp,Q2,fi] + tk4f[nq,nupp] tA2d[nq,4+nupp,Q2,fi];

  (* Shorter notation: *)
  AN4[nq_, M_,nupp_,Q2_,fi_]:=A2dN4[nq,M,nupp,Q2,fi];

  (* For integer index nupp=1, at N^4 LO ('N4') order: *)

A5N4[nq_,Q2_,fi_]:=tA2d[nq,5,Q2,fi];  (* This is the analog of a(Q^2)^5 in N^4LO approximation *)
A4N4[nq_,Q2_,fi_]:=tA4[nq,Q2,fi]+ tk1f[nq,4] tA2d[nq,5,Q2,fi];
   (* This is the analog of a(Q^2)^4 in N^4 LO *)
A3N4[nq_,Q2_,fi_]:=tA3[nq,Q2,fi]+ tk1f[nq,3] tA4[nq,Q2,fi]+ tk2f[nq,3] tA2d[nq,5,Q2,fi];
   (* This is the analog of a(Q^2)^3 in N^4 LO *)
A2N4[nq_,Q2_,fi_]:=tA2[nq,Q2,fi]+ tk1f[nq,2] tA3[nq,Q2,fi]+ tk2f[nq,2] tA4[nq,Q2,fi] + tk3f[nq,2] tA2d[nq,5,Q2,fi];
  (* This is the analog of a(Q^2)^2 in N^4 LO *)
		 A1N4[nq_,Q2_,fi_]:=A1[nq,Q2,fi];
  (*=AN4[0,1,Q2,fi]; this is the analog of a(Q^2) *)
		 (* tkjf[n] -> tkjf[nq,n], on Nov.25, 2018 *) 
(*-----------------------------------------------------*)
 (* Now the case of N^3LO truncation ('N3', i.e., up to THREE correction terms; recommended application of these N3-coefficients: if we have available N^3LO series of an observable: D(Q2)_{pt} = a^{nu} + d1 a^{nu+1} + d2 a^{nu+2} + d3 a^{nu+3}): *)
 (* This is, for example, the case of the presently available Adler function perturbation expansion *)

  (* A2dN3[Nf,M,nupp,Q2,fi] = A[Nf,M+nupp,Q2,fi] = A_{nupp+M}(Q^2;Nf) in N^3LO truncation approach; where  Q^2 = Q2 Exp[i fi]; Q2=|Q^2|(= u LL2[Nf]) in GeV^2 *)
  (* In practice, we take M=0,1,2,3; and usually 0<nupp<2 (then n=M+nupp noninteger) *)

A2dN3[nq_,3,nupp_,Q2_,fi_]:=tA2d[nq,3+nupp,Q2,fi];
A2dN3[nq_,2,nupp_,Q2_,fi_]:=tA2d[nq,2+nupp,Q2,fi] + tk1f[nq,2+nupp] tA2d[nq,3+nupp,Q2,fi];
A2dN3[nq_,1,nupp_,Q2_,fi_]:=tA2d[nq,1+nupp,Q2,fi] + tk1f[nq,1+nupp] tA2d[nq,2+nupp,Q2,fi] + tk2f[nq,1+nupp] tA2d[nq,3+nupp,Q2,fi];
A2dN3[nq_,0,nupp_,Q2_,fi_]:=tA2d[nq,nupp,Q2,fi] + tk1f[nq,nupp] tA2d[nq,1+nupp,Q2,fi] + tk2f[nq,nupp]  tA2d[nq,2+nupp,Q2,fi]+ tk3f[nq,nupp] tA2d[nq,3+nupp,Q2,fi];

  (* Shorter notation: *)
  AN3[nq_, M_,nupp_,Q2_,fi_]:=A2dN3[nq,M,nupp,Q2,fi];

  (* For integer index nupp=1, at N^3 LO ('N3') order: *)
   A4N3[nq_,Q2_,fi_]:=tA4[nq,Q2,fi];  (* This is the analog of a(Q^2)^4 *)
		 A3N3[nq_,Q2_,fi_]:=tA3[nq,Q2,fi]+ tk1f[nq,3] tA4[nq,Q2,fi];
   (* This is the analog of a(Q^2)^3 *)
		 A2N3[nq_,Q2_,fi_]:=tA2[nq,Q2,fi]+ tk1f[nq,2] tA3[nq,Q2,fi]+ tk2f[nq,2] tA4[nq,Q2,fi];
   (* This is the analog of a(Q^2)^2 *)
   A1N3[nq_,Q2_,fi_]:=A1[nq,Q2,fi];
  (*=AN3[0,1,Q2,fi]; this is the analog of a(Q^2) *)
		 (* tkjf[n] -> tkjf[nq,n], on Nov.25, 2018 *) 
		 (*-----------------------------------------------------*)
   (* Now the case of N^2LO truncation ('N2', i.e., up to TWO correction terms; recommended application of these N2-coefficients: if we have available N^2LO series of an observable: D(Q2)_{pt} = a^{nu} + d1 a^{nu+1} + d2 a^{nu+2}): *)
 (* This is, for example, the case of presently available DIS structure function perturbation expansion *) 

  (* A2dN2[Nf,M,nupp,Q2,fi] = A[Nf,M+nupp,Q2,fi] = A_{nupp+M}(Q^2) in N^2LO truncation approach; where  Q^2 = Q2 Exp[i fi]; Q2=|Q^2| (= u LL2[Nf]) in GeV^2 *)
  (* In practice, we take M=0,1,2; and usually 0<nupp<2 (then n=M+nupp noninteger) *)

A2dN2[nq_,2,nupp_,Q2_,fi_]:=tA2d[nq,2+nupp,Q2,fi];
A2dN2[nq_,1,nupp_,Q2_,fi_]:=tA2d[nq,1+nupp,Q2,fi] + tk1f[nq,1+nupp] tA2d[nq,2+nupp,Q2,fi];
A2dN2[nq_,0,nupp_,Q2_,fi_]:=tA2d[nq,nupp,Q2,fi] + tk1f[nq,nupp] tA2d[nq,1+nupp,Q2,fi] + tk2f[nq,nupp]  tA2d[nq,2+nupp,Q2,fi];

 (* Shorter notation: *)
  AN2[nq_, M_,nupp_,Q2_,fi_]:=A2dN2[nq,M,nupp,Q2,fi];  

  (* For integer index nupp=1, at N^2 LO ('N2') order: *)
  A3N2[nq_,Q2_,fi_]:=tA3[nq,Q2,fi];  (* This is the analog of a(Q^2)^3 *)
		 A2N2[nq_,Q2_,fi_]:=tA2[nq,Q2,fi]+tk1f[nq,2] tA3[nq,Q2,fi];
   (* This is the analog of a(Q^2)^2 *)
		 (* tk1f[2] -> tk1f[nq,2], on Nov.25, 2018 *) 
   A1N2[nq_,Q2_,fi_]:=A1[nq,Q2,fi]; (*=AN2[0,1,Q2,fi]; this is the analog of a(Q^2) *)
  (*-----------------------------------------------------------------*)
 (*-----------------------------------------------------*)
 (* Now the case of N^1LO truncation ('N1', i.e., up to ONE correction term; recommended application of these N1-coefficients: if we have available N^1LO series of an observable: D(Q2)_{pt} = a^{nu} + d1 a^{nu+1}: *)
 (* This is, for example, the case of presently available DIS structure function perturbation expansion *) 

  (* A2dN1[Nf,M,nupp,Q2,fi] = A[Nf,M+nupp,Q2,fi] = A_{nupp+M}(Q^2) in N^1LO truncation approach; where  Q^2 = Q2 Exp[i fi]; Q2=|Q^2| (= u LL2[Nf]) in GeV^2 *)
  (* In practice, we take M=0,1; and usually 0<nupp<2 (then n=M+nupp noninteger) *)

  A2dN1[nq_,1,nupp_,Q2_,fi_]:=tA2d[nq,1+nupp,Q2,fi];
  A2dN1[nq_,0,nupp_,Q2_,fi_]:=tA2d[nq,nupp,Q2,fi] + tk1f[nq,nupp] tA2d[nq,1+nupp,Q2,fi];

   (* Shorter notation: *)
  AN1[nq_, M_,nupp_,Q2_,fi_]:=A2dN1[nq,M,nupp,Q2,fi];  

(* For integer index nupp=1, at N LO ('N1') order: *)
  A2N1[nq_,Q2_,fi_]:=tA2[nq,Q2,fi];  (* This is the analog of a(Q^2)^2 *)
  A1N1[nq_,Q2_,fi_]:=A1[nq,Q2,fi]; (*=AN1[0,1,Q2,fi]; this is the analog of a(Q^2) *)
  (*-----------------------------------------------------------------*)


 (* Note that the difference {tilde A}_n(Q^2)-{tilde a_pt}_n(Q^2) is for ALL n (even n large) roughly the same (~ (L2/Q^2)^5 ); which means that the "man" series of an observable D(Q^2) in {tilde A}_n and the "mpt" series in {tilde a]_n for any observable will differ from each other, more and more when the number of terms increases (because the coefficients {tilde k}_n (at {tilde a_{pt}}_n) increase rapidly (as ~ n!) when n increases: D_{man}^{[N]}-D{mpt}^{[N]} increases (rapidly) when the truncation number N increases *)

  (* This may also explain why the truncated "power" series D_{an}^{[N]} - D_{pt}^{[N]} also increases, even more rapidly, when the truncation number N increases -- note that, due to our truncated construction of A_{nu}'s on the basis of tA_{nu+m}'s, D_{man}^{[N]}=D_{an}^{[N]}, and this may make the difference D_{an}^{[N]}-D_{pt}^{[N]} = D_{man}^{[N]}-D_{pt}^{[N]} sometimes even greater than D_{man}^{[N]}-D{mpt}^{[N]} *)
  (*------------------------------*)
  (* Save["2danQCDcouplrt0201al01189.m", {s0, s1s0, s2s0, f1ss0, f2ss0, M0, LLamb, LL2, be0, c1, c22d, c32d,c42d, c52d, c2MS,c3MS, btP11, betaMS4l,kappa2d,  alphaMSMZ, apert, aMSnf3, A1, tA2, tA3, tA4, tA, AN4, AN3, A5N4, A4N4, A3N4, A2N4, A1N4, A4N3, A3N3, A2N3, A1N3, AN2, A3N2, A2N2, A1N2, AN1, A2N1, A1N1,alphaMSMZ,rtaudBG22,rtauTPS,rtauLBbLB,d1bar,d2bar,d3bar,td1bar,td2bar,td3bar,d1,d2,d3,td1,td2,td3,td1bLB,td2bLB,td3bLB}] *)
 (* Save["2danQCDcouplrt0201al01185.m", {s0, s1s0, s2s0, f1ss0, f2ss0, M0, LLamb, LL2, be0, c1, c22d, c32d,c42d, c52d, c2MS,c3MS, btP11, betaMS4l,kappa2d,  alphaMSMZ, apert, aMSnf3, A1, tA2, tA3, tA4, tA, AN4, AN3, A5N4, A4N4, A3N4, A2N4, A1N4, A4N3, A3N3, A2N3, A1N3, AN2, A3N2, A2N2, A1N2, AN1, A2N1, A1N1,alphaMSMZ,rtaudBG22,rtauTPS,rtauLBbLB,d1bar,d2bar,d3bar,td1bar,td2bar,td3bar,d1,d2,d3,td1,td2,td3,td1bLB,td2bLB,td3bLB}] *)
 (* Save["2danQCDcouplrt0201al01181.m", {s0, s1s0, s2s0, f1ss0, f2ss0, M0, LLamb, LL2, be0, c1, c22d, c32d,c42d, c52d, c2MS,c3MS, btP11, betaMS4l,kappa2d,  alphaMSMZ, apert, aMSnf3, A1, tA2, tA3, tA4, tA, AN4, AN3, A5N4, A4N4, A3N4, A2N4, A1N4, A4N3, A3N3, A2N3, A1N3, AN2, A3N2, A2N2, A1N2, AN1, A2N1, A1N1,alphaMSMZ,rtaudBG22,rtauTPS,rtauLBbLB,d1bar,d2bar,d3bar,td1bar,td2bar,td3bar,d1,d2,d3,td1,td2,td3,td1bLB,td2bLB,td3bLB}] *)

  
