(* ***

Examples:
=========

<<asMOM.m

PiMSml
PiMSsma
PiMSlar
betaMOMml

AsMOMrun[0.1324,91.19,500,2,6]
AsMOMrun[0.1331,91.19,500,3,6]
AsMOMrun[0.1328,91.19,4,2,5]
AsMOMrun[0.1332,91.19,4,3,5]


Notation:
=========

lqq = \log(-q^2/mu^2)
api = alpha_s(mu)/Pi
z   = q^2/4/Mq^2

ca = 3
cf = 4/3
tr = 1/2

z2 = \zeta(2)
z3 = \zeta(3)


Included expressions and Modules:
=================================

betaMOMml: beta function in MOM scheme; all nl quarks are massless

PiMSml:  \Pi(q^2) in MS-bar scheme; massless contribution
PiMSsma: \Pi(q^2) in MS-bar scheme; contribution from heavy quark 
         for q^2 << Mq^2 up to order z^3
PiMSlar: \Pi(q^2) in MS-bar scheme; contribution from heavy quark 
         for q^2 >> Mq^2 up to order 1/z^5

Module: AsMOMrun[alphasMOMmuini, muini, muend, l, nfMOM_]:
  Compute alpha_s^MOM(muend) from alpha_s^MOM(muini)=alphasMOMmuini
  using the l-loop approximation of the beta function.
  nfMOM determines the quark flavour which is considered as massive.
  (nfMOM=6 -> top, nfMOM=5 -> bottom, ...)

*** *)

<<asMOMexpr.m;

numcol = {ca->3,cf->4/3,tr->1/2};
numprec = 20;
cut[x_, n_Integer] := (x^p_Integer /; p > n -> 0);
numcf = {tr->1/2,ca->3,cf->4/3,nl2->nl};
Co[ex__] := Coefficient[ex];
Col[ex__] := Collect[ex];

(* massless beta function: *)

(betaMOMml = -api^2*(1/4*betCSmlMOM[0] * api
		    +1/4^2*betCSmlMOM[1] * api^2
		    +1/4^3*betCSmlMOM[2] * api^3));


(* \Pi(q^2) in MS-bar scheme: *)

(* massless: *)
(PiMSml = (api*Col[-1/4^1*pi1CSml/.numcf,nl,Expand]
	  +api^2*Col[-1/4^2*pi2CSml/.numcf,nl,Expand]
	  +api^3*Col[-1/4^3*pi3CSml/.numcf,{nl,lqq},Expand])
);

(* z->0: *)
ex3ltmp = Expand[-1/4^3*pi3CSs/.numcf];
(PiMSsma = ((api*Co[Expand[    -1/4^1*pi1CSs/.numcf],nh1]/.cut[z,3])
	    +(api^2*Col[Co[Expand[-1/4^2*pi2CSs/.numcf],nh1]/.cut[z,3],z])
	    +(api^3*Col[(Co[ex3ltmp,nh1]+Co[ex3ltmp,nh1,2])/.cut[z,3],{z,nl,Log[x__],lmm}]))
);

(* z->-infinity: *)
ex3ltmp = Expand[-1/4^3*pi3CSl/.numcf];
(PiMSlar = ((api*Co[Expand[    -1/4^1*pi1CSl/.numcf],nh1]/.{z->1/iz}/.cut[iz,5]/.{iz->1/z})
	    +(api^2*Col[Co[Expand[-1/4^2*pi2CSl/.numcf],nh1]/.{z->1/iz}/.cut[iz,5],iz]/.{iz->1/z})
	    +(api^3*Col[(Co[ex3ltmp,nh1]+Co[ex3ltmp,nh1,2])/.{z->1/iz}/.cut[iz,5],
		       {iz,nl,Log[x__],lmm}]/.{iz->1/z}))
);

(* Module: running of alpha_s in the MOM scheme *)

(AsMOMrun[alphasmu0_,mu0_,muend_,loops_,nfMOM_] := Module[
  {rhs,apiMOM,x,bb0MOM,bb1MOM,bb2MOM,nlMOM,massq,massc,massb,masst}, 
  If[(loops<1)||(loops>3), Print["AsMOMrun: PROCEDURE IS NOT IMPLEMENTED FOR ",
				 loops," LOOPS."];
     Return[]];
  masst = 175.;
  massb = 4.70;
  massc = 1.65;
  nlMOM = nfMOM-1;
  If[nfMOM==6, massq=masst];
  If[nfMOM==5, massq=massb];
  If[nfMOM==4, massq=massc];
  bb0MOM[x_] := 1/4*(11/3*ca - 4/3*tr*(nlMOM * 1
				       + b0MOM[-1/4*x/massq^2])
		     )/.numcol;
  bb1MOM[x_] := 1/16*(34*ca^2/3 - tr*(nlMOM * ((20*ca)/3 + 4*cf)
				      + b1MOM[-1/4*x/massq^2])
		      )/.numcol;
  bb2MOM[x_] := 1/64*((betCSmlMOM[2]/.nl->nlMOM) +(betCSMOM2[1/4*x/massq^2,nlMOM])
		      /.{nh->1}/.numcol/.num2 ); 
  If[ mu0 == muend , Return[alphasmu0] ];
  If[ loops==1, rhs = -apiMOM[x]^2*(bb0MOM[x]) ];
  If[ loops==2, rhs = -apiMOM[x]^2*(bb0MOM[x] + bb1MOM[x]*apiMOM[x]) ];
  If[ loops==3, rhs = -apiMOM[x]^2*(bb0MOM[x] + bb1MOM[x]*apiMOM[x]
				    + bb2MOM[x]*apiMOM[x]^2) ];

  sol = N[(Pi*apiMOM[x]/.
	   NDSolve[{apiMOM[mu0^2]== alphasmu0/Pi,
		    x*apiMOM'[x] == rhs}, 
		   apiMOM[x], 
		   {x,muend^2,mu0^2}]
	   /.x->muend^2)[[1]],numprec];
  Return[sol];
]);
