(* ::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:: *)
Print["Copyright(C)2017 Yoshifumi Hyakutake"]
Print["All rights are reserved."];


(* ::Input::Initialization:: *)
NSProd /: NSProd[T1_,T2_,o_]:=
Module[{S1,S2,n1,n2,res},
S1[0]=T1 /. \[Epsilon]->0;
Do[S1[n1]=T1;Do[S1[n1]=D[S1[n1],\[Epsilon]],{n2,1,n1}];S1[n1]=Expand[1/(n1)! (S1[n1] /. \[Epsilon]->0)];,{n1,1,o}];
S2[0]=T2 /. \[Epsilon]->0;
Do[S2[n1]=T2;Do[S2[n1]=D[S2[n1],\[Epsilon]],{n2,1,n1}];S2[n1]=Expand[1/(n1)! (S2[n1] /. \[Epsilon]->0)];,{n1,1,2}];
res=Sum[\[Epsilon]^n1 Sum[S1[n2].S2[n1-n2],{n2,0,n1}],{n1,0,o}];
Expand[res]
];


(* ::Input::Initialization:: *)
VielbeinPt /: VielbeinPt[o_]:=
Module[{T,tvbPt,ivPt,tivPt,lvbPt,gdPt,guPt,\[Mu],a},

tvbPt=Expand[Transpose[vbPt]];                    (* tvb=Subscript[e, \[Mu]]^a *)
ivPt=Expand[Normal[Series[Inverse[vbPt],{\[Epsilon],0,o}]]];                           (* iv=Subscript[e^\[Mu], a] *)
tivPt=Expand[Transpose[ivPt]];                    (* tiv=Subscript[e, a]^\[Mu] *)
lvbPt=Expand[eta.vbPt];                                    (* lvb=Subscript[e, a\[Mu]] *)
gdPt=Expand[NSProd[tvbPt.eta,vbPt,o]];                         (* gd=Subscript[g, \[Mu]\[Nu]] *)
guPt=Expand[NSProd[ivPt.eta,tivPt,o]];                                  (* gu=g^\[Mu]\[Nu] *)

(* Do[T=Sum[vbPt[[a,\[Mu]]]dX[[\[Mu]]],{\[Mu],1,dim}];
     Print[ee[[a]]," = ",T];
,{a,1,dim}];      *)                                 (* Subscript[e^a, \[Mu]]dx^\[Mu] *)

{tvbPt,ivPt,tivPt,lvbPt,gdPt,guPt}
];


(* ::Input::Initialization:: *)
ConnectionPt /: ConnectionPt[o_]:=
Module[{dePt,T1,T2,T3,SdPt,SlPt,\[Mu],a,b,c},

dePt=Expand[Table[D[lvbPt,cd[[\[Mu]]]],{\[Mu],1,dim}]];                                      (* de=\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Mu]\)]
\*SubscriptBox[\(e\), \(a\[Nu]\)]\) *)
T1=NSProd[dePt,ivPt,o];
T2=NSProd[tivPt,dePt,o];
T3=NSProd[T2,ivPt,o];
SdPt=Expand[1/2Transpose[T1,{1,3,2}]-1/2Transpose[T2,{2,3,1}]-1/2T1+1/2Transpose[T2,{3,2,1}]-1/2NSProd[tvbPt,Transpose[T3,{2,1,3}],o]+1/2NSProd[tvbPt,Transpose[T3,{3,1,2}],o]];                                              (* Sd=Subscript[\[Omega], \[Mu]ab] *)
SlPt=Expand[NSProd[tivPt,SdPt,o]];                                                                                   (* Sl=Subscript[\[Omega], cab] *)

(* Do[T=Sum[SlPt[[c,a,b]]ee[[c]],{c,1,dim}];
     If[T=!=0,Print[Subscript[\[Omega],StringJoin[id[[a]],id[[b]]]]," = ",T] ];
,{a,1,dim},{b,a+1,dim}];  *)                                        (* Subscript[\[Omega], cab]e^c *)

{SdPt,SlPt}
];


(* ::Input::Initialization:: *)
RiemannPt /: RiemannPt[o_]:=
Module[{T,domgPt,omg2Pt,RdPt,RlPt,RcPt,RsPt,\[Mu],\[Nu],a,b,c,d},

T=Expand[Table[D[SdPt,cd[[\[Mu]]]],{\[Mu],1,dim}]];                                                     (* T=\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Mu]\)]
\*SubscriptBox[\(\[Omega]\), \(\[Nu]ab\)]\) *)
domgPt=Expand[Transpose[T,{3,4,1,2}]-Transpose[T,{4,3,1,2}]];   (* domg=part of Subscript[R, ab\[Mu]\[Nu]]=\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Mu]\)]
\*SubscriptBox[\(\[Omega]\), \(\[Nu]ab\)]\)-\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Nu]\)]
\*SubscriptBox[\(\[Omega]\), \(\[Mu]ab\)]\) *)
T=Expand[Transpose[SdPt,{2,1,3}]];(* Subscript[T, a\[Mu]b]=Subscript[\[Omega], \[Mu]ab] *)
omg2Pt=Expand[Transpose[NSProd[SdPt.eta,T,o],{3,1,4,2}]-Transpose[NSProd[SdPt.eta,T,o],{4,1,3,2}]];                       (* omg2=part of Subscript[R, ab\[Mu]\[Nu]]=Subscript[\[Omega], \[Mu]a]^cSubscript[T, c\[Nu]b]-Subscript[\[Omega], \[Nu]a]^cSubscript[T, c\[Mu]b] *)

RdPt=Expand[domgPt+omg2Pt];                                                                                                                             (* Rd=Subscript[R, ab\[Mu]\[Nu]] *)
RlPt=Expand[Transpose[NSProd[NSProd[tivPt,Transpose[RdPt,{2,3,1,4}],o],ivPt,o],{3,1,2,4}]];      (* Rl=Subscript[R, abcd] *)
T=Expand[RlPt.eta];
RcPt=Table[Sum[T[[a,c,b,c]],{c,1,dim}],{a,1,dim},{b,1,dim}];   (* Rc=Subscript[R, ab] *)
T=RcPt.eta;
RsPt=Sum[T[[a,a]],{a,1,dim}];                                                                                        (* Rs=R *)

(* Do[If[RlPt[[a,b,c,d]]=!=0,
Print[Subscript[R,StringJoin[id[[a]],id[[b]],id[[c]],id[[d]]]]," = ",RlPt[[a,b,c,d]] ]];
,{a,1,dim},{b,a+1,dim},{c,1,dim},{d,c+1,dim}];
Do[If[RcPt[[a,b]]=!=0,Print[Subscript[R,StringJoin[id[[a]],id[[b]]]]," = ",RcPt[[a,b]] ]];
,{a,1,dim},{b,a,dim}];
Print[R," = ",RsPt]; *)

{RdPt,RlPt,RcPt,RsPt}
];


(* ::Input::Initialization:: *)
WeylPt /: WeylPt[Tl_]:=
Module[{T,Tc,Ts,T1,T2,WlPt,a,b,c},

T=Tl.eta;
Tc=Table[Sum[T[[a,c,b,c]],{c,1,dim}],{a,1,dim},{b,1,dim}];   (* Subscript[T, ab] *)
T=Tc.eta;
Ts=Sum[T[[a,a]],{a,1,dim}];                                                                                     (* T *)
T1=Outer[Times,eta,Tc];
T2=Outer[Times,eta,eta];

WlPt=Expand[Tl-1/(dim-2) (Transpose[T1,{1,3,2,4}]-Transpose[T1,{1,4,2,3}]-Transpose[T1,{2,3,1,4}]+Transpose[T1,{2,4,1,3}])+1/((dim-1)(dim-2)) (Transpose[T2,{1,3,2,4}]-Transpose[T2,{1,4,2,3}])Ts];                                                    (* Subscript[W, abcd] *)

WlPt
];


(* ::Input::Initialization:: *)
CovDPt /: CovDPt[T_,id_,o_]:=
Module[{len,T1,T2,res,nb,P1,rot,Om},

len=Length[id];
T1=Table[D[T,cd[[\[Mu]]]],{\[Mu],1,dim}];                                   (* T1=\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Mu]\)]
\*SubscriptBox[
SuperscriptBox[\(T\), 
SubscriptBox[\(b\), \(1\)]], \(\(
\*SubscriptBox[\(b\), \(2\)]
\*SubscriptBox[\(b\), \(3\)] ... \)
\*SubscriptBox[\(b\), \(n\)]\)]\) *)
res=NSProd[tivPt,T1,o];                          (* res=\!\(
\*SubscriptBox[\(\[PartialD]\), \(a\)]
\*SubscriptBox[
SuperscriptBox[\(T\), 
SubscriptBox[\(b\), \(1\)]], \(\(
\*SubscriptBox[\(b\), \(2\)]
\*SubscriptBox[\(b\), \(3\)] ... \)
\*SubscriptBox[\(b\), \(n\)]\)]\) *)

nb=Table[u,{u,1,len}];
(* part of Subscript[D, a]Subscript[T^Subscript[b, 1], Subscript[b, 2]Subscript[b, 3]...Subscript[b, n]]=Subscript[Subscript[\[Omega], a]^Subscript[b, 1], c]Subscript[T^c, Subscript[b, 2]Subscript[b, 3]...Subscript[b, n]]+Subscript[\[Omega], Subscript[ab, 2]]^cSubscript[T^Subscript[b, 1], Subscript[cb, 3]...Subscript[b, n]]+... *)
Do[P1=Table[v,{v,1,len}] /. {1->u,u->1};
    T1=Transpose[T,P1];
    If[id[[u]]==1,
         Om=Transpose[Transpose[SlPt,{1,3,2}].eta,{1,3,2}];
         T2=Expand[NSProd[Om,T1,o]];,                    (*  Subscript[Subscript[\[Omega], a]^Subscript[b, 1], c]=Subscript[X, acd]\[Eta]^Subscript[db, 1]=Subscript[Subscript[\[Omega], a], dc]\[Eta]^Subscript[db, 1] *)
         T2=Expand[NSProd[SlPt.eta,T1,o]];
  ];
    P1=Table[v,{v,1,len+1}] /. {2->u+1,u+1->2};
    T2=Transpose[T2,P1];
    res=Expand[res+T2];
,{u,1,len}];

res
];


(* ::Input::Initialization:: *)
TensorProdPt /: TensorProdPt[{T1_,T2_},Met_,o_]:=
Module[{Tsr,Isv,les,Iall,led,S,Ict,leI,jd,Ps,leP,Rt,tensor,u},

Tsr={T1,T2};(* {{X,{a,c,d}},{Y,{b,e,c}}} *)
{Isv,Iall,S}=Table[{0,0},{u,1,3}];
Ict=Intersection[T1[[2]],T2[[2]]];(* {c} *)
leI=Length[Ict];(* 1 *)
jd=Length[T1[[2]]]Length[T2[[2]]];                   (* 9; note jd=0 for a scalar *)

If[leI==0&&jd==0,
     tensor=Normal[Series[T1[[1]]T2[[1]],{\[Epsilon],0,o}]];
     Isv=Join[T1[[2]],T2[[2]]];
];

If[leI==0&&jd=!=0,
     tensor=Normal[Series[Outer[Times,T1[[1]],T2[[1]]],{\[Epsilon],0,o}]];
     Isv=Join[T1[[2]],T2[[2]]];
];

If[leI=!=0,
Do[Ps={};Do[Ps=Join[Ps,Position[Tsr[[u,2]],Ict[[v]]]],{v,1,leI}];(* {2} *)
     Isv[[u]]=Delete[Tsr[[u,2]],Ps];(* {a,d} *)
     Iall[[u]]=Join[Ict,Isv[[u]]];(* {c,a,d} *)
     Rt={};Do[Rt=Join[Rt,Position[Iall[[u]],Tsr[[u,2,v]]]],{v,1,Length[Iall[[u]]]}];
     Rt=Flatten[Rt,1];(* {2,1,3} *)
     S[[u]]=Transpose[Tsr[[u,1]],Rt]; (* exchange indices. contracted ones are front *)
,{u,1,2}];(* Subscript[X', cad]=Subscript[X, acd] *)

Do[Rt=Table[v,{v,1,Length[Iall[[2]]]}] /. {1->u,u->1};
     S[[2]]=NSProd[Met,Transpose[S[[2]],Rt],o];
     S[[2]]=Transpose[S[[2]],Rt];(* lift contracted indices in the second tensor *)
,{u,1,leI}];

Do[S[[u]]=Flatten[S[[u]],leI-1];,{u,1,2}];
Rt=RotateRight[Table[u,{u,1,Length[Isv[[1]]]+1}],1];(* {2,1} *)
S[[1]]=Transpose[S[[1]],Rt];
tensor=Expand[NSProd[S[[1]],S[[2]],o]];                               (* contract indices *)
Isv=Flatten[Isv,1];
];

{tensor,Isv}
];


(* ::Input::Initialization:: *)
TensorProdPt /: TensorProdPt[{T1_,T2_,T3_},Met_,o_]:=
Module[{T,len=0,nb={1,2},lentmp,od},

T={T1,T2,T3};
Do[lentmp=Length[Intersection[T[[u,2]],T[[v,2]]]];
    If[len<lentmp,len=lentmp;nb={u,v};];
,{u,1,2},{v,u+1,3}];

od=Join[nb,Complement[{1,2,3},nb]];
T={T[[od[[1]]]],T[[od[[2]]]],T[[od[[3]]]]};

TensorProdPt[Join[{TensorProdPt[{T[[1]],T[[2]]},Met,o]},{T[[3]]}],Met,o]
];


(* ::Input::Initialization:: *)
TensorProdPt /: TensorProdPt[{T1_,T2_,T3_,T4_},Met_,o_]:=
Module[{T,len=0,nb={1,2},lentmp,od},

T={T1,T2,T3,T4};
Do[lentmp=Length[Intersection[T[[u,2]],T[[v,2]]]];
    If[len<lentmp,len=lentmp;nb={u,v};];
,{u,1,3},{v,u+1,4}];

od=Join[nb,Complement[{1,2,3,4},nb]];
T={T[[od[[1]]]],T[[od[[2]]]],T[[od[[3]]]],T[[od[[4]]]]};

TensorProdPt[Join[{TensorProdPt[{T[[1]],T[[2]]},Met,o]},{TensorProdPt[{T[[3]],T[[4]]},Met,o]}],Met,o]
];


(* ::Input::Initialization:: *)
TensorProdPt /: TensorProdPt[{T1_,T2_,T3_,T4_,T5_},Met_,o_]:=
Module[{T,len=0,nb={1,2},lentmp,od},

T={T1,T2,T3,T4,T5};
Do[lentmp=Length[Intersection[T[[u,2]],T[[v,2]]]];
    If[len<lentmp,len=lentmp;nb={u,v};];
,{u,1,4},{v,u+1,5}];

od=Join[nb,Complement[{1,2,3,4,5},nb]];
T={T[[od[[1]]]],T[[od[[2]]]],T[[od[[3]]]],T[[od[[4]]]],T[[od[[5]]]]};

TensorProdPt[Join[{TensorProdPt[{T[[1]],T[[2]]},Met,o]},{TensorProdPt[{T[[3]],T[[4]],T[[5]]},Met,o]}],Met,o]
];


(* ::Input::Initialization:: *)
XtensorsPt /: XtensorsPt[{T1_,T2_},S_,o_]:=
Module[{Tg,T,Rt={},X,DX,DDX,res,\[Mu],aa,bb,cc,dd,ii,jj,vv,ww,xx,yy},

Do[Switch[Length[T1[[u]]],
    4,Tg[u]={S,T1[[u]]};,
    2,Tg[u]={Rc,T1[[u]]};,
    0,Tg[u]={Rs,T1[[u]]};
 ];
,{u,1,3}];
T=TensorProdPt[{Tg[1],Tg[2],Tg[3]},eta,o];

Switch[Length[T2],
4,Tg[4]=T2;,
2,Tg[4]={T2[[1]],xx,T2[[2]],yy};
      T[[1]]=Outer[Times,T[[1]],eta];
      T[[2]]=Join[T[[2]],{xx,yy}];,
0,Tg[4]={vv,xx,ww,yy};
      T[[1]]=T[[1]]*Outer[Times,eta,eta];
      T[[2]]={vv,ww,xx,yy};
];

Do[Rt=Join[Rt,Position[Tg[4],T[[2,u]]]],{u,1,Length[T[[2]]]}];Rt=Flatten[Rt,1];
T=Transpose[T[[1]],Rt];
X=Expand[(Transpose[T,{1,2,3,4}]-Transpose[T,{2,1,3,4}]-Transpose[T,{1,2,4,3}]
           +Transpose[T,{2,1,4,3}]+Transpose[T,{3,4,1,2}]-Transpose[T,{3,4,2,1}]
           -Transpose[T,{4,3,1,2}]+Transpose[T,{4,3,2,1}])/8];
If[S==Wl,X=Weyl[X];];

T=CovDPt[X,{-1,-1,-1,-1},o];
DX=Expand[TensorProdPt[{{T,{aa,cc,ii,jj,bb}},{eta,{aa,bb}}},eta,o][[1]]];     (* Subscript[D, a]Subscript[X, cij]^a *)

T=CovDPt[T,{-1,-1,-1,-1,-1},o];
DDX=Expand[TensorProdPt[{{T,{aa,cc,bb,ii,jj,dd}},{eta,{aa,bb}},{eta,{cc,dd}}},eta,o][[1]]];                (* Subscript[D, a]Subscript[D, b]Subscript[X^a, ij]^b *)

{X,DX,DDX}
];


(* ::Input::Initialization:: *)
W4EOMPt /: W4EOMPt[S_,T_,o_]:=
Module[{Tg,R4part,X,DX,DDX,RXpart=0,DDXpart=0,res,a,b,c,i,j,u},

Do[Switch[Length[S[[u]]],
    4,Tg[u]={T,S[[u]]};,
    2,Tg[u]={Rc,S[[u]]};,
    0,Tg[u]={Rs,S[[u]]};
 ];
,{u,1,4}];

R4part=TensorProdPt[{Tg[1],Tg[2],Tg[3],Tg[4]},eta,o][[1]];
Do[{X,DX,DDX}=XtensorsPt[{Delete[S,u],S[[u]]},T,o];
     RXpart=RXpart+2TensorProdPt[{{Rl,{i,a,b,c}},{X,{j,a,b,c}}},eta,o][[1]];
     DDXpart=DDXpart-2DDX-2Transpose[DDX];
,{u,1,4}];

res=Expand[-(R4part/2)eta+1/2 RXpart+1/2 DDXpart];
Simplify[res,TimeConstraint->1]
];
