<<Gravity.m;<<Gravity2.m;<<GravityPt.m;

(* 1. Background (6) *)

(* Vielbein Formalism *)
id={"t","(Subscript[x, 1])","(Subscript[x, 2])","(Subscript[x, 3])","(Subscript[y, 1])","(Subscript[y, 2])","(Subscript[y, 3])","(Subscript[y, 4])","(Subscript[y, 5])","(Subscript[y, 6])","(Subscript[y, 7])"};(* local Lorentz indices *)
dim=Length[id];(* dimension *)
cd={t,(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[y, 1]),(Subscript[y, 2]),(Subscript[y, 3]),(Subscript[y, 4]),(Subscript[y, 5]),(Subscript[y, 6]),(Subscript[y, 7])};(* general coordinates *)
dX={dt,(Subscript[dx, 1]),(Subscript[dx, 2]),(Subscript[dx, 3]),(Subscript[dy, 1]),(Subscript[dy, 2]),(Subscript[dy, 3]),(Subscript[dy, 4]),(Subscript[dy, 5]),(Subscript[dy, 6]),(Subscript[dy, 7])};(* 1-forms of the general coordinates *)
ee=Table[(e^(id[[u]])),{u,1,dim}];(* 1-forms of the vielbein *)
eta=DiagonalMatrix[Join[{-1},Table[1,{u,1,dim-1}]]];(* local flat metric *)
vb=(({{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0,  a[t], 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, a[t] , 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, a[t] , 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, b[t], 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, b[t], 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, b[t], 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, b[t], 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, b[t], 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, b[t], 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, b[t]}}));
Phi=0;
{tvb,iv,tiv,lvb,gd,gu}=Vielbein;
{Sd,Sl}=Connection;
{Rd,Rl,Rc,Rs}=Riemann;
Wl=Weyl[Rl];
Eg=Rc-(1/2)eta Rs;

(* 2. Tensor Perturbation (15) *)

(* Vielbein Formalism *)
id={"t","(Subscript[x, 1])","(Subscript[x, 2])","(Subscript[x, 3])","(Subscript[y, 1])","(Subscript[y, 2])","(Subscript[y, 3])","(Subscript[y, 4])","(Subscript[y, 5])","(Subscript[y, 6])","(Subscript[y, 7])"};(* local Lorentz indices *)
dim=Length[id];(* dimension *)
cd={t,(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[y, 1]),(Subscript[y, 2]),(Subscript[y, 3]),(Subscript[y, 4]),(Subscript[y, 5]),(Subscript[y, 6]),(Subscript[y, 7])};(* general coordinates *)
dX={dt,(Subscript[dx, 1]),(Subscript[dx, 2]),(Subscript[dx, 3]),(Subscript[dy, 1]),(Subscript[dy, 2]),(Subscript[dy, 3]),(Subscript[dy, 4]),(Subscript[dy, 5]),(Subscript[dy, 6]),(Subscript[dy, 7])};(* 1-forms of the general coordinates *)
ee=Table[(e^(id[[u]])),{u,1,dim}];(* 1-forms of the vielbein *)
eta=DiagonalMatrix[Join[{-1},Table[1,{u,1,dim-1}]]];(* local flat metric *)
vbPt=(({{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0,  a[t](1+(ϵ/2) hp[t,(Subscript[x, 3])]), 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, a[t]ϵ hc[t,(Subscript[x, 3])], a[t](1-(ϵ/2) hp[t,(Subscript[x, 3])]) , 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, a[t] , 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, b[t], 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, b[t], 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, b[t], 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, b[t], 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, b[t], 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, b[t], 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, b[t]}}));
{tvbPt,ivPt,tivPt,lvbPt,gdPt,guPt}=VielbeinPt[1];
{SdPt,SlPt}=ConnectionPt[1];
{RdPt,RlPt,RcPt,RsPt}=RiemannPt[1];
WlPt=WeylPt[RlPt];
EgPt=Expand[RcPt-(1/2)eta RsPt];

gdPt /. {f_[t,(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[y, 1]),(Subscript[y, 2]),(Subscript[y, 3]),(Subscript[y, 4]),(Subscript[y, 5]),(Subscript[y, 6]),(Subscript[y, 7])]->f,f_[t]->f} // MatrixForm

(* 3. Evaluate (Subscript[X, abcd]) (4) *)

XTns0={TensorProdPt[{{WlPt,{a5,a6,a7,a8}},{WlPt,{a5,a6,a7,a8}}},eta,1],TensorProdPt[{{WlPt,{a4,a6,a7,a8}},{WlPt,{a5,a6,a7,a8}}},eta,1],TensorProdPt[{{WlPt,{a3,a4,a7,a8}},{WlPt,{a5,a6,a7,a8}}},eta,1],TensorProdPt[{{WlPt,{a2,a6,a4,a8}},{WlPt,{a5,a6,a7,a8}}},eta,1]};

XTensor /: XTensor:=
Module[{T1,Cf,eta2,T2,Cf2,Rt,Xprime,Xtns,DX,DDX,RXpart},
      
T1={{Normal[Series[XTns0[[1,1]]WlPt,{ϵ,0,1}]],{a1,a2,a3,a4}},TensorProdPt[{{WlPt,{a1,a2,a3,a5}},{XTns0[[2,1]],XTns0[[2,2]]}},eta,1],TensorProdPt[{{WlPt,{a1,a2,a5,a6}},{XTns0[[3,1]],XTns0[[3,2]]}},eta,1],TensorProdPt[{{WlPt,{a1,a5,a3,a7}},{XTns0[[4,1]],XTns0[[4,2]]}},eta,1],TensorProdPt[{{WlPt,{a1,a2,a5,a7}},{XTns0[[4,1]],(XTns0[[4,2]] /. {a2->a3,a4->a5,a5->a4,a7->a7})}},eta,1],TensorProdPt[{{WlPt,{a5,a2,a6,a4}},{XTns0[[3,1]],(XTns0[[3,2]] /. {a3->a1,a4->a5,a5->a3,a6->a6})}},eta,1],TensorProdPt[{{WlPt,{a1,a2,a5,a6}},{XTns0[[3,1]],(XTns0[[3,2]] /. {a3->a3,a4->a5,a5->a4,a6->a6})}},eta,1]};

Cf=96{1,-16,2,16,-16,-16,8};
Xprime=0;
Do[Rt={};
     Do[Rt=Join[Rt,Position[{a1,a2,a3,a4},T1[[u,2,v]]]],{v,1,4}];
     Rt=Flatten[Rt,1];
     Xprime=Xprime+Cf[[u]]Transpose[T1[[u,1]],Rt];
,{u,1,Length[Cf]}];

Xtns=Expand[(1/8)(Xprime-Transpose[Xprime,{2,1,3,4}]-Transpose[Xprime,{1,2,4,3}]+Transpose[Xprime,{2,1,4,3}]+Transpose[Xprime,{3,4,1,2}]-Transpose[Xprime,{4,3,1,2}]-Transpose[Xprime,{3,4,2,1}]+Transpose[Xprime,{4,3,2,1}])];
Xtns
];

(* 4. Evaluate (Subscript[Y, abcd]) (5) *)

Xtns=XTensor;

WXlPt=WeylPt[Xtns];    (* W(Subscript[(X), abcd]) up to O((ϵ^2)) *)
WXl=WXlPt /. ϵ->0;      (* W(Subscript[(X), abcd])(Subscript[|, ϵ=0]) *)

D[WXlPt[[1,2,1,2]],ϵ]/. {f_[t,(Subscript[x, 3])]->f,f_[t]->f}

(* 5. Evaluate (Subscript[δE, ab]) (13) *)

δRc=D[RcPt,ϵ];δRs=D[RsPt,ϵ];δRl=D[RlPt,ϵ];δWXl=D[WXlPt,ϵ];

ECEOMW4Pt /: ECEOMW4Pt:=
Module[{δR,δW4,δRW,T1,DPtW,DDPtW,DWPt,DDWPt,δDDWX,res},

       δR=δRc-(1/2)eta δRs;

δW4=TensorProd[{{δRl,{a,b,c,d}},{WXl,{a,b,c,d}}},eta][[1]];
δRW=TensorProd[{{δRl,{a,b,c,i}},{WXl,{a,b,c,j}}},eta][[1]]+TensorProd[{{Rl,{a,b,c,i}},{δWXl,{a,b,c,j}}},eta][[1]];

T1=CovDPt[WXl,{-1,-1,-1,-1},1];
DPtW=TensorProdPt[{{T1,{b,a,i,j,c}},{eta,{b,c}}},eta,1][[1]];(* (Subscript[D, b])W((Subscript[(X), aij])^b) *)
T1=CovDPt[DPtW,{-1,-1,-1},1];
DDPtW=TensorProdPt[{{T1,{b,a,i,j}},{eta,{a,b}}},eta,1][[1]];(* (Subscript[D, a])(Subscript[D, b])W((Subscript[(((X))^a), ij])^b) *)

T1=CovD[δWXl,{-1,-1,-1,-1}];
DWPt=TensorProd[{{T1,{b,a,i,j,c}},{eta,{b,c}}},eta][[1]];(* (Subscript[D, b])W((Subscript[(δX), aij])^b) *)
T1=CovD[DWPt,{-1,-1,-1}];
DDWPt=TensorProd[{{T1,{b,a,i,j}},{eta,{a,b}}},eta][[1]];(* (Subscript[D, a])(Subscript[D, b])W((Subscript[(((δX))^a), ij])^b) *)

δDDWX=D[DDPtW+Transpose[DDPtW],ϵ]+DDWPt+Transpose[DDWPt];

res=Expand[{δR,δW4,-(1/2)eta δW4+δRW-δDDWX}]
];

ECEOMW4Pt>>"TensorPt.m";

EOMW4Pt=<<"TensorPt.m";

(* 6. Change variables *)

rep1={a'[t]->a[t]H[t],a''[t]->a[t]((H[t])^2)+a[t]H'[t],Derivative[3][a][t]->a[t]((H[t])^3)+3a[t]H[t]H'[t]+a[t]H''[t],Derivative[4][a][t]->a[t] ((H[t])^4)+6 a[t] ((H[t])^2) H'[t]+3 a[t] ((H'[t])^2)+4 a[t] H[t] H''[t]+a[t] Derivative[3][H][t],b'[t]->b[t]G[t],b''[t]->b[t]((G[t])^2)+b[t]G'[t],Derivative[3][b][t]->b[t]((G[t])^3)+3b[t]G[t]G'[t]+b[t]G''[t],Derivative[4][b][t]->b[t] ((G[t])^4)+6 b[t] ((G[t])^2) G'[t]+3 b[t] ((G'[t])^2)+4 b[t] G[t] G''[t]+b[t] Derivative[3][G][t]};

EOMPt=Expand[D[EgPt,ϵ]+Γ EOMW4Pt[[3]] //. rep1];

rep2=Flatten[Join[Table[{Derivative[i,0][f_][t,(Subscript[x, 3])]->Derivative[i][f][t]},{i,1,4}],Table[{Derivative[i,2][f_][t,(Subscript[x, 3])]->-(k^2)Derivative[i][f][t],Derivative[i,4][f_][t,(Subscript[x, 3])]->(k^4)Derivative[i][f][t]},{i,0,4}]]];

field={k,((a[t])^(n_)),((b[t])^(n_)),hp[t],hc[t],hp'[t],hc'[t],hp''[t],hc''[t],Derivative[3][hp][t],Derivative[3][hc][t],Derivative[4][hp][t],Derivative[4][hc][t],Γ};

rep3={Derivative[4][H][t]-> (((((-1+(Sqrt[21]))/2)(Subscript[H, 0])))^4)(24((Subscript[H, 0])/(τ^5))+5040Γ ((Subscript[H, 0])^7)((13824 (477087-97732 (Sqrt[21])))/(8575 (τ^11)))),Derivative[4][G][t]->((-7+(Sqrt[21]))/14)(((((-1+(Sqrt[21]))/2)(Subscript[H, 0])))^4)(24((Subscript[H, 0])/(τ^5))+5040Γ ((Subscript[H, 0])^7)((20736 (200843-34423 (Sqrt[21])))/(8575 (τ^11)))),Derivative[3][H][t]-> (((((-1+(Sqrt[21]))/2)(Subscript[H, 0])))^3)(-6((Subscript[H, 0])/(τ^4))-504Γ ((Subscript[H, 0])^7)((13824 (477087-97732 (Sqrt[21])))/(8575 (τ^10)))),Derivative[3][G][t]->((-7+(Sqrt[21]))/14)(((((-1+(Sqrt[21]))/2)(Subscript[H, 0])))^3)(-6((Subscript[H, 0])/(τ^4))-504Γ ((Subscript[H, 0])^7)((20736 (200843-34423 (Sqrt[21])))/(8575 (τ^10)))),H''[t]-> (((((-1+(Sqrt[21]))/2)(Subscript[H, 0])))^2)(2((Subscript[H, 0])/(τ^3))+56Γ ((Subscript[H, 0])^7)((13824 (477087-97732 (Sqrt[21])))/(8575 (τ^9)))),G''[t]->((-7+(Sqrt[21]))/14)(((((-1+(Sqrt[21]))/2)(Subscript[H, 0])))^2)(2((Subscript[H, 0])/(τ^3))+56Γ ((Subscript[H, 0])^7)((20736 (200843-34423 (Sqrt[21])))/(8575 (τ^9)))),H'[t]-> ((-1+(Sqrt[21]))/2)(Subscript[H, 0])(-((Subscript[H, 0])/(τ^2))-7Γ ((Subscript[H, 0])^7)((13824 (477087-97732 (Sqrt[21])))/(8575 (τ^8)))),G'[t]->((-7+(Sqrt[21]))/14)((-1+(Sqrt[21]))/2)(Subscript[H, 0])(-((Subscript[H, 0])/(τ^2))-7Γ ((Subscript[H, 0])^7)((20736 (200843-34423 (Sqrt[21])))/(8575 (τ^8)))),H[t]->((Subscript[H, 0])/τ)+Γ ((Subscript[H, 0])^7)((13824 (477087-97732 (Sqrt[21])))/(8575 (τ^7))),G[t]->((-7+(Sqrt[21]))/14)(((Subscript[H, 0])/τ)+Γ ((Subscript[H, 0])^7)((20736 (200843-34423 (Sqrt[21])))/(8575 (τ^7))))};

(* 7. EOM (16) *)

Do[If[EOMPt[[n1,n2]]=!=0,Print[{n1,n2}];],{n1,1,11},{n2,1,11}]

Expand[EOMPt[[2,2]]+EOMPt[[3,3]]]

EOMPt[[2,3]]

Expand[EOMPt[[2,2]]-(EOMPt[[2,3]] /. hc->hp)]

E1=Collect[Expand[2EOMPt[[2,2]] //. rep2 /. {f_[t,(Subscript[x, 3])]->f[t]}],field];
E1=Sum[Collect[Simplify[E1[[i]]],field],{i,1,Length[E1]}]

E1S=Collect[Sum[Collect[FullSimplify[Normal[Series[Expand[E1[[i]] //. rep3],{Γ,0,1}]]],field],{i,1,Length[E1]}],field]

E2=Collect[Expand[2EOMPt[[2,3]] //. rep2 /. {f_[t,(Subscript[x, 3])]->f[t]}],field];
E2=Sum[Collect[Simplify[E2[[i]]],field],{i,1,Length[E2]}]

E2S=Collect[Sum[Collect[FullSimplify[Normal[Series[Expand[E2[[i]] //. rep3],{Γ,0,1}]]],field],{i,1,Length[E2]}],field]