(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



(* ::Input::Initialization:: *)
\[Gamma] /: \[Gamma]^x_/;x>1 :=0;
PBWBasis={W,u,M};

B[U@u,U@M]=-(B[U@M,U@u]=U@M);
B[U@W,U@u]=-(B[U@u,U@W]=U@W);

B[U@W,U@M]=-(B[U@M,U@W]=(1-t)U[]-\[Gamma](1+t) U[u]);


(* ::Input::Initialization:: *)
UU[l___,x_^n_,r___]:= UU[l,Sequence@@Table[x,{n}],r];
UU[l___,1,r___]:= UU[l,r];
UU[]=U[]; UU[l_,r___]:=U[l]**UU[r];
Subscript[U, i_][\[ScriptCapitalE]_] := \[ScriptCapitalE] /. {t->Subscript[t, i], u_U :> UU@@Replace[u, x_ :> Subscript[x, i],1]};


(* ::Input::Initialization:: *)
B[x_,x_]=0;
B[U[Subscript[(x_), i_]],U[Subscript[(y_), i_]]]:=Subscript[U, i][B[U@x,U@y]];
B[U[Subscript[(x_), i_]],U[Subscript[(y_), j_]]] /; i=!=j := 0;
B[x_,y_] := x**y-y**x;


(* ::Input::Initialization:: *)
x_\[PrecedesEqual]y_:=OrderedQ[{y,x}]; x_\[Precedes]y_:=!OrderedQ[{x,y}]; (*Anti-alphabetic ordering*)
Simp[\[ScriptCapitalE]_]:= Collect[\[ScriptCapitalE], _U,Together];


(* ::Input::Initialization:: *)
Unprotect[NonCommutativeMultiply];
NonCommutativeMultiply[x_] := x;
0**_=_**0=0;
x_**U[]:=x; U[]**x_:=x;
(a_*x_U)**(b_*y_U) := If[a b === 0, 0, Simp[a b (x**y)]];
(a_*x_U)**y_ := Simp[a(x**y)]; x_**(a_*y_U) := Simp[a(x**y)];
(x_Plus)**y_ := (#**y)& /@ x; x_**(y_Plus) := (x**#)& /@ y;
U[x_]**U[y_] := (*U[x]**U[y] =*) If[x\[Precedes]y, U[x,y], U[y,x]+B[U@x,U@y]];
U[x_] ** U[y1_,yy__] := (*U[x] ** U[y1,yy] =*)If[x\[PrecedesEqual]y1,U[x,y1,yy],(U@x**U@y1)**U@yy];
U[xx__,xn_]**U[yy__] := (*U[xx,xn]**U[yy] =*) U@xx**(U@xn**U@yy);


(* ::Input::Initialization:: *)
\[Sigma][i_,j_][\[ScriptCapitalE]_] := \[ScriptCapitalE] /. {Subscript[W, i]->Subscript[W, j],Subscript[t, i]->Subscript[t, j],Subscript[M, i]->Subscript[M, j],Subscript[c, i]->Subscript[c, j],Subscript[u, i]->Subscript[u, j]};


(* ::Input::Initialization:: *)
mul[i_,j_][\[ScriptCapitalE]_] := Simp[\[ScriptCapitalE] /. x_U :> DeleteCases[x,Subscript[_, j]]**U@@Cases[x,Subscript[y_, j]:>Subscript[y, i]] /. {Subscript[c, j]->Subscript[c, i], Subscript[t, j]->Subscript[t, i]}];
mul[i_,j_,k_][\[ScriptCapitalE]_] := \[ScriptCapitalE] // mul[i,j] // \[Sigma][i,k];


(* ::Input::Initialization:: *)
\[CapitalDelta][i_,j_,k_][\[ScriptCapitalE]_] := Simp[\[ScriptCapitalE]/. {
z_. U[]:>(z /. {Subscript[c, i]->Subscript[c, j]+Subscript[c, k],Subscript[t, i]->Subscript[t, j] Subscript[t, k]})U[],
z_. x_U :> (z /. {Subscript[c, i]->Subscript[c, j]+Subscript[c, k],Subscript[t, i]->Subscript[t, j] Subscript[t, k]})NonCommutativeMultiply@@(x /. {
Subscript[u, i]->U@Subscript[u, j]+U@Subscript[u, k],
 Subscript[W, i]-> U[Subscript[W, j]]+\[Gamma]/2 U[Subscript[W, j],Subscript[u, k]]+U[Subscript[W, k]]-\[Gamma]/2 U[Subscript[W, k],Subscript[u, j]],
Subscript[M, i]-> Subscript[t, k]U[Subscript[M, j]]+Subscript[t, k] \[Gamma]/2 U[Subscript[u, k],Subscript[M, j]]+U[Subscript[M, k]]-\[Gamma]/2 U[Subscript[u, j],Subscript[M, k]],
 Subscript[y_, l_]:>U@Subscript[y, l]
})
}]


(* ::Input::Initialization:: *)
S[i_][\[ScriptCapitalE]_] := Simp[\[ScriptCapitalE] /. {z_. x_U:>(z /. {Subscript[c, i]->-Subscript[c, i],Subscript[t, i]->Subscript[t, i]^-1})S[i][x]}];
S[i_][U[]]=U[];
S[i_][U[Subscript[y_, j_],more___]] /; i!=j := U[Subscript[y, j]]**S[i][U[more]];   (*Careful! if mathematica cannot decide i\[Equal]j or not then we get an error*)
S[i_][U[Subscript[u, i_],more___]] := S[i][U[more]]**(-U@Subscript[u, i]);
S[i_][U[Subscript[M, i_],more___]] := S[i][U[more]]**(-Subscript[t, i]^-1(1-\[Gamma]/2)U[Subscript[M, i]]);S[i_][U[Subscript[W, i_],more___]] := S[i][U[more]]**( -(1+\[Gamma]/2)U[Subscript[W, i]]);


(* ::Input::Initialization:: *)
Sinv[i_][\[ScriptCapitalE]_] := Simp[\[ScriptCapitalE] /. {z_. x_U:>(z /. {Subscript[c, i]->-Subscript[c, i],Subscript[t, i]->Subscript[t, i]^-1})Sinv[i][x]}];
Sinv[i_][U[]]=U[];
Sinv[i_][U[Subscript[y_, j_],more___]] /; i!=j := U[Subscript[y, j]]**Sinv[i][U[more]];   (*Careful! if mathematica cannot decide i\[Equal]j or not then we get an error*)
Sinv[i_][U[Subscript[u, i_],more___]] := Sinv[i][U[more]]**(-U@Subscript[u, i]);
Sinv[i_][U[Subscript[M, i_],more___]] := Sinv[i][U[more]]**(-Subscript[t, i]^-1(1+\[Gamma]/2)U[Subscript[M, i]]);Sinv[i_][U[Subscript[W, i_],more___]] := Sinv[i][U[more]]**( -(1-\[Gamma]/2)U[Subscript[W, i]]);


(* ::Input::Initialization:: *)
CoUnit[i_][\[ScriptCapitalE]_] :=Simp[\[ScriptCapitalE]/. {z_. x_U:>(z/. {t->1,c->0,Subscript[c, i]->0,Subscript[t, i]->1})CoUnit[i][x]}] ;
CoUnit[i_][U[]]=U[];
CoUnit[i_][U[Subscript[y_, j_],more___]] /; i!=j := U[Subscript[y, j]]**CoUnit[i][U[more]];
CoUnit[i_][U[Subscript[y_, i_],more___]] /;y!=a:=0;
CoUnit[i_][U[Subscript[a, i_],more___]]:= CoUnit[i][U[more]];


(* ::Input::Initialization:: *)
UExp[\[ScriptCapitalE]_,n_]:=Module[{t=U[],k}, U[]+Sum[(t=t**\[ScriptCapitalE])/k!, {k,n}]]//Simp


(* ::Input::Initialization:: *)
ToDegree[n_][\[ScriptCapitalE]_] := (Simp[\[ScriptCapitalE]] /. {\[Gamma]->\[HBar] \[Gamma], Subscript[c, i_]:>\[HBar] Subscript[c, i], Subscript[t, i_]:>E^(\[HBar] Subscript[c, i]),c->\[HBar] c,  t:>E^(\[HBar] c),x_U:>\[HBar]^(Count[x,W|Subscript[W, _]]+Count[x,u|Subscript[u, _]]) x} /. a_. x_U :>Normal[Series[a,{\[HBar],0,n}]]*x)/.\[HBar]->1


(* ::Input::Initialization:: *)
R[i_,j_,d_]:= Sum[UU[Subscript[W, i]^a]**(U[]+\[Gamma]/2 a U[Subscript[u, i]])**UU[Subscript[u, i]^b]**(Subscript[c, j]^b U[]+\[Gamma] b Subscript[c, j]^(b-1) U[Subscript[u, j]])**(U[]+\[Gamma]/2 a U[Subscript[u, j]])**UU[Subscript[M, j]^a]/(a!b!) (1+\[Gamma] 1/4 (a-1)a),{b,0,d},{a,0,d-b}]//ToDegree[d]



(* ::Input::Initialization:: *)
Rinv[i_,j_,d_]:=S[i][R[i,j,d]]//Expand


(* ::Input::Initialization:: *)
(*\[Eta] /: \[Eta]^n_/;n>1 :=0; *) 
\[Eta]=(t+1)/(t-1) \[Epsilon];
q = 1+\[Eta];
qI[k_]:=k(1+\[Eta] (k-1)/2)//Expand
qFac[n_]:=n!(1+\[Eta] (n-1)n/4)//Expand
InvqFac[n_]:=(1-\[Eta] (n-1)n/4)/n!//Expand
qBin[n_,k_]:= Binomial[n,k](1+\[Eta] k (n-k) /2)//Expand




