(* 1.) Type the INPUT value of alpha_s(mtau^2;MSbar) (=alsin): alsin=0.... (for example: alsin=0.320); *)
(* 2.) Call this program: <</home/gorazd/3danQCD.dir/aMSQ2complS_almtauinput.m; *)
(* Now this gives: aMSpert[Q2,phi]=a(Q^2=Q2 Exp[I phi]; MSbar; nf=3)  (Q2=|Q^2|>0) *)

(* A practical way to make a function out of this program: *)
(* aMSp[Q2f_, fif_, alsinf_] := Module[{Q2ff = Q2f, fiff = fif, alsinff = alsinf}, Clear[alsin] ; alsin = alsinff; <</home/gorazd/3danQCD.dir/aMSQ2complS_almtauinput.m; aMSpert[Q2ff, fiff]]; *)
(* aMSp[0.8, 0.4, 0.320] is then a(Q^2=0.8 Exp[I 0.4]; MSbar) with initial condition alpha_s(mtau^2;MSbar)=0.320  (all with nf=3) *)
(*----------------------------------------------------------_*)
(* aMSQ2complS_almtauinput.m: as  aMSQ2complexStable.m, but independent of outher programs for input; the input is a(mtau^2)=amutau; it gives aMSpert[Q2,phi]=a(Q^2=Q2 Exp[I phi]; MSbar; nf=3)  (Q2=|Q^2|>0) *)

(* aMSQ2complexStable.m: as aMSQ2complex.m, but for Q^2 such that |Q^2|<Q2_{br.point} (=appr 0.403628 GeV^2 when alpha_s(MZ^2)=0.1181) we have STABLE results because we integrate along phi=Arg[Q^2] with initial condition NOT at phi0=0 (which would be on the landau cut for |Q^2|<Q2_{br.point}), but at phi0=1 (and integrate then up to phi=+Pi and down to phi=+0); for phi<0 the results are cc of those at phi>0 *)
  (*----------------------------------------------------*)
(* The program calculates a(Q^2) in MSbar scheme, for nf=3; uses 5-loop beta function *)
(*----------------------------------------------------*)
(* INPUT: alsin=alpha_s(mtau^2;MSbar): *)
 (* alsin=0.320; *)
(* end of INPUT *)
(*--------------*)
 mu2in=1.77686^2;   (* mu2in=(mu_in)^2=mtau^2 *)
 ain=alsin/Pi;
(*--------------*)

be0 = 9/4; ce1 = 16/9; ce2MS = 3863/864; ce3MS = 20.990247170138886;  
ce4MS = (8157455/16 - (9801*Pi^4)/20 + 81*(1205/2916 - (152*Zeta[3])/81) + (621885*Zeta[3])/2 + 9*(25960913/1944 - (5263*Pi^4)/405 + (698531*Zeta[3])/81 - (381760*Zeta[5])/81) - 288090*Zeta[5] + 27*(-630559/5832 + (809*Pi^4)/1215 - (48722*Zeta[3])/243 + (460*Zeta[5])/9) + 3*(-336460813/1944 + (6787*Pi^4)/108 - (4811164*Zeta[3])/81 + (1358995*Zeta[5])/27))/2304;
Nf = 3;
  

  (* This gives us be0, ce1, ce2MS, ce3MS, ce4MS, for nf=3; *)
  (* we will need aMSnf3[mu2in] for mu2in=mtau^2; as the initial condition (mu2in=(mtau)^2 = 1.77686^2 GeV^2) *)
 
beta[x_]:=-be0 x^2 (1 + ce1 x + ce2MS x^2 + ce3MS x^3 + ce4MS x^4);
(* ?aMSnf3 *)
(*--------------------------------*)
(*--------------------------------*)
reimbt[u,v]=ComplexExpand[{Re[beta[u + I v]], Im[beta[u + I v]]}];
rebtt=reimbt[u,v][[1]];
imbtt=reimbt[u,v][[2]];
rebt[u0_,v0_]:=rebtt/.{u->u0,v->v0};
imbt[u0_,v0_]:=imbtt/.{u->u0,v->v0};
(*-----------------------------------*)
   (* Now solving the RGE: *)
bt[x_]:=beta[x];

Needs["NumericalCalculus`"]

(* ain=a[mtau^2;MSbar]; *) 
(* Now, at |Q^2|=mu2in=(mu_in)^2 (= mtau^2), Q^2=mu2in Exp[I phi], we integrate along the angle phi from phi=0 to phi0 (=1): *)
phi0=1.;
aabephi=NDSolve[{(-I) yy'[phi]==beta[yy[phi]],yy[0]==ain},yy,{phi,0,phi0+0.1}];
a1mu2inphi[phi_]:=yy[phi]/.aabephi[[1]]; 
(* ' a1mu2inphi[phi]= a[Q^2) with Q^2=mu2in Exp[I phi] *)
ainphi0=yy[phi0]/.aabephi[[1]]; 
  (* ainphi0 = a(mu2in Exp[I phi0]); this will be our inital value for integration along the axis parallel to the real axis in Q^2 plane: a(|Q^2| Exp[I phi0]); these values will be used as initial values for integration along phi, for a(|Q^2| Exp[I phi])  *)

(*---------------------------------------------------*)
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
(*----------------------------------------------*)
   (* Now the following notations: *)
   (* Note: a = u + I v; u=u(x,y), v=v(x,y) where Log[Q^2/mu_in^2]=x + I y, i.e., Q^2 = mu_in^2 Exp[x] Exp[I y]; note: mu_in^2 = mu2in *)


	      (* initial condition at (x,y)=(0,phi0), i.e., Q^2=mu2in Exp[I phi0]: *)
	      (* u[0,0]=ain; v[0,0]=0 *)
	      (* u0m[x] = u[x,0]; v0m[x] = v[x,0] = 0 *)
              (* u[0,phi0]=Re[ainphi0]; v[0,phi]=Im[ainphi0] *)
              (* Now integration over x first, at y=phi0 *)
              (* uphi0m[x] = u[x,phi0]; vphi0m[x] = v[x,phi0] *)
   xmin=-50.; xmax=1000000.; (* We need large xmax for the integration in rtau -- to obtain rho(t) for large t *)
uvxyphi0m=NDSolve[{
    D[uphi0m[x], x] == rebt[uphi0m[x],vphi0m[x]],
    D[vphi0m[x], x] == imbt[uphi0m[x],vphi0m[x]],
    uphi0m[0] == Re[ainphi0], vphi0m[0] == Im[ainphi0]},
  {uphi0m,vphi0m}, {x,xmin, xmax}];

uphi0mf[x0_]:=Re[Evaluate[uphi0m[x0] /. uvxyphi0m][[1]]];  (* uphi0mf[x]=u[x,phi0] *)
vphi0mf[x0_]:=Re[Evaluate[vphi0m[x0] /. uvxyphi0m][[1]]];   (* vphi0mf[x0]=v[x,phi0]=0 *)

(* Integration over y at strictly fixed x=x0: u1[y]=u[x0,y];v1[y]=v[x0,y]: *)


	      uvxy1m[x0_]:=NDSolve[{ 
  D[u1m[y], y] == -imbt[u1m[y],v1m[y]],
  D[v1m[y], y] ==  rebt[u1m[y],v1m[y]],
	      u1m[phi0] == uvxyphi0m[[1,1,2]][x0], 
	      v1m[phi0] == uvxyphi0m[[1,2,2]][x0]}, {u1m, v1m}, {y,0, Pi}];
(*	      u1m[phi0] ==  uphi0mf[x0], v1m[phi0] == vphi0mf[x0] }, *)


u1mf[x_,y_]:=If[y>0,Re[Evaluate[u1m[y] /. uvxy1m[x]][[1]]],Re[Evaluate[u1m[-y] /. uvxy1m[x]][[1]]]]; 
v1mf[x_,y_]:=If[y>0, Re[Evaluate[v1m[y] /. uvxy1m[x]][[1]]], -Re[Evaluate[v1m[-y] /. uvxy1m[x]][[1]]]];

uu1mf[x_,y_]:= If[y>0,uvxy1m[x][[1,1,2]][y],uvxy1m[x][[1,1,2]][-y]];
(* Note: uu1mf[x,y]=u1mf[x,y]=u[x,y] =Re[a(Q^2=mu2in Exp[x + I y])] *)
vv1mf[x_,y_]:= If[y>0,uvxy1m[x][[1,2,2]][y],-uvxy1m[x][[1,2,2]][-y]];
(* Note: vv1mf[x,y]=v1mf[x,y]=v[x,y] =Im[a(Q^2=mu2in Exp[x + I y])] *)

  (* When alpha_s(MZ;MSbar)=0.1181, we get v1mf[Log[|Q^2|/mu2in], 10.^-5] nonnzero for |Q^2| < |Q^2|_{br.pt} =appr 0.40364 *)
  (*--------------------------------------------------------*)
(* Now the usual definition of a(Q^2=Q2 Exp[I phi];MSbar) = aMSpert[Q2,phi]: *)

ssol = NSolve[(1/be0) 1/(xmax + Y) == u1mf[xmax, 0] + I v1mf[xmax, 0]][[1]];
CPT = Re[Y /. ssol[[1]]];

  aMSpert[Q2_,phi_]:=If[mu2in Exp[xmin] < Q2 && Q2 < mu2in Exp[xmax],u1mf[Log[Q2/mu2in],phi]+ I v1mf[Log[Q2/mu2in],phi],If[Q2 <=  mu2in Exp[xmin],0,If[Q2> mu2in Exp[xmax], 1/(be0 (Log[Q2/mu2in] + CPT))]]];

(* Save["aMSpertComplexnf3almtau0XXXX.save",{u1mf,v1mf,aMSpert}] *)
