(* ::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:: *)
AffineConnection /: AffineConnection[Gd0_]:=
Module[{J1,Gu0,dG0,\[CapitalGamma]d0,\[CapitalGamma]u0},

Gu0=Simplify[Inverse[Gd0]];                                                            (* Gu=G^JK *)
dG0=Table[D[Gd0,cd[[J1]]],{J1,1,dim}];                             (* dG=\!\(
\*SubscriptBox[\(\[PartialD]\), \(J\)]
\*SubscriptBox[\(G\), \(KL\)]\) *)
\[CapitalGamma]d0=1/2Transpose[dG0,{2,1,3}]+1/2Transpose[dG0,{3,1,2}]-1/2dG0;(* \[CapitalGamma]d=Subscript[\[CapitalGamma], JKL]=1/2\!\(
\*SubscriptBox[\(\[PartialD]\), \(K\)]
\*SubscriptBox[\(G\), \(JL\)]\)+1/2\!\(
\*SubscriptBox[\(\[PartialD]\), \(L\)]
\*SubscriptBox[\(G\), \(JK\)]\)-1/2\!\(
\*SubscriptBox[\(\[PartialD]\), \(J\)]
\*SubscriptBox[\(G\), \(KL\)]\) *)
\[CapitalGamma]u0=Gu0.\[CapitalGamma]d0;                                                                                                                  (* \[CapitalGamma]u=Subscript[\[CapitalGamma]^J, KL] *)
{Gu0,\[CapitalGamma]d0,\[CapitalGamma]u0}=Simplify[PowerExpand[{Gu0,\[CapitalGamma]d0,\[CapitalGamma]u0}],TimeConstraint->1];

{Gu0,\[CapitalGamma]d0,\[CapitalGamma]u0}
];


(* ::Input::Initialization:: *)
AffineRiemann /: AffineRiemann[Gd0_,\[CapitalGamma]u0_]:=
Module[{Gu0,T,domg,omg2,Ru0,Rd0,Rc0,Rs0,i1,i2,i,j,a1,a2,a3,a4},

Gu0=Simplify[Inverse[Gd0]];
T=Table[D[\[CapitalGamma]u0,cd[[\[Mu]1]]],{\[Mu]1,1,dim}];                                                       (* T=\!\(
\*SubscriptBox[\(\[PartialD]\), \(i\)]\*
SubscriptBox[
SuperscriptBox[
StyleBox["\[CapitalGamma]", "Subsubtitle"], "k"], "jl"]\) *)
domg=Transpose[T,{3,1,4,2}]-Transpose[T,{4,1,3,2}];           (* Subscript[d\[CapitalGamma]^k, lij]=\!\(
\*SubscriptBox[\(\[PartialD]\), \(i\)]\*
SubscriptBox[
SuperscriptBox[
StyleBox["\[CapitalGamma]", "Subsubtitle"], "k"], "jl"]\)-\!\(
\*SubscriptBox[\(\[PartialD]\), \(j\)]\*
SubscriptBox[
SuperscriptBox[
StyleBox["\[CapitalGamma]", "Subsubtitle"], "k"], "il"]\) *)
omg2=Transpose[\[CapitalGamma]u0.\[CapitalGamma]u0,{1,3,4,2}]-Transpose[\[CapitalGamma]u0.\[CapitalGamma]u0,{1,4,3,2}];           (* Subscript[\[CapitalGamma]\[CapitalGamma]^k, lij]=Subscript[\[CapitalGamma]^k, im]Subscript[\[CapitalGamma]^m, jl]-Subscript[\[CapitalGamma]^k, jm]Subscript[\[CapitalGamma]^m, il] *)

Ru0=domg+omg2;                                                            (* Ru=Subscript[R^k, lij] *)
Rd0=Gd0.Ru0;                                                                  (* Rd=Subscript[R, klij] *)
Rc0=Expand[TensorProd[{{Gd0,{i1,i2}},{Rd0,{i1,i,i2,j}}},Gu0][[1]]]; (* Rc=Subscript[R, ij] *)
Rs0=Expand[TensorProd[{{Gd0,{i1,i2}},{Rc0,{i1,i2}}},Gu0][[1]]];             (* Rs=R *)
{Ru0,Rd0,Rc0,Rs0}=Simplify[PowerExpand[{Ru0,Rd0,Rc0,Rs0}],TimeConstraint->1];

Do[If[Rd0[[a1,a2,a3,a4]]=!=0,
Print[Subscript[R,StringJoin[id[[a1]],id[[a2]],id[[a3]],id[[a4]]]]," = ",Rd0[[a1,a2,a3,a4]] ]];
,{a1,1,dim},{a2,a1+1,dim},{a3,1,dim},{a4,a3+1,dim}];
Do[If[Rc0[[a1,a2]]=!=0,Print[Subscript[R,StringJoin[id[[a1]],id[[a2]]]]," = ",Rc0[[a1,a2]] ]];
,{a1,1,dim},{a2,a1,dim}];
Print[R," = ",Rs0];

{Ru0,Rd0,Rc0,Rs0}
];


(* ::Input::Initialization:: *)
AffineWeyl /: AffineWeyl[Rd0_]:=
Module[{Rc0,Rs0,Wd0,i1,i2,i,j},

Rc0=Expand[TensorProd[{{Gd,{i1,i2}},{Rd0,{i1,i,i2,j}}},Gu][[1]]];
Rs0=Expand[TensorProd[{{Gd,{i1,i2}},{Rc0,{i1,i2}}},Gu][[1]]];
 
Wd0=Table[Rd0[[a,b,c,d]]-(Gd[[a,c]]Rc0[[b,d]]-Gd[[a,d]]Rc0[[b,c]]-Gd[[b,c]]Rc0[[a,d]]+Gd[[b,d]]Rc0[[a,c]])/(dim-2)+(Gd[[a,c]]Gd[[b,d]]-Gd[[a,d]]Gd[[b,c]])Rs0/((dim-1)(dim-2)),{a,1,dim},{b,1,dim},{c,1,dim},{d,1,dim}];

Wd0=Simplify[PowerExpand[Wd0],TimeConstraint->1];
Wd0
];


(* ::Input::Initialization:: *)
CovD\[CapitalGamma] /: CovD\[CapitalGamma][\[CapitalGamma]u0_,T_,id_]:=
Module[{len,T1,T2,res,nb,P1,\[Mu]0,u0,v0},

len=Length[id];
res=Table[D[T,cd[[\[Mu]0]]],{\[Mu]0,1,dim}];                                                 (* res=\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Mu]\)]
\*SubscriptBox[
SuperscriptBox[\(T\), 
SubscriptBox[\(\[Nu]\), \(1\)]], \(\(
\*SubscriptBox[\(\[Nu]\), \(2\)]
\*SubscriptBox[\(\[Nu]\), \(3\)] ... \)
\*SubscriptBox[\(\[Nu]\), \(n\)]\)]\) *)
nb=Table[u0,{u0,1,len}];

(* part of Subscript[D, \[Mu]]Subscript[T^Subscript[\[Nu], 1], Subscript[\[Nu], 2]Subscript[\[Nu], 3]...Subscript[\[Nu], n]]=Subscript[\[CapitalGamma]^Subscript[\[Nu], 1], \[Mu]\[Rho]]Subscript[T^\[Rho], Subscript[\[Nu], 2]Subscript[\[Nu], 3]...Subscript[\[Nu], n]]-Subscript[\[CapitalGamma]^\[Rho], Subscript[\[Mu]\[Nu], 2]]Subscript[T^Subscript[\[Nu], 1], Subscript[\[Rho]\[Nu], 3]...Subscript[\[Nu], n]]-... *)
Do[P1=Table[v0,{v0,1,len}] /. {1->u0,u0->1};
    T1=Transpose[T,P1];
    If[id[[u0]]==1,T2=\[CapitalGamma]u0.T1;,T2=-Transpose[\[CapitalGamma]u0,{3,2,1}].T1;];
    P1=Table[v0,{v0,1,len+1}] /. {1->u0+1,2->1,u0+1->2};(*(\[Mu],Subscript[\[Nu], 1],Subscript[\[Nu], 2],...,Subscript[\[Nu], u],...,Subscript[\[Nu], n])=(Subscript[\[Nu], u],\[Mu],Subscript[\[Nu], 2],...,Subscript[\[Nu], 1],...,Subscript[\[Nu], n])*)
    T2=Transpose[T2,P1];
    res=res+T2;
,{u0,1,len}];

res=Simplify[res,TimeConstraint->1];
res
];


(* ::Input::Initialization:: *)
AffineXtensors /: AffineXtensors[{T1_,T2_}]:=
Module[{Tg,T,Rt={},X0,DDX0,i1,i2,i3,i4,i5,i6},

Do[Switch[Length[T1[[u]]],
4,Tg[u]={ARd,T1[[u]]};,
2,Tg[u]={ARc,T1[[u]]};,
0,Tg[u]={ARs,T1[[u]]};
 ];
,{u,1,3}];
T=TensorProd[{Tg[1],Tg[2],Tg[3]},Gu];

(* X tensor Subscript[X, \[Mu]\[Nu]\[Rho]\[Sigma]], Subscript[X, \[Mu]\[Nu]], X *)
Switch[Length[T2],
4,Do[Rt=Join[Rt,Position[T2,T[[2,u]]][[1]]],{u,1,4}];
      X0=Transpose[T[[1]],Rt];
      X0=Expand[(X0-Transpose[X0,{2,1,3,4}]-Transpose[X0,{1,2,4,3}]
  +Transpose[X0,{2,1,4,3}]+Transpose[X0,{3,4,1,2}]-Transpose[X0,{4,3,1,2}]
  -Transpose[X0,{3,4,2,1}]+Transpose[X0,{4,3,2,1}])/8];,
2,X0=Expand[(T[[1]]+Transpose[T[[1]]])/2];,
0,X0=T[[1]];
];

(* DDX tensor Subscript[D, \[Rho]]Subscript[D, \[Sigma]]Subscript[X^\[Rho], \[Mu]\[Nu]]^\[Sigma], Subscript[D, \[Rho]]Subscript[D, \[Sigma]]Subscript[X, \[Mu]\[Nu]], Subscript[D, \[Rho]]Subscript[D, \[Sigma]]X *)
Switch[Length[T2],
4,T=CovD\[CapitalGamma][\[CapitalGamma]u,X0,{-1,-1,-1,-1}];
      T=CovD\[CapitalGamma][\[CapitalGamma]u,T,{-1,-1,-1,-1,-1}];
DDX0=TensorProd[{{Gu,{i1,i3}},{Gu,{i2,i6}},{T,{i1,i2,i3,i4,i5,i6}}},IdentityMatrix[dim]][[1]];,
2,T=CovD\[CapitalGamma][\[CapitalGamma]u,X0,{-1,-1}];
      DDX0=CovD\[CapitalGamma][\[CapitalGamma]u,T,{-1,-1,-1}];,
0,T=Table[D[X0,cd[[u]]],{u,1,dim}];
      DDX0=CovD\[CapitalGamma][\[CapitalGamma]u,T,{-1}];
];

{X0,DDX0}
];


(* ::Input::Initialization:: *)
AffineR4EOM  /: AffineR4EOM[T1_]:=
Module[{Tg,ind,Ps,S,T,R4part,RXpart=0,DDXpart=0,res},

Do[Switch[Length[T1[[u]]],
4,Tg[u]={ARd,T1[[u]]};,
2,Tg[u]={ARc,T1[[u]]};,
0,Tg[u]={ARs,T1[[u]]};
 ];
,{u,1,4}];

(* variation of volume factor *)
R4part=-TensorProd[{Tg[1],Tg[2],Tg[3],Tg[4]},Gu][[1]]*Gd/2;

(* variation of metric *)
ind=Union[Flatten[T1,1]];
Do[Ps=Position[T1,ind[[u]]];
    S=ReplacePart[T1,{Ps[[1]]->\[Mu],Ps[[2]]->\[Nu]}];
T=TensorProd[{{Tg[1][[1]],S[[1]]},{Tg[2][[1]],S[[2]]},{Tg[3][[1]],S[[3]]},{Tg[4][[1]],S[[4]]}},Gu][[1]];
    RXpart=RXpart+(T+Transpose[T])/2;
,{u,1,Length[ind]}];

(* variation of Riemann tensor *)
Do[S=AffineXtensors[{Delete[T1,u],T1[[u]]}];
    Switch[Length[T1[[u]]],
    4,T=TensorProd[{{ARd,{\[Mu],a1,a2,a3}},{S[[1]],{\[Nu],a1,a2,a3}}},Gu][[1]];
          RXpart=RXpart-(T+Transpose[T])/2;
          DDXpart=DDXpart-S[[2]]-Transpose[S[[2]]];,
    2,T=-ARc.Gu.S[[1]]/2
           +TensorProd[{{ARd,{\[Mu],a1,\[Nu],a2}},{S[[1]],{a1,a2}}},Gu][[1]]/2;
          RXpart=RXpart+(T+Transpose[T])/2;
          T=TensorProd[{{Gd,{a1,a2}},{S[[2]],{a1,a2,\[Mu],\[Nu]}}},Gu][[1]]/2
           -TensorProd[{{Gd,{a1,a2}},{S[[2]],{a1,\[Mu],a2,\[Nu]}}},Gu][[1]]/2
           -TensorProd[{{Gd,{a1,a2}},{S[[2]],{\[Mu],a1,a2,\[Nu]}}},Gu][[1]]/2
           +TensorProd[{{Gd,{a1,a2}},{Gd,{a3,a4}},{S[[2]],{a1,a3,a2,a4}}},Gu][[1]]/2*Gd;
          DDXpart=DDXpart+(T+Transpose[T])/2;,
    0,RXpart=RXpart+S[[1]]ARc;
          RXpart=RXpart-(S[[2]]+Transpose[S[[2]]])/2
                      +TensorProd[{{Gd,{a1,a2}},{S[[2]],{a1,a2}}},Gu][[1]]*Gd;
  ];
,{u,1,4}];

res=Simplify[PowerExpand[R4part+RXpart+DDXpart],TimeConstraint->1];
res
];


(* ::Input::Initialization:: *)
CompareLG /: CompareLG[Ltns_,Gtns_]:=
Module[{nb,T1,Rt0,Rt,T2,res},

nb=TensorRank[Ltns];
T1=Ltns;
Rt0=Table[u,{u,1,nb}];
Do[Rt=Rt0 /. {1->u,u->1};
    T1=tvb.Transpose[T1,Rt];
    T1=Transpose[T1,Rt];
,{u,1,nb}];

T1=Simplify[PowerExpand[T1],TimeConstraint->1];
res=Simplify[PowerExpand[T1-Gtns],TimeConstraint->1];
res=FullSimplify[res];
T2=Table[0,{u,1,dim}];
Do[T2=Table[T2,{u,1,dim}];,{u,1,nb-1}];
If[res==T2,Print["No Problem. Tensors are exactly the same."];,Print["Something is wrong."];];
];


(* ::Input::Initialization:: *)
Compare2LG /: Compare2LG[Ltns_,Gtns_,rep_]:=
Module[{nb,T1,Rt0,Rt,T2,res},

nb=TensorRank[Ltns];
T1=Ltns;
Rt0=Table[u,{u,1,nb}];
Do[Rt=Rt0 /. {1->u,u->1};
    T1=tvb.Transpose[T1,Rt];
    T1=Transpose[T1,Rt];
,{u,1,nb}];

T1=Simplify[PowerExpand[T1],TimeConstraint->1];
res=Simplify[PowerExpand[T1-Gtns],TimeConstraint->1];
res=FullSimplify[res /. rep];
res=res /. rep;
T2=Table[0,{u,1,dim}];
Do[T2=Table[T2,{u,1,dim}];,{u,1,nb-1}];
If[res==T2,Print["No Problem. Tensors are exactly the same."];,Print["Something is wrong."];];
];
