(* ::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:: *)
SetDirectory["C:\\drorbn\\AcademicPensieve\\2017-08"];


(* ::Input::Initialization:: *)
(*$TD=\[Infinity]; \[HBar] /: \[HBar]^d_. /;d>$TD :=0;*)
$T\[Epsilon]D=2; \[Epsilon] /: \[Epsilon]^d_. /;d>$T\[Epsilon]D :=0;


(* ::Input::Initialization:: *)
q:=Sum[(\[HBar] \[Gamma] \[Epsilon])^k/k!, {k,0,Min[$T\[Epsilon]D]}];
AlgebraAtom = y|a|x;
PBWRule={y->1,a->2,x->3};
B[U@a,U@y]=-\[Gamma] U@y; B[U@x,U@a]=-\[Gamma] U@x;
(*B[U@x,U@y]=(q-1)UU[y,x]+UU@Sum[-(((t-2 a \[Epsilon])^(k+1)\[HBar]^k)/(k+1)!) ,{k,0,$TD}];*)
B[U@x,U@y]:=(q-1)UU[y,x]+UU[\[HBar]^-1 (1-Sum[T (-2 a \[Epsilon] \[HBar])^k/k! ,{k,0,$T\[Epsilon]D}])];


(* ::Input::Initialization:: *)
x_\[PrecedesEqual]y_:=OrderedQ[{x,y}/.PBWRule]; x_\[Precedes]y_:=!OrderedQ[{y,x}/.PBWRule];
Simp[\[ScriptCapitalE]_]:= Collect[\[ScriptCapitalE], _U, Expand];


(* ::Input::Initialization:: *)
Subscript[U, i_][\[ScriptCapitalE]_] := \[ScriptCapitalE] /. {t->Subscript[t, i],T->Subscript[T, i],u_U :> Replace[u, x_ :> Subscript[x, i],1]};
B[U[Subscript[(x_), i_]],U[Subscript[(y_), i_]]]:=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_,x_]=0;
B[U[y_],U[x_]]=Simp[-B[U[x],U[y]]];
B[x_,y_] := x**y-y**x;


(* ::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;


(* ::Input::Initialization:: *)
U[xx___,x_] ** U[y_,yy___] := If[x\[PrecedesEqual]y,U[xx,x,y,yy], U@xx**(U@y**U@x+B[U@x,U@y])**U@yy];


(* ::Input::Initialization:: *)
UU[c_.*(l:AlgebraAtom)^n_,r___]/;FreeQ[c, AlgebraAtom]:= Expand[c UU[Sequence@@Table[l,{n}],r]];
UU[c_.*l:AlgebraAtom,r___] := Expand[c U[l]**UU[r]];
UU[c_,r___]/;FreeQ[c, AlgebraAtom]:= Expand[c UU[r]];
UU[]=U[];
UU[l_Plus, r___] := UU[#,r]& /@ l;
UU[l_,r___] := UU[Expand[l],r];


(* ::Input::Initialization:: *)
\[DoubleStruckCapitalO][poly_, specs___] := Module[{vs,us,z},
vs=Join@@(First /@ {specs});
us=Join@@({specs} /. (l_->s_) :> (l/.Subscript[x_, i_]:>Subscript[x, s]));
Simp@Total[CoefficientRules[Normal@Series[poly, {\[HBar],0,$TD}], vs] /. (p_->c_) :> c UU@@(us^p)]
]


(* ::Input::Initialization:: *)
Subscript[\[CapitalLambda], k_][h_,e_,l_,f_,\[Alpha]_,\[Beta]_,\[Delta]_]:=Module[
{\[Rho]h,\[Rho]e,\[Rho]l,\[Rho]f,eqn,a,b,c,sol,\[Lambda],q,\[Nu]},
\[Rho]h=({
 {-1, 0},
 {0, -1}
});  \[Rho]e=({
 {0, 0},
 {-\[Epsilon], 0}
});  \[Rho]l=({
 {-(1+1/\[Epsilon])/2, 0},
 {0, (1-1/\[Epsilon])/2}
}); \[Rho]f=({
 {0, 1},
 {0, 0}
});
eqn=MatrixExp[\[Alpha] \[Rho]f].MatrixExp[\[Beta] \[Rho]e]==MatrixExp[a \[Rho]e].MatrixExp[c (\[Rho]h-2\[Epsilon] \[Rho]l)].MatrixExp[b \[Rho]f];
Echo[sol=Solve[Thread[Flatten/@eqn], {a,b,c}][[1]] /. C[1]->0];
\[Lambda]=Simplify[E^(-f \[Alpha]-e \[Beta]+h \[Alpha] \[Beta]) Normal@Series[E^(c h + a e - 2\[Epsilon] c l + b f)/.sol, {\[Epsilon],0,k}]];
q=E^(\[Nu](f \[Alpha]+e \[Beta]-h \[Alpha] \[Beta]+e f \[Delta]));
Collect[q^-1 Subscript[DP, \[Alpha]->Subscript[D, f],\[Beta]->Subscript[D, e]][\[Lambda]][q] /. \[Nu]->(1+h \[Delta])^-1,\[Epsilon],Simplify]];
