(* RGErundown5l.m: as RGErundown.m, but uses 5-loop RGE (instead of 4-loop) and 4-loop quark threshold matching (instead of 3-loop) *)
(*---------------------------------------------------------------*)
(* As prg. RGErun.m, but ut uses RGE-running downwards, from a given value of alpha(MZ^2,MSbar); *)
  (* in MSbar, 4-loop; with 3-loop matching *)
(*---------------------------------------------------------------*)
  (* -- the program is based on the results of paper Chetyrkin, Kniehl and Steinhauser, PRL79,2184 *)
   (* alphaMSMZ= alpha_s(MZ^2;MSbar) is the EXTERNAL INPUT value; at nf=5; *)
(* Another  EXTERNAL INPUT are the MSbar values mbbar(mbbar)=mb, and mcbar(mcbar)=mc; and kappa: note, the thresholds at kappa*mb and kappa*mc *)
   (* Another input: Q --> at which we want to obtain a(Q^2;nf=3); we assume that Q =< kappa*mc; and kappa=2 *)
(*---------------------------------------------------------------*)
(* Input: Q (in GeV) *)
(* Q: We want to obtain a[Q^2,nf=3]  *)
(* usually: Q=kappa mc *)
(* Q=kappa mc; *)

  (* aa=al/Pi; *)

MZ=91.1876;

   (* mc=1.27; input *)
   (* mb=4.20; input *)

 (* mb=4.20: Note, we use the low value of mb, as suggested by our analysis with Pineda (2013/2014) *)
(* mb=mbbar(mbbar^2);  mc=mcbar(mcbar^2) *)

   (* kappa=2.; the usual input *)

(* Threshold nf=5 -> nf=4 at kappa*mb *)
(* Threshold nf=4 -> nf=3 at kappa*mc *)
  (* Starting Q^2 = MZ^2, in GeV^2 *)

QMZ=MZ;
aaMZ=alphaMSMZ/Pi;
(* aaMZ=a(QMZ^2) *)

   Q5=kappa mb; (* Q4=mu_thr(nf=5->4) *)
   Q4=kappa mc; (* Q3=mu_thr(nf=4->3) *)

   (* x will be 2 Log[Q(nf)/Qnf], wher Q(nf) is the running Q (> Qnf)  *)
   (*-----------------------------------------------------------------*)
   (* First the running in the nf=5 regime from MZ down to Q5=kappa*mb: *)
 bbe0[nq_]:=(11 -2 nq/3)/4;
cc1[nq_]:=(102-38 nq/3)/(11 -2 nq/3)/4; 
bbe2MS[nq_]:=(2857/2 - 5033 nq/18 + 325 nq^2/54)/64;

be3MS[nq_]:=(1/4^4)( (149753/6) + 3564 Zeta[3] - nq ( (1078361/162) + (6508/27) Zeta[3] )  + nq^2 ( (50065/162) + (6472/81) Zeta[3] )  + nq^3 (1093/729) );
bbe3MS[nq_]:=(29243.0 - 6946.30 nq + 405.089 nq^2 + 1.49931 nq^3)/256; 

bbe4MS[nq_]:=(1/4^5) ( (8157455/16) + (621885/2) Zeta[3] - (88209/2) Zeta[4] - 288090 Zeta[5] + nq ( -(336460813/1944) - (4811164/81) Zeta[3] + (33935/6) Zeta[4] + (1358995/27) Zeta[5] ) + nq^2 ( (25960913/1944) + (698531/81) Zeta[3] - (10526/9) Zeta[4] - (381760/81) Zeta[5] ) + nq^3 (-(630559/5832) - (48722/243) Zeta[3] + (1618/27) Zeta[4] + (460/9) Zeta[5] ) + nq^4 ( (1205/2916) - (152/81) Zeta[3] ) );
(* see 1606.08659, Baikov, Chetyrkin and Kuhn *)

cc2MS[nq_]:=bbe2MS[nq]/bbe0[nq]; cc3MS[nq_]:=bbe3MS[nq]/bbe0[nq];
cc4MS[nq_]:=bbe4MS[nq]/bbe0[nq]; 

nff=5;
Clear[Bta];
Bta[z_]:=- bbe0[nff] z^2 (1 + cc1[nff] z + cc2MS[nff] z^2 + cc3MS[nff] z^3 + cc4MS[nff] z^4);

xx5=2*Log[QMZ/Q5];
del=0.01;
sol1a = NDSolve[{y'[x]==Bta[y[x]], y[xx5]==aaMZ}, y, {x,-del,+del}];
aa=First[y[0.]/.sol1a];
aaQ5plus=aa;
(* aa ' is now a(Q5+0)=a(kappa*mb+0; nff=5) *)
  (*-----------------------------------------------------------------*)
   (* Now transition from  a(Q5+0) to  a(Q5-0)=a(kappa*mb-0; nf=5-1): *)
<<RGE4lmatchdown.m;
aaQ5low=aa;
(* aaQ5low = a(Q5-0;nf=5-1) *)
   (*-----------------------------------------------------------------*)
  (* Now the running in the nf=4 regime from Q5=kappa*mb to Q4=kappa*mc: *)

nff=4;
Clear[Bta];
Bta[z_]:=- bbe0[nff] z^2 (1 + cc1[nff] z + cc2MS[nff] z^2 + cc3MS[nff] z^3 + cc4MS[nff] z^4);


xx4=2*Log[Q5/Q4];
del=0.01;
sol1a = NDSolve[{y'[x]==Bta[y[x]], y[xx4]==aaQ5low}, y, {x,-del,+del}];
aa=First[y[0.]/.sol1a];
(* aa ' is now a(Q4+0)=a(kappa*mc+0; nf=4) *)
aaQ4plus=aa;
  (*-----------------------------------------------------------------*)
   (* Now transition from  a(Q4+0) to  a(Q4-0)=a(kappa*mc-0; nf=4-1): *)
<<RGE4lmatchdown.m;
aaQ4low=aa;
(* aaQ4low = a(Q4-0;nf=4-1) *)
  (*-----------------------------------------------------------------*)
   (* 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_]:=- bbe0[nff] z^2 (1 + cc1[nff] z + cc2MS[nff] z^2 + cc3MS[nff] z^3 + cc4MS[nff] z^4);


  xx3=2*Log[Q/Q4];  (* 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 *)
