Tetrad approach to General Relativity

Choose whether to perturb

(* Perturb=False is standard choice; set Perturb=True if you want to compute perturbations in Minkowski or cosmological spacetimes *)
Clear[Perturb];
Perturb=False

Algorithm for simplifying anything; this one involves taking all equations to linear order in the small parameter ϵ

Clear[S];
If[Perturb,
S[x__]:=Simplify[Normal[Series[Together[x, Trig->True],{ϵ,0,1}]],TimeConstraint->Infinity]
,(*else*)
S[x__]:=Simplify[Together[x, Trig->True],TimeConstraint->Infinity]
];

(* bug in Mathematica 11: Series[(ϵ^2),{ϵ,0,0}] gives (ϵ^2) instead of 0 *)
Series[(ϵ^2),{ϵ,0,0}]

Choose number of spacetime dimensions

(* Choose number of time and space dimensions, usually 1 and 3 *)
Clear[TimeDimension,SpaceDimension,Dimension];
TimeDimension=1;
SpaceDimension=3;

(* Number of spacetime dimensions *)
Clear[Dimension];
Dimension=TimeDimension+SpaceDimension

Choose coordinates

(* Options for time coordinate *)
Clear[Time,Conformal];
(* Choose time coordinate *)
Clear[TimeOption];
TimeOption=Time

(* Options for spatial coordinates *)
Clear[Cartesian,Polar,ComovingPolar,Cylindrical];
(* Choose spatial coordinates *)
Clear[SpaceOption];
SpaceOption=Polar

Clear[FirstIndex,LastIndex];
(* FirstIndex=0 for a time dimension, or FirstIndex=1 for only spatial dimensions *)
FirstIndex=If[TimeDimension>0,0,1];
LastIndex=SpaceDimension;

Clear[t,η,x,y,z,w,v,u,𝓍,𝓎,𝓏,𝓌,𝓋,𝓊,r,θ,ϕ,χ,ψ];
Clear[Coordinate];
Which[
TimeOption===Time,
Coordinate[0]=t;
,TimeOption===Conformal,
Coordinate[0]=η;
];
Which[
SpaceOption===Cartesian,
Coordinate[1]=x;
Coordinate[2]=y;
If[LastIndex≥3,Coordinate[3]=z];
If[LastIndex≥4,Coordinate[4]=w];
If[LastIndex≥5,Coordinate[5]=v];
If[LastIndex≥6,Coordinate[6]=u];
If[LastIndex≥7,Coordinate[7]=𝓍];
If[LastIndex≥8,Coordinate[8]=𝓎];
If[LastIndex≥9,Coordinate[9]=𝓏];
If[LastIndex≥10,Coordinate[10]=𝓌];
If[LastIndex≥11,Coordinate[11]=𝓋];
If[LastIndex≥12,Coordinate[12]=𝓊];
,SpaceOption===Polar,
Coordinate[1]=r;
Coordinate[2]=θ;
If[LastIndex≥3,Coordinate[3]=ϕ];
If[LastIndex≥4,Coordinate[4]=χ];
If[LastIndex≥5,Coordinate[5]=ψ];
,SpaceOption===ComovingPolar,
Coordinate[1]=x;
Coordinate[2]=θ;
If[LastIndex≥3,Coordinate[3]=ϕ];
,SpaceOption===Cylindrical,
Coordinate[1]=r;
Coordinate[2]=ϕ;
Coordinate[3]=z;
];
Table[Coordinate[i],{i,FirstIndex,LastIndex}]

Clear[FunctionCoordinates];
FunctionCoordinates[f_]:=(f[Table[Coordinate[i],{i,FirstIndex,LastIndex}]]/.{x__}->x)

Choose metric sign (-1 for +---, or +1 for -+++)

Clear[MetricSign];
MetricSign=1

Choose the tetrad metric

(* Options for tetrad metric *)
Clear[Minkowski,Spinor,NewmanPenrose,(*for traditional Arnowitt-Deser-Misner*)Spatial];
(* Choose tetrad metric *)
Clear[TetradMetricOption];
TetradMetricOption=Minkowski

(* Tetrad metric tensor (Subscript[γ, mn]) *)
Clear[TetradMetricMatrix];
TetradMetricMatrix=Switch[TetradMetricOption
(* Minkowski *)
,Minkowski,
Table[Which[i==0&&j==0,-1,i==j,1,True,0],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}]
(* spinor *)
,Spinor,
(({{-1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 0}}))
(* Newman-Penrose *)
,NewmanPenrose,
Table[If[i-FirstIndex<2IntegerPart[(Dimension/2)]&&j-FirstIndex<2IntegerPart[(Dimension/2)],Which[(i==0&&j==1)||(i==1&&j==0),-1,(Mod[i-FirstIndex,2]==0&&j==i+1)||(Mod[i-FirstIndex,2]==1&&j==i-1),1,True,0],
(*else*)
Which[i==0&&j==0,-1,i==j,1,True,0]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}]
(* spatial *)
,Spatial,
Table[Which[i==0&&j==0,-1,i==0||j==0,0,i>j,FunctionCoordinates[γ[j,i]],True,FunctionCoordinates[γ[i,j]]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}]
];

(* Define simplification rules for printing out quantities *)
Clear[SPrint];
Which[SpaceOption===Cartesian,
SPrint[expr__]:=Simplify[expr/.FunctionCoordinates[func_]->func//.{(x^2)+(y^2)+(z^2)->(r^2),((((x^2)+(y^2)+(z^2)))^(n_Integer))/;EvenQ[n]->(r^n)},Assumptions->r≥0,TimeConstraint->Infinity];
,True,
SPrint[expr__]:=Simplify[expr/.FunctionCoordinates[func_]->func,TimeConstraint->Infinity];
]

(* change sign of tetrad metric to have desired signature *)
If[MetricSign===-1,TetradMetricMatrix=-TetradMetricMatrix];

(* whether or not tetrad metric is constant in spacetime *)
Clear[TetradMetricIsConstant];
TetradMetricIsConstant=Switch[TetradMetricOption,
Minkowski,True,
NewmanPenrose,True,
_,False]

Clear[TetradMetric];
TetradMetric[i_,j_]:=TetradMetric[i,j]=TetradMetricMatrix[[i+1-FirstIndex,j+1-FirstIndex]];
Table[SPrint[TetradMetric[i,j]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}]//MatrixForm

Clear[InverseTetradMetricMatrix];
InverseTetradMetricMatrix=S[Inverse[TetradMetricMatrix]];

Clear[InverseTetradMetric];
InverseTetradMetric[i_,j_]:=InverseTetradMetric[i,j]=InverseTetradMetricMatrix[[i+1-FirstIndex,j+1-FirstIndex]];
Table[SPrint[InverseTetradMetric[i,j]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}]//MatrixForm

(* matrix that transforms from Minkowski to desired tetrad metric *)
Clear[TetradTransformMatrix,InverseTetradTransformMatrix];
TetradTransformMatrix=Switch[TetradMetricOption
,Minkowski,Table[KroneckerDelta[i,j],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}]
,Spinor,(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, (1/(Sqrt[2])), (I/(Sqrt[2]))}, {0, 0, (1/(Sqrt[2])), -(I/(Sqrt[2]))}}))
,NewmanPenrose,
Table[If[i-FirstIndex<2IntegerPart[(Dimension/2)]&&j-FirstIndex<2IntegerPart[(Dimension/2)],(1/(Sqrt[2]))Which[Mod[j-FirstIndex,2]==0&&(i==j||i==j+1),1,Mod[j-FirstIndex,2]==1&&i==j-1,If[FirstIndex==0&&i==0,1,I],Mod[j-FirstIndex,2]==1&&i==j,-If[FirstIndex==0&&i==1,1,I],True,0],(*else*)Which[i==0&&j==0,-1,i==j,1,True,0]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}]
,Spatial,Table[KroneckerDelta[i,j],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}]
];
InverseTetradTransformMatrix=S[Inverse[TetradTransformMatrix]];
TetradTransformMatrix//MatrixForm

Clear[DeterminantTetradTransformMatrix];
DeterminantTetradTransformMatrix=S[Det[TetradTransformMatrix]]

(* Totally antisymmetric (Levi-Civita) tensor (Subscript[ε, ijkl]) *)
Clear[TotallyAntisymmetric];
(* sign of (Subscript[ε, 0123]) in locally inertial frame *)
SignTotallyAntisymmetric=-1;
(* FIX: this is not correct for TetradMetricOption===Spatial *)
TotallyAntisymmetric[i_Integer,j_Integer,k_Integer,l_Integer]:=SignTotallyAntisymmetric DeterminantTetradTransformMatrix Signature[{i,j,k,l}]

Choose the inverse vierbein tensor (choose inverse vierbein rather than vierbein for historical reasons)

(* Options for vierbein *)
Clear[
(* Minkowski *)
Minkowski,
(* Friedmann-Lemaître-Robertson-Walker *)
FriedmannLemaitreRobertsonWalker,
(* Reissner-Nordström *)
ReissnerNordstrom,
(* Gullstrand-Painlevé *)
GullstrandPainleve,
(* general spherical geometry *)
GeneralSpherical,
(* Boyer-Lindquist *)
BoyerLindquist,
(* Doran *)
Doran,
(* Arnowitt-Deser-Misner (3+1) *)
ArnowittDeserMisner,
(* Dupin cyclide *)
DupinCyclide
];
(* Choose vierbein option *)
Clear[VierbeinOption];
VierbeinOption=GeneralSpherical

(* InverseVierbein tensor ((Subscript[e, m])^μ) *)
Clear[InverseVierbeinMatrix];
Which[
(* Minkowski *)
VierbeinOption===Minkowski,
Which[TimeOption===Time&&SpaceOption===Cartesian,
InverseVierbeinMatrix=(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}));
];

(* Friedmann-Lemaitre-Robertson-Walker *)
,VierbeinOption===FriedmannLemaitreRobertsonWalker,
(* cosmic scale factor a and FLRW curvature constant κ *)
Clear[a,κ];
Which[TimeOption===Time&&SpaceOption===Cartesian,
(* spatially flat universe *)
InverseVierbeinMatrix=(({{1, 0, 0, 0}, {0, (1/(a[t])), 0, 0}, {0, 0, (1/(a[t])), 0}, {0, 0, 0, (1/(a[t]))}}));
,TimeOption===Conformal&&SpaceOption===Cartesian,
InverseVierbeinMatrix=(1/(a[η]))(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}));
,TimeOption===Time&&SpaceOption===ComovingPolar,
(* Here x is the comoving circumferential radius *)
InverseVierbeinMatrix=(({{1, 0, 0, 0}, {0, ((Sqrt[1-κ (x^2)])/(a[t])), 0, 0}, {0, 0, (1/(a[t]x)), 0}, {0, 0, 0, (1/(a[t]x Sin[θ]))}}));
,TimeOption===Conformal&&SpaceOption===ComovingPolar,
InverseVierbeinMatrix=(1/(a[η]))(({{1, 0, 0, 0}, {0, (Sqrt[1-κ (x^2)]), 0, 0}, {0, 0, (1/x), 0}, {0, 0, 0, (1/(x Sin[θ]))}}));
,TimeOption===Time&&SpaceOption===Polar,
(* Here r is the comoving geodesic radius *)
Clear[ConformalPoint];
If[κ===0,
InverseVierbeinMatrix=(({{1, 0, 0, 0}, {0, (1/(a[t])), 0, 0}, {0, 0, (1/(a[t]r)), 0}, {0, 0, 0, (1/(a[t]r Sin[θ]))}}))
,(*else*)
ConformalPoint={κ->1,r->(π/2),θ->(π/2)};
InverseVierbeinMatrix=(({{1, 0, 0, 0}, {0, (1/(a[t])), 0, 0}, {0, 0, ((Sqrt[κ])/(a[t]Sin[(Sqrt[κ]) r])), 0}, {0, 0, 0, ((Sqrt[κ])/(a[t]Sin[(Sqrt[κ]) r]Sin[θ]))}}))
];
,TimeOption===Conformal&&SpaceOption===Polar,
Clear[ConformalPoint];
If[κ===0,
InverseVierbeinMatrix=(1/(a[η]))(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, (1/r), 0}, {0, 0, 0, (1/(r Sin[θ]))}}))
,(*else*)
ConformalPoint={κ->1,r->(π/2),θ->(π/2)};
InverseVierbeinMatrix=(1/(a[η]))(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, ((Sqrt[κ])/(Sin[(Sqrt[κ]) r])), 0}, {0, 0, 0, ((Sqrt[κ])/(Sin[(Sqrt[κ]) r]Sin[θ]))}}))
];
];

(* Reissner-Nordström *)
,VierbeinOption===ReissnerNordstrom,
Clear[B,M,Q];
(*Q=0;*)
B[r_]:=1-((2M)/r)+((Q^2)/(r^2));
Which[TimeOption===Time&&SpaceOption===Polar,
InverseVierbeinMatrix=(({{(1/(Sqrt[B[r]])), 0, 0, 0}, {0, (Sqrt[B[r]]), 0, 0}, {0, 0, (1/r), 0}, {0, 0, 0, (1/(r Sin[θ]))}}));
];

(* Gullstrand-Painlevé *)
,VierbeinOption===GullstrandPainleve,
Clear[β,M,Q];
Q=0;
(* Reissner-Nordström *)
β[r_]:=-(Sqrt[((2M)/r)-((Q^2)/(r^2))]);
Which[TimeOption===Time&&SpaceOption===Cartesian,
(* use script 𝓇 to avoid conflict with r *)
Clear[𝓇];
𝓇=(Sqrt[(x^2)+(y^2)+(z^2)]);
InverseVierbeinMatrix=(({{1, ((β[𝓇]x)/𝓇), ((β[𝓇]y)/𝓇), ((β[𝓇]z)/𝓇)}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}));
,TimeOption===Time&&SpaceOption===Polar,
InverseVierbeinMatrix=(({{1, β[r], 0, 0}, {0, 1, 0, 0}, {0, 0, (1/r), 0}, {0, 0, 0, (1/(r Sin[θ]))}}));
];

(* general spherical *)
,VierbeinOption===GeneralSpherical,
Which[FirstIndex==0,
Which[TimeOption===Time&&SpaceOption===Polar,
Clear[A,βr,βt];
A=1;
InverseVierbeinMatrix=(({{(1/(α[t,r])), βt[t,r], 0, 0}, {0, βr[t,r], 0, 0}, {0, 0, (1/r), 0}, {0, 0, 0, (1/(r Sin[θ]))}}));
,TimeOption===Time&&SpaceOption===ComovingPolar,
Clear[α,A,λt,λx,r,𝓇];
A=1;
r=𝓇[t,x];
InverseVierbeinMatrix=(({{(1/(α[t,x])), λt[t,x], 0, 0}, {0, λx[t,x], 0, 0}, {0, 0, (1/r), 0}, {0, 0, 0, (1/(r Sin[θ]))}}));
,TimeOption===Conformal&&SpaceOption===ComovingPolar,
Clear[a,A,βx,v,ξr,ξt];
a=(E^(v η));
A=a x;
InverseVierbeinMatrix=(1/A)(({{(1/(ξη[x])), -((x βx[x]ξx[x])/(ξη[x])), 0, 0}, {0, x βx[x], 0, 0}, {0, 0, 1, 0}, {0, 0, 0, (1/(Sin[θ]))}}));
];
,FirstIndex==1,
Which[SpaceOption===Polar,
Clear[βr];
InverseVierbeinMatrix=(({{βr[r], 0, 0}, {0, (1/r), 0}, {0, 0, (1/(r Sin[θ]))}}));
];
];

(* Boyer-Lindquist *)
,VierbeinOption===BoyerLindquist,
Clear[a,β,M,Q,R,ρ];
(*Q=0;*)
Which[TimeOption===Time&&SpaceOption===Polar,
R[r_]:=(Sqrt[(r^2)+(a^2)]);
ρ[r_,θ_]:=(Sqrt[(r^2)+(a^2)((Cos[θ])^2)]);
β[r_]:=-((Sqrt[2M r-(Q^2)])/(R[r]))/;(r∈Integers)=!=True;
InverseVierbeinMatrix=(({{((R[r])/(ρ[r,θ](Sqrt[1-((β[r])^2)]))), 0, 0, (a/(R[r]ρ[r,θ](Sqrt[1-((β[r])^2)])))}, {0, ((R[r](Sqrt[1-((β[r])^2)]))/(ρ[r,θ])), 0, 0}, {0, 0, (1/(ρ[r,θ])), 0}, {((a Sin[θ])/(ρ[r,θ])), 0, 0, (1/(ρ[r,θ]Sin[θ]))}}));
];

(* Doran *)
,VierbeinOption===Doran,
Clear[a,α,β,M,Q,R,ρ];
(*Q=0;*)
(* Kerr-Newman *)
(*β[r_]:=-((Sqrt[2M r-(Q^2)])/(R[r]))/;(r∈Integers)=!=True;*)
Which[TimeOption===Time&&SpaceOption===Cartesian,
(* (r^2) is NOT (x^2)+(y^2)+(z^2) *)
(*r=(Sqrt[(1/2)((x^2)+(y^2)+(z^2)-(a^2)+(Sqrt[((((x^2)+(y^2)+(z^2)-(a^2)))^2)+4 (a^2) (z^2)]))]);
R=(Sqrt[(1/2)((x^2)+(y^2)+(z^2)+(a^2)+(Sqrt[((((x^2)+(y^2)+(z^2)-(a^2)))^2)+4 (a^2) (z^2)]))]);
ρ=((((((x^2)+(y^2)+(z^2)-(a^2)))^2)+4 (a^2) (z^2))^(1/4));*)
R[r_]:=(Sqrt[(r^2)+(a^2)]);
ρ[r_,z_]:=(Sqrt[(r^2)+(((a^2)(z^2))/(r^2))]);
Derivative[1,0][rtr][x_,y_]:=(x/(rtr[x,y]));
Derivative[0,1][rtr][x_,y_]:=(y/(rtr[x,y]));
Derivative[1,0][r][rtr_,z_]:=rtr ((r[rtr,z])/((ρ[r[rtr,z],z])^2));
Derivative[0,1][r][rtr_,z_]:=((z ((R[r[rtr,z]])^2))/(r[rtr,z]((ρ[r[rtr,z],z])^2)));
Derivative[1][ζ][r_]:=((a β[r])/((R[r])^2));
(* (Subscript[α, m]) *)
α[0]=1;
α[1]=((y a)/((R[r[rtr[x,y],z]])^2));
α[2]=-((x a)/((R[r[rtr[x,y],z]])^2));
α[3]=0;
(* (β^μ) *)
β[0]=0;
Do[β[μ]=β[r[rtr[x,y],z]](∂_(Coordinate[μ]))r[rtr[x,y],z],{μ,1,3}];
InverseVierbeinMatrix=Table[KroneckerDelta[i,j]+α[i]β[j],{i,0,3},{j,0,3}];
,TimeOption===Time&&SpaceOption===Polar,
R[r_]:=(Sqrt[(r^2)+(a^2)]);
ρ[r_,θ_]:=(Sqrt[(r^2)+(a^2)((Cos[θ])^2)]);
InverseVierbeinMatrix=(({{1, ((((R[r])^2) β[r])/((ρ[r,θ])^2)), 0, 0}, {0, ((R[r])/(ρ[r,θ])), 0, 0}, {0, 0, (1/(ρ[r,θ])), 0}, {0, -((a R[r] Sin[θ] β[r])/((ρ[r,θ])^2)), 0, (1/(R[r]Sin[θ]))}}));
];

(* Arnowitt-Deser-Misner (3+1) *)
,VierbeinOption===ArnowittDeserMisner,
Clear[α,β];
InverseVierbeinMatrix=Table[Which[m==0&&μ==0,(1/(FunctionCoordinates[α])),m==0,((FunctionCoordinates[β[μ]])/(FunctionCoordinates[α])),True,KroneckerDelta[m,μ]],{m,FirstIndex,LastIndex},{μ,FirstIndex,LastIndex}];

(* Dupin cyclide *)
,VierbeinOption===DupinCyclide,
If[SpaceOption===Cartesian,
Clear[fxf,fyf,fzf,fwf,fx,fy,fz,fw,ηf,ξf,ρf,xf,yf,zf,wf,η,ξ,ρ,fDupinCyclide];
Clear[rx,ry,rz,rw];
Clear[x0,y0,z0,w0,v0];
Switch[LastIndex
(* 2 spatial dimensions *)
,2,
ξf=ξ[x,y];
fxf=fx;fyf=fy[x];
xf=(1/rx);yf=(1/(ry+rx Cos[x]));
ρf=(((rx Sin[x]-x0))^2)+((((ry+rx Cos[x])Sin[y]-y0))^2)+((((ry+rx Cos[x])Cos[y]-z0))^2);
(* {x,y,z} ∈ [0,2π]; [π,π,π] is at centre of geometry, [0,0,0] at antipodes *)
fDupinCyclide={ξf->Log[ρf],(∂_({x,ix_Integer},{y,iy_Integer}))ξf:>(∂_({x,ix},{y,iy}))Log[Evaluate[ρf]],fxf->xf,fyf->yf,(∂_({x,ix_Integer}))fyf:>(∂_({x,ix}))yf};
InverseVierbeinMatrix=Switch[FirstIndex
(* time and space *)
,0,(({{1, 0, 0}, {0, (E^(ξf))fxf, 0}, {0, 0, (E^(ξf))fyf}}))
(* just space, no time *)
,1,(({{(E^(ξf))fxf, 0}, {0, (E^(ξf))fyf}}))
]
(* 3 spatial dimensions *)
,3,
ξf=ξ[x,y,z];
fxf=fx;fyf=fy[x];fzf=fz[x,y];
xf=(1/rx);yf=(1/(ry+rx Cos[x]));zf=(1/(rz+(ry+rx Cos[x])Cos[y]));
ρf=(((rx Sin[x]-x0))^2)+((((ry+rx Cos[x])Sin[y]-y0))^2)+((((rz+(ry+rx Cos[x])Cos[y])Sin[z]-z0))^2)+((((rz+(ry+rx Cos[x])Cos[y])Cos[z]-w0))^2);
(* {x,y,z} ∈ [0,2π]; [π,π,π] is at centre of geometry, [0,0,0] at antipodes *)
fDupinCyclide={ξf->Log[ρf],(∂_({x,ix_Integer},{y,iy_Integer},{z,iz_Integer}))ξf:>(∂_({x,ix},{y,iy},{z,iz}))Log[Evaluate[ρf]],fxf->xf,fyf->yf,(∂_({x,ix_Integer}))fyf:>(∂_({x,ix}))yf,fzf->zf,(∂_({x,ix_Integer},{y,iy_Integer}))fzf:>(∂_({x,ix},{y,iy}))zf};
InverseVierbeinMatrix=Switch[FirstIndex
(* time and space *)
,0,
(({{1, 0, 0, 0}, {0, (E^(ξf))fxf, 0, 0}, {0, 0, (E^(ξf))fyf, 0}, {0, 0, 0, (E^(ξf))fzf}}))
(* just space, no time *)
,1,(({{(E^(ξf))fxf, 0, 0}, {0, (E^(ξf))fyf, 0}, {0, 0, (E^(ξf))fzf}}))
]
(* 4 spatial dimensions *)
,4,
ξf=ξ[x,y,z,w];
fxf=fx;fyf=fy[x];fzf=fz[x,y];fwf=fw[x,y,z];
xf=(1/rx);yf=(1/(ry+rx Cos[x]));zf=(1/(rz+(ry+rx Cos[x])Cos[y]));wf=(1/(rw+(rz+(ry+rx Cos[x])Cos[y])Cos[z]));
ρf=(((rx Sin[x]-x0))^2)+((((ry+rx Cos[x])Sin[y]-y0))^2)+((((rz+(ry+rx Cos[x])Cos[y])Sin[z]-z0))^2)+((((rw+(rz+(ry+rx Cos[x])Cos[y])Cos[z])Sin[w]-w0))^2)+((((rw+(rz+(ry+rx Cos[x])Cos[y])Cos[z])Cos[w]-v0))^2);
(* {x,y,z,w} ∈ [0,2π]; [π,π,π,π] is at centre of geometry, [0,0,0,0] at antipodes *)
fDupinCyclide={ξf->Log[ρf],(∂_({x,ix_Integer},{y,iy_Integer},{z,iz_Integer},{w,iw_Integer}))ξf:>(∂_({x,ix},{y,iy},{z,iz},{w,iw}))Log[Evaluate[ρf]],fxf->xf,fyf->yf,(∂_({x,ix_Integer}))fyf:>(∂_({x,ix}))yf,fzf->zf,(∂_({x,ix_Integer},{y,iy_Integer}))fzf:>(∂_({x,ix},{y,iy}))zf,fwf->wf,(∂_({x,ix_Integer},{y,iy_Integer},{z,iz_Integer}))fwf:>(∂_({x,ix},{y,iy},{z,iz}))wf};
InverseVierbeinMatrix=Switch[FirstIndex
(* time and space *)
,0,
(({{1, 0, 0, 0, 0}, {0, (E^(ξf))fxf, 0, 0, 0}, {0, 0, (E^(ξf))fyf, 0, 0}, {0, 0, 0, (E^(ξf))fzf, 0}, {0, 0, 0, 0, (E^(ξf))fwf}}))
(* just space, no time *)
,1,(({{(E^(ξf))fxf, 0, 0, 0}, {0, (E^(ξf))fyf, 0, 0}, {0, 0, (E^(ξf))fzf, 0}, {0, 0, 0, (E^(ξf))fwf}}))
]
];
];
];

(* Re-define simplification rules for printing out quantities *)
Which[
VierbeinOption===Doran,
Clear[SPrint];
SPrint[f__]:=Simplify[f//.{rtr[x,y]->rtr,r[rtr,z]->r,(x^(m_Integer))/;EvenQ[m]->(((((rtr[x,y])^2)-(y^2)))^(m/2)),(rtr^(n_Integer))/;EvenQ[n]->(((((((r^2)+(a^2))((r^2)-(z^2)))/(r^2))))^(n/2))}];
]

(* transform vierbein to desired tetrad metric *)
InverseVierbeinMatrix=TetradTransformMatrix.InverseVierbeinMatrix;

Clear[InverseVierbein];
InverseVierbein[m_Integer,μ_Integer]:=InverseVierbein[m,μ]=InverseVierbeinMatrix[[m+1-FirstIndex,μ+1-FirstIndex]];
Table[SPrint[InverseVierbein[m,μ]],{m,FirstIndex,LastIndex},{μ,FirstIndex,LastIndex}]//MatrixForm

(* InverseVierbein determinant |((Subscript[e, m])^μ)| *)
Clear[DeterminantInverseVierbein];
DeterminantInverseVierbein=S[Det[InverseVierbeinMatrix]];
SPrint[DeterminantInverseVierbein]

If[Perturb,
Clear[UnperturbedInverseVierbeinMatrix,UnperturbedInverseVierbein,UnperturbedDeterminantInverseVierbein];
UnperturbedInverseVierbeinMatrix=InverseVierbeinMatrix;
UnperturbedInverseVierbein[m_Integer,μ_Integer]:=InverseVierbein[m,μ];
UnperturbedDeterminantInverseVierbein=DeterminantInverseVierbein;
]

(* Explicit expression for (β^i) in Doran-Cartesian tetrad *)
If[VierbeinOption===Doran,
If[TimeOption===Time&&SpaceOption===Cartesian,
Switch[TetradMetricOption,
Cartesian,
Module[{temp},Do[temp=SPrint[S[β[i]-(∂_(Coordinate[i]))(Integrate[β[r],r]/.r->r[rtr[x,y],z])]];If[temp=!=0,Print["Difference[",i,"] = ",temp]],{i,FirstIndex,LastIndex}]]
,_,Null
]
]
]

Compute the vierbein tensor

(* Vierbein tensor (Subscript[(e^m), μ]) *)
Clear[VierbeinMatrix];
VierbeinMatrix=S[Transpose[Inverse[InverseVierbeinMatrix]]];

Clear[Vierbein];
Vierbein[m_Integer,μ_Integer]:=Vierbein[m,μ]=VierbeinMatrix[[m+1-FirstIndex,μ+1-FirstIndex]];
Table[SPrint[Vierbein[m,μ]],{m,FirstIndex,LastIndex},{μ,FirstIndex,LastIndex}]//MatrixForm

If[Perturb,
Clear[UnperturbedVierbeinMatrix,UnperturbedVierbein];
UnperturbedVierbeinMatrix=VierbeinMatrix;
UnperturbedVierbein[m_Integer,μ_Integer]:=Vierbein[m,μ];
]

Choose the perturbation of the vierbein

If[Perturb,
(* Options for perturbation *)
Clear[NoGauge,Newtonian,ScalarNewtonian];
(* Choose perturbation *)
Clear[GaugeOption];
GaugeOption=Newtonian
]

(* InverseVierbein perturbation matrix (Subscript[φ, mn]) *)
If[Perturb,
Clear[CovariantPerturbationMatrix];

Which[
GaugeOption===NoGauge,
Clear[φ,φf];
φf[i_Integer,j_Integer]:=FunctionCoordinates[φ[i,j]];
CovariantPerturbationMatrix=ϵ Table[φf[i,j],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}];

,GaugeOption===Newtonian,
Clear[hc,hp,wy,wz,Φ,Ψ];
Clear[hcf,hpf,Wf,Φf,Ψf,hijf];
Φf=Φ[Coordinate[0],Coordinate[3]];
Ψf=Ψ[Coordinate[0],Coordinate[3]];
Wf[3]:=0;
Wf[i_Integer]:=W[i][Coordinate[0],Coordinate[3]];
(* hp = (SubPlus[h]) = (Subscript[h, xx]) = -(Subscript[h, yy]) *)
hpf=hp[Coordinate[0],Coordinate[3]];
(* hc = (Subscript[h, ×]) = (Subscript[h, xy]) = (Subscript[h, yx]) *)
hcf=hc[Coordinate[0],Coordinate[3]];
hijf[i_Integer,j_Integer]:=hijf[j,i]/;i>j;
hijf[1,1]=hpf;
hijf[2,2]=-hpf;
hijf[1,2]=hcf;
CovariantPerturbationMatrix=ϵ(({{Ψf, Wf[1], Wf[2], 0}, {0, Φf+hpf, hcf, 0}, {0, hcf, Φf-hpf, 0}, {0, 0, 0, Φf}}));

,GaugeOption===ScalarNewtonian,
Clear[Φ,Ψ];
Clear[Φf,Ψf];
Φf=Φ[Coordinate[0],Coordinate[3]];
Ψf=Ψ[Coordinate[0],Coordinate[3]];
CovariantPerturbationMatrix=ϵ(({{Ψf, 0, 0, 0}, {0, Φf, 0, 0}, {0, 0, Φf, 0}, {0, 0, 0, Φf}}));
];
]

CovariantPerturbationMatrix=S[TetradTransformMatrix.CovariantPerturbationMatrix.InverseTetradTransformMatrix];

If[Perturb,
Clear[CovariantPerturbation];
CovariantPerturbation[m_Integer,μ_Integer]:=CovariantPerturbation[m,μ]=CovariantPerturbationMatrix[[m+1-FirstIndex,μ+1-FirstIndex]];
Table[SPrint[CovariantPerturbation[m,μ]],{m,FirstIndex,LastIndex},{μ,FirstIndex,LastIndex}]//MatrixForm
]

If[Perturb,
Clear[UnitMatrix];
UnitMatrix=Table[KroneckerDelta[i,j],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}];
InverseVierbeinMatrix=(UnitMatrix+CovariantPerturbationMatrix.InverseTetradMetricMatrix).UnperturbedInverseVierbeinMatrix;
VierbeinMatrix=(UnitMatrix-InverseTetradMetricMatrix.(CovariantPerturbationMatrix^\[Transpose])).UnperturbedVierbeinMatrix;
]

If[Perturb,
Clear[InverseVierbein];
InverseVierbein[m_Integer,μ_Integer]:=InverseVierbein[m,μ]=InverseVierbeinMatrix[[m+1-FirstIndex,μ+1-FirstIndex]];
Table[SPrint[InverseVierbein[m,μ]],{m,FirstIndex,LastIndex},{μ,FirstIndex,LastIndex}]//MatrixForm
]

If[Perturb,
Clear[Vierbein];
Vierbein[m_Integer,μ_Integer]:=Vierbein[m,μ]=VierbeinMatrix[[m+1-FirstIndex,μ+1-FirstIndex]];
Table[SPrint[Vierbein[m,μ]],{m,FirstIndex,LastIndex},{μ,FirstIndex,LastIndex}]//MatrixForm
]

(* should be unit matrix *)
If[Perturb,
S[(VierbeinMatrix^\[Transpose]).InverseVierbeinMatrix]//MatrixForm
]

Compute the metric tensor

(* Metric tensor (Subscript[g, μν]) = (Subscript[γ, mn])(Subscript[(e^m), μ])(Subscript[(e^n), ν]) *)
Clear[Metric];
Metric[μ_Integer,ν_Integer]:=Metric[ν,μ]/;μ>ν;
Metric[μ_Integer,ν_Integer]:=Metric[μ,ν]=Module[{m,n},S[(∑_(m=FirstIndex)^LastIndex)(∑_(n=FirstIndex)^LastIndex)TetradMetric[m,n]Vierbein[m,μ]Vierbein[n,ν]]];
Table[SPrint[Metric[μ,ν]],{μ,FirstIndex,LastIndex},{ν,FirstIndex,LastIndex}]//MatrixForm

(* Inverse metric tensor (g^(μν)) = (γ^mn)((Subscript[e, m])^μ)((Subscript[e, n])^ν) *)
(*Clear[InverseMetric];
InverseMetric[μ_Integer,ν_Integer]:=InverseMetric[ν,μ]/;μ>ν;
InverseMetric[μ_Integer,ν_Integer]:=InverseMetric[μ,ν]=Module[{m,n},S[(∑_(m=FirstIndex)^LastIndex)(∑_(n=FirstIndex)^LastIndex)InverseTetradMetric[m,n]InverseVierbein[m,μ]InverseVierbein[n,ν]]];
Table[SPrint[InverseMetric[μ,ν]],{μ,FirstIndex,LastIndex},{ν,FirstIndex,LastIndex}]//MatrixForm*)

Define the directed derivative

(* Directed derivative (∂_m) = ((Subscript[e, m])^μ)(∂/(∂(x^μ))) *)
Clear[Ddirected];
Ddirected[m_,x__]:=Module[{μ},S[(∑_(μ=FirstIndex)^LastIndex)InverseVierbein[m,μ](∂_(Coordinate[μ]))x]];

If[Perturb,
Clear[UnperturbedDdirected];
UnperturbedDdirected[m_,x__]:=Module[{μ},S[(∑_(μ=FirstIndex)^LastIndex)UnperturbedInverseVierbein[m,μ](∂_(Coordinate[μ]))x]];
]

Compute vierbein derivatives

Clear[DVierbein,CovariantDVierbein];
(* Vierbein derivatives (Subscript[(d^k), mn]) = -(Subscript[(e^k), κ])(∂_n)((Subscript[e, m])^κ) = ((Subscript[e, m])^μ)(∂_n)(Subscript[(e^k), μ]) *)
Module[{k,l,m,n,κ},
(*DVierbein[k_,m_,n_]:=DVierbein[k,m,n]=S[-(∑_(κ=FirstIndex)^LastIndex)Vierbein[k,κ]Ddirected[n,InverseVierbein[m,κ]]];*)
DVierbein[k_,m_,n_]:=DVierbein[k,m,n]=S[(∑_(μ=FirstIndex)^LastIndex)InverseVierbein[m,μ]Ddirected[n,Vierbein[k,μ]]];
(* Covariant vierbein derivatives (Subscript[d, lmn]) = (Subscript[γ, lk])(Subscript[(d^k), mn]) *)
CovariantDVierbein[l_,m_,n_]:=CovariantDVierbein[l,m,n]=S[(∑_(k=FirstIndex)^LastIndex)TetradMetric[l,k]DVierbein[k,m,n]];
];

(*Module[{temp},Do[temp=SPrint[CovariantDVierbein[i,j,k]];If[temp=!=0,Print["CovariantDVierbein[",i,",",j,",",k,"] = ",temp]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex},{k,FirstIndex,LastIndex}]]*)

Compute the connection coefficients (all indices down)

(* Connection coefficients (Subscript[Γ, ijk]) *)
Clear[CovariantConnection];
If[TetradMetricIsConstant,
CovariantConnection[i_,i_,k_]:=0;
CovariantConnection[i_Integer,j_Integer,k_]:=-CovariantConnection[j,i,k]/;i>j;
];
CovariantConnection[i_Integer,j_Integer,k_Integer]:=CovariantConnection[i,j,k]=S[(1/2)(Ddirected[k,TetradMetric[i,j]]+Ddirected[j,TetradMetric[i,k]]-Ddirected[i,TetradMetric[j,k]]-(CovariantDVierbein[i,j,k]-CovariantDVierbein[i,k,j]+CovariantDVierbein[j,k,i]- CovariantDVierbein[j,i,k]+CovariantDVierbein[k,j,i]-CovariantDVierbein[k,i,j]))];

Module[{temp},Do[temp=SPrint[CovariantConnection[i,j,k]];If[temp=!=0,Print["CovariantConnection[",i,",",j,",",k,"] = ",temp]],{i,FirstIndex,LastIndex},{j,If[TetradMetricIsConstant,i,FirstIndex],LastIndex},{k,FirstIndex,LastIndex}]]

(* Explicit solution for tetrad connections in perturbed FLRW in Newtonian or PlaneWave gauge *)
If[Perturb&&VierbeinOption===FriedmannLemaitreRobertsonWalker&&TimeOption===Conformal&&SpaceOption===Cartesian,
Clear[Ff];
Switch[GaugeOption
,Newtonian,
Ff=(((∂_η)a[η])/(a[η]))Ψf+(∂_η)Φf;
Do[Print[0,i,0," : ",S[CovariantConnection[0,i,0]-(ϵ/(a[η]))(-(∂_(Coordinate[i]))Ψf)]],{i,1,LastIndex}];
Do[Print[0,i,j," : ",S[CovariantConnection[0,i,j]-(-(((∂_η)a[η])/((a[η])^2))KroneckerDelta[i,j]+(ϵ/(a[η]))(Ff KroneckerDelta[i,j]-(1/2)((∂_(Coordinate[i]))Wf[j]+(∂_(Coordinate[j]))Wf[i])+(∂_η)hijf[i,j]))]],{i,1,LastIndex},{j,1,LastIndex}];
Do[Print[i,j,0," : ",S[CovariantConnection[i,j,0]-(ϵ/(a[η]))((1/2)((∂_(Coordinate[i]))Wf[j]-(∂_(Coordinate[j]))Wf[i]))]],{i,1,LastIndex},{j,1,LastIndex}];
Do[Print[i,j,0," : ",S[CovariantConnection[i,j,k]-(ϵ/(a[η]))(KroneckerDelta[j,k](∂_(Coordinate[i]))Φf-KroneckerDelta[i,k](∂_(Coordinate[j]))Φf+(∂_(Coordinate[i]))hijf[j,k]-(∂_(Coordinate[j]))hijf[i,k])]],{i,1,LastIndex},{j,1,LastIndex},{k,1,LastIndex}];
];
]

(* Explicit solution for tetrad connections for general spherical metric *)
If[VierbeinOption===GeneralSpherical,
Which[FirstIndex==0&&TimeOption===Time&&(SpaceOption===Polar||SpaceOption===ComovingPolar),
Switch[TetradMetricOption
,Cartesian,
Clear[β,g,h];
β[0]=Ddirected[0,r];
β[1]=Ddirected[1,r];
Switch[SpaceOption
,Polar,
g=Ddirected[1,Log[α[t,r]]];
h=(∂_r)β[0]+β[0](∂_r)Log[α[t,r]]-Ddirected[0,Log[β[1]]];
,ComovingPolar,
g=Ddirected[1,Log[α[t,x]]];
h=(∂_x)λt[t,x]+λt[t,x](∂_x)Log[α[t,x]]-Ddirected[0,Log[λx[t,x]]];
];
Module[{temp},Do[temp=SPrint[S[CovariantConnection[i,j,k]+MetricSign Switch[{i,j,k},{0,1,0},g,{0,1,1},h,{0,2,2}|{0,3,3},((β[0])/r),{1,2,2}|{1,3,3},((β[1])/r),{2,3,3},((Cot[θ])/r),_,0]]];If[temp=!=0,Print["Difference[",i,",",j,",",k,"] = ",temp]],{i,FirstIndex,LastIndex},{j,If[TetradMetricIsConstant,i,FirstIndex],LastIndex},{k,FirstIndex,LastIndex}]];
,NewmanPenrose,
Clear[β,g,h];
β[0]=Ddirected[0,r];
β[1]=Ddirected[1,r];
Switch[SpaceOption
,Polar,
g=βr[t,r](∂_r)Log[α[t,r]];
h=(∂_r)βt[t,r]+βt[t,r](∂_r)Log[α[t,r]]-((1/(α[t,r]))(∂_t)Log[βr[t,r]]+βt[t,r](∂_r)Log[βr[t,r]]);
,ComovingPolar,
g=((Ddirected[0,Log[α[t,x]]]-Ddirected[1,Log[α[t,x]]])/(Sqrt[2]));
h=(∂_x)λt[t,x]+λt[t,x](∂_x)Log[α[t,x]]-((Ddirected[0,Log[λx[t,x]]]+Ddirected[1,Log[λx[t,x]]])/(Sqrt[2]));
];
Module[{temp},Do[temp=SPrint[S[CovariantConnection[i,j,k]+MetricSign Switch[{i,j,k},{0,1,0},-((g+h)/(Sqrt[2])),{0,1,1},-((g-h)/(Sqrt[2])),{0,2,3}|{0,3,2},((β[0])/r),{1,2,3}|{1,3,2},((β[1])/r),{2,3,2},((Cot[θ])/((Sqrt[2])r)),{2,3,3},-((Cot[θ])/((Sqrt[2])r)),_,0]]];If[temp=!=0,Print["Difference[",i,",",j,",",k,"] = ",temp]],{i,FirstIndex,LastIndex},{j,If[TetradMetricIsConstant,i,FirstIndex],LastIndex},{k,FirstIndex,LastIndex}]],
_,Null
];
,FirstIndex==0&&TimeOption===Conformal&&SpaceOption===ComovingPolar,
Switch[TetradMetricOption
,Cartesian,
Clear[β,g,h];
β[0]=Ddirected[0,A];
β[1]=Ddirected[1,A];
g=A Ddirected[1,Log[A ξη[x]]];
h=A Ddirected[0,Log[A ξx[x]]];
Module[{temp},Do[temp=SPrint[S[CovariantConnection[i,j,k]+MetricSign(1/A)Switch[{i,j,k},{0,1,0},g,{0,1,1},h,{0,2,2}|{0,3,3},β[0],{1,2,2}|{1,3,3},β[1],{2,3,3},Cot[θ],_,0]]];If[temp=!=0,Print["Difference[",i,",",j,",",k,"] = ",temp]],{i,FirstIndex,LastIndex},{j,If[TetradMetricIsConstant,i,FirstIndex],LastIndex},{k,FirstIndex,LastIndex}]];
,_,Null
]
]
]

(* Explicit solution for tetrad connections for Doran-Cartesian tetrad *)
If[VierbeinOption===Doran,
If[TimeOption===Time&&SpaceOption===Cartesian,
Switch[TetradMetricOption,
Minkowski,
Module[{temp},Do[temp=SPrint[S[CovariantConnection[i,j,k]+MetricSign(∂_(Coordinate[k]))(α[i]β[j]-α[j]β[i]+Signature[{0,i,j,3}]ζ[r[rtr[x,y],z]])]];If[temp=!=0,Print["Difference[",i,",",j,",",k,"] = ",temp]],{i,FirstIndex,LastIndex},{j,If[TetradMetricIsConstant,i,FirstIndex],LastIndex},{k,FirstIndex,LastIndex}]]
,_,Null
]
]
]

(* Explicit solution for tetrad connections in Arnowitt-Deser-Misner formalism *)
If[VierbeinOption===ArnowittDeserMisner,
Module[{temp},
Clear[Gravity,ExtrinsicCurvature];
Gravity[0]=0;
Gravity[i_]:=-CovariantDVierbein[0,0,i];
ExtrinsicCurvature[i_,0]=0;
ExtrinsicCurvature[0,j_]=0;
ExtrinsicCurvature[i_Integer,j_Integer]:=ExtrinsicCurvature[j,i]/;i>j;
ExtrinsicCurvature[i_,j_]:=ExtrinsicCurvature[i,j]=S[(1/2)(Ddirected[0,TetradMetric[i,j]]-(-CovariantDVierbein[i,j,0]-CovariantDVierbein[j,i,0]+CovariantDVierbein[i,0,j]+CovariantDVierbein[j,0,i]))];
Do[temp=SPrint[S[CovariantConnection[i,j,k]-Which[i==0&&j==0,0,i==0&&k==0,-Gravity[j],j==0&&k==0,Gravity[i],k==0,ExtrinsicCurvature[i,j]-CovariantDVierbein[i,j,0]+CovariantDVierbein[i,0,j],j==0,ExtrinsicCurvature[i,k],i==0,-ExtrinsicCurvature[j,k],True,(1/2)Switch[TetradMetricOption
,Minkowski,-(CovariantDVierbein[i,j,k]-CovariantDVierbein[i,k,j]+CovariantDVierbein[j,k,i]-CovariantDVierbein[j,i,k]+CovariantDVierbein[k,j,i]-CovariantDVierbein[k,i,j])
,Spatial,(∂_(Coordinate[k]))TetradMetric[i,j]+(∂_(Coordinate[j]))TetradMetric[i,k]-(∂_(Coordinate[i]))TetradMetric[j,k]]]]];If[temp=!=0,Print["Difference[",i,",",j,",",k,"] = ",temp]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex},{k,FirstIndex,LastIndex}];
];
]

Coordinate gauge transformation of the inverse vierbein

If[Perturb,
(* (Subscript[φ, ij]) -> (Subscript[φ, ij]) - (∂_i)(Subscript[ϵ, j]) + ((Subscript[Γ, kji]) + (Subscript[Γ, jik]))(ε^k) *)
Clear[CovariantCoordinateShift,CovariantCoordinateShiftf,CovariantVierbeinShift];
CovariantCoordinateShiftf[i_Integer]:=ϵ CovariantCoordinateShift[i][Coordinate[0],Coordinate[1],Coordinate[2],Coordinate[3]];
CovariantVierbeinShift[i_Integer,j_Integer]:=CovariantVierbeinShift[i,j]=Module[{k},S[-UnperturbedDdirected[i,CovariantCoordinateShiftf[j]]+(∑_(k=FirstIndex)^LastIndex)(((CovariantConnection[k,j,i]/.ϵ->0)+(CovariantConnection[j,i,k]/.ϵ->0))CovariantCoordinateShiftf[k])]]
]

If[Perturb,
Do[Print["CovariantVierbeinShift[",i,",",j,"]= ",SPrint[CovariantVierbeinShift[i,j]]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}]
]

Compute the components of the completely covariant Riemann tensor (all indices down).  This uses all the symmetries.

(* Covariant Riemann tensor (Subscript[R, ijkl]) *)
Clear[CovariantRiemannf];
CovariantRiemannf[i_Integer,j_Integer,k_Integer,l_Integer]:=Module[{m,n},S[Ddirected[i,CovariantConnection[k,l,j]]-Ddirected[j,CovariantConnection[k,l,i]]+(∑_(m=FirstIndex)^LastIndex)(∑_(n=FirstIndex)^LastIndex)InverseTetradMetric[m,n](CovariantConnection[m,k,j]CovariantConnection[n,l,i]-CovariantConnection[m,k,i]CovariantConnection[n,l,j]+(CovariantConnection[m,i,j]-CovariantConnection[m,j,i])CovariantConnection[k,l,n])]];

(* Covariant Riemann tensor invoking symmetries *)
Clear[CovariantRiemann];
CovariantRiemann[i_,i_,k_,l_]:=0;
CovariantRiemann[i_,j_,k_,k_]:=0;
CovariantRiemann[i_Integer,j_Integer,k_,l_]:=-CovariantRiemann[j,i,k,l]/;i>j;
CovariantRiemann[i_,j_,k_Integer,l_Integer]:=-CovariantRiemann[i,j,l,k]/;k>l;
CovariantRiemann[i_Integer,j_Integer,k_Integer,l_Integer]:=CovariantRiemann[k,l,i,j]/;i<j&&k<l&&(i>k||i==k&&j>l);
CovariantRiemann[i_Integer,j_Integer,k_Integer,l_Integer]:=-CovariantRiemann[i,l,j,k]-CovariantRiemann[i,k,l,j]/;i<k<l<j;
CovariantRiemann[i_Integer,j_Integer,k_Integer,l_Integer]:=CovariantRiemann[i,j,k,l]=CovariantRiemannf[i,j,k,l];

Module[{temp},Do[If[(i<j&&k<l&&(i<k||(i==k&&j≤l))),temp=SPrint[CovariantRiemann[i,j,k,l]];If[temp=!=0,Print["CovariantRiemann[",i,",",j,",",k,",",l,"] = ",temp]]],{i,FirstIndex,LastIndex},{j,i,LastIndex},{k,i,LastIndex},{l,k,LastIndex}]]

(* Explicit solution for Riemann tensor for general spherical metric *)
If[VierbeinOption===GeneralSpherical,
If[FirstIndex==0&&((TimeOption===Time&&(SpaceOption===Polar||SpaceOption===ComovingPolar))||(TimeOption===Conformal&&SpaceOption===ComovingPolar)),
Switch[TetradMetricOption,
Cartesian,
Clear[M,R0101,R0202,R1212,R0212,R1202];
Which[TimeOption===Time&&(SpaceOption===Polar||SpaceOption===ComovingPolar),
M=S[(r/2)(1+((β[0])^2)-((β[1])^2))];
R0101=S[Ddirected[0,h]-Ddirected[1,g]+(h^2)- (g^2)];
R0202=S[(1/r)(Ddirected[0,β[0]]-β[1]g)];
R1212=S[(1/r)(Ddirected[1,β[1]]-β[0]h)];
R0212=S[(1/r)(Ddirected[0,β[1]]-β[0]g)];
R1202=S[(1/r)(Ddirected[1,β[0]]-β[1]h)];
R2323=-((2M)/(r^3));
,TimeOption===Conformal&&SpaceOption===ComovingPolar,
M=S[(1/2)(1+((β[0])^2)-((β[1])^2))];
R0101=S[(A^2) (Ddirected[0,(h/A)]-Ddirected[1,(g/A)])+(h^2)-(g^2)];
R0202=S[A Ddirected[0,β[0]]-β[1]g];
R1212=S[A Ddirected[1,β[1]]-β[0]h];
R0212=S[A Ddirected[0,β[1]]-β[0]g];
R1202=S[A Ddirected[1,β[0]]-β[1]h];
R2323=-2M;
];
Module[{temp},Do[If[(i<j&&k<l&&(i<k||(i==k&&j≤l))),temp=SPrint[S[CovariantRiemann[i,j,k,l]+MetricSign(1/(A^2))Switch[{i,j,k,l},{0,1,0,1},R0101,{0,2,0,2}|{0,3,0,3},R0202,{1,2,1,2}|{1,3,1,3},R1212,{0,2,1,2}|{0,3,1,3},R0212,{2,3,2,3},R2323,_,0]]];If[temp=!=0,Print["Difference[",i,",",j,",",k,",",l,"] = ",temp]]],{i,FirstIndex,LastIndex},{j,i,LastIndex},{k,i,LastIndex},{l,k,LastIndex}]]
,_,Null
]
]
]

(* Explicit solution for Riemann tensor in Arnowitt-Deser-Misner formalism *)
If[VierbeinOption===ArnowittDeserMisner,
(* (Subscript[D, i]^((3)))(Subscript[D, j]^((3)))α *)
Clear[DDCovariantSpaceLapse];
Module[{k,l},
DDCovariantSpaceLapse[i_Integer,j_Integer]:=DDCovariantSpaceLapse[j,i]/;i>j;
DDCovariantSpaceLapse[i_,j_]:=DDCovariantSpaceLapse[i,j]=-MetricSign(1/(FunctionCoordinates[α]))((∂_(Coordinate[i]))(∂_(Coordinate[j]))FunctionCoordinates[α]-(∑_(k=1)^LastIndex)(∑_(l=1)^LastIndex)InverseTetradMetric[k,l]CovariantConnection[k,j,i](∂_(Coordinate[l]))FunctionCoordinates[α])];
(* (Subscript[D, t])(Subscript[K, ij]) *)
Clear[DCovariantTimeExtrinsicCurvature];
Module[{k,l},
DCovariantTimeExtrinsicCurvature[i_Integer,j_Integer]:=DCovariantTimeExtrinsicCurvature[j,i]/;i>j;
DCovariantTimeExtrinsicCurvature[i_,j_]:=DCovariantTimeExtrinsicCurvature[i,j]=S[Ddirected[0,ExtrinsicCurvature[i,j]]-(∑_(k=FirstIndex)^LastIndex)(∑_(l=FirstIndex)^LastIndex)InverseTetradMetric[k,l](CovariantConnection[k,i,0]ExtrinsicCurvature[j,l]+CovariantConnection[k,j,0]ExtrinsicCurvature[i,l])]];
(* (Subscript[D, i]^((3)))(Subscript[K, jk]) *)
Clear[DCovariantSpaceExtrinsicCurvature];
Module[{l,m},
DCovariantSpaceExtrinsicCurvature[i_,j_Integer,k_Integer]:=DCovariantSpaceExtrinsicCurvature[i,k,j]/;j>k;
DCovariantSpaceExtrinsicCurvature[i_,j_,k_]:=DCovariantSpaceExtrinsicCurvature[i,j,k]=S[(∂_(Coordinate[i]))ExtrinsicCurvature[j,k]-(∑_(l=1)^LastIndex)(∑_(m=1)^LastIndex)InverseTetradMetric[l,m](CovariantConnection[l,j,i]ExtrinsicCurvature[k,m]+CovariantConnection[l,k,i]ExtrinsicCurvature[j,m])]];
(* (Subscript[R, ijkl]^((3))) *)
Clear[CovariantRiemannSpace,CovariantRiemannSpacef];Module[{m,n},CovariantRiemannSpacef[i_Integer,j_Integer,k_Integer,l_Integer]:=S[(∂_(Coordinate[i]))CovariantConnection[k,l,j]-(∂_(Coordinate[j]))CovariantConnection[k,l,i]+(∑_(m=1)^LastIndex)(∑_(n=1)^LastIndex)InverseTetradMetric[m,n](CovariantConnection[m,k,j]CovariantConnection[n,l,i]-CovariantConnection[m,k,i]CovariantConnection[n,l,j]+(CovariantConnection[m,i,j]-CovariantConnection[m,j,i])CovariantConnection[k,l,n])]];
CovariantRiemannSpace[i_,j_,k_,l_]:=0/;i===0||j===0||k===0||l===0;
CovariantRiemannSpace[i_,i_,k_,l_]:=0;
CovariantRiemannSpace[i_,j_,k_,k_]:=0;
CovariantRiemannSpace[i_Integer,j_Integer,k_,l_]:=-CovariantRiemannSpace[j,i,k,l]/;i>j;
CovariantRiemannSpace[i_,j_,k_Integer,l_Integer]:=-CovariantRiemannSpace[i,j,l,k]/;k>l;
CovariantRiemannSpace[i_Integer,j_Integer,k_Integer,l_Integer]:=CovariantRiemannSpace[k,l,i,j]/;i<j&&k<l&&(i>k||i==k&&j>l);
CovariantRiemannSpace[i_Integer,j_Integer,k_Integer,l_Integer]:=-CovariantRiemannSpace[i,l,j,k]-CovariantRiemannSpace[i,k,l,j]/;i<k<l<j;
CovariantRiemannSpace[i_Integer,j_Integer,k_Integer,l_Integer]:=CovariantRiemannSpace[i,j,k,l]=CovariantRiemannSpacef[i,j,k,l];
Module[{m,n,temp},
Do[If[(i<j&&k<l&&(i<k||(i==k&&j≤l))),temp=SPrint[S[CovariantRiemann[i,j,k,l]-Which[
i==0&&k==0,
-DCovariantTimeExtrinsicCurvature[j,l]-(∑_(m=1)^LastIndex)(∑_(n=1)^LastIndex)InverseTetradMetric[m,n]ExtrinsicCurvature[j,m]ExtrinsicCurvature[l,n]-DDCovariantSpaceLapse[j,l]
(*-Ddirected[0,ExtrinsicCurvature[j,l]+-(∑_(m=1)^LastIndex)(∑_(n=1)^LastIndex)InverseTetradMetric[m,n](ExtrinsicCurvature[l,m]CovariantDVierbein[n,0,j]+ExtrinsicCurvature[j,m]CovariantDVierbein[n,0,l])+(∑_(m=1)^LastIndex)(∑_(n=1)^LastIndex)InverseTetradMetric[m,n]ExtrinsicCurvature[j,m]ExtrinsicCurvature[l,n]-DDCovariantSpaceLapse[j,l]*)
,i==0,
DCovariantSpaceExtrinsicCurvature[l,j,k]-DCovariantSpaceExtrinsicCurvature[k,j,l]
,True,
-MetricSign(ExtrinsicCurvature[i,l]ExtrinsicCurvature[j,k]-ExtrinsicCurvature[i,k]ExtrinsicCurvature[j,l])+CovariantRiemannSpace[i,j,k,l]
]]];
If[True,Print["Difference[",i,",",j,",",k,",",l,"] = ",temp]]],{i,FirstIndex,LastIndex},{j,i,LastIndex},{k,FirstIndex,LastIndex},{l,k,LastIndex}]];
]

The Ricci tensor

(* Ricci tensor (Subscript[R, ij]) *)
Clear[CovariantRicci];
CovariantRicci[i_Integer,j_Integer]:=CovariantRicci[j,i]/;i>j;
CovariantRicci[i_Integer,j_Integer]:=CovariantRicci[i,j]=Module[{k,l},S[(∑_(k=FirstIndex)^LastIndex)(∑_(l=FirstIndex)^LastIndex)InverseTetradMetric[k,l]CovariantRiemann[k,i,l,j]]];

(*Module[{temp},Do[temp=SPrint[CovariantRicci[i,j]];If[temp=!=0,Print["CovariantRicci[",i,",",j,"] = ",temp]],{i,FirstIndex,LastIndex},{j,i,LastIndex}]]*)

(* Explicit solution for Ricci tensor for general spherical metric *)
If[VierbeinOption===GeneralSpherical,
If[FirstIndex==0&&((TimeOption===Time&&(SpaceOption===Polar||SpaceOption===ComovingPolar))||(TimeOption===Conformal&&SpaceOption===ComovingPolar)),
Switch[TetradMetricOption,
Cartesian,
Clear[F,Falt,R,Ralt,P,Palt,Pperp,Pperpalt];
Which[TimeOption===Time&&(SpaceOption===Polar||SpaceOption===ComovingPolar),
(* (Subscript[R, tθrθ]) *)
(*F=R0212;*)
F=S[(1/r)(Ddirected[0,β[1]]-β[0]g)];
(* (Subscript[R, rθtθ]) *)
(*Falt=R1202;*)
Falt=S[(1/r)(Ddirected[1,β[0]]-β[1]h)];
(* - (Subscript[R, rθrθ]) - (1/2)(Subscript[R, θϕθϕ]) *)
(*R=S[-R1212-(1/2)R2323];*)
R=S[(1/r)(-Ddirected[1,β[1]]+β[0]h+(M/(r^2)))];
Ralt=S[(1/(β[1]))((1/(r^2))Ddirected[1,M]-β[0]F)];
(* - (Subscript[R, tθtθ]) + (1/2)(Subscript[R, θϕθϕ]) *)
(*P=S[-R0202+(1/2)R2323];*)
P=S[(1/r)(-Ddirected[0,β[0]]+β[1]g-(M/(r^2)))];
Palt=S[-(1/(β[0]))((1/(r^2))Ddirected[0,M]+β[1]F)];
(* - (Subscript[R, trtr]) - (Subscript[R, tθtθ]) + (Subscript[R, rθrθ]) *)
Pperp=S[(1/2)(-(1/r)Ddirected[0,r h+β[0]]+(1/r)Ddirected[1,r g+β[1]]-(h^2)+(g^2))];
Pperpalt[1]=S[-(r/(2β[0]))(Ddirected[0,R]+Ddirected[1,F]+(2((β[0])/r)+h)R+2(((β[1])/r)+g)F+h P)];
Pperpalt[2]=S[(r/(2β[1]))(Ddirected[1,P]+Ddirected[0,F]+(2((β[1])/r)+g)P+2(((β[0])/r)+h)F+g R)];
,TimeOption===Conformal&&SpaceOption===ComovingPolar,
(* (A^2)(Subscript[R, tθrθ]) *)
F=S[A Ddirected[0,β[1]]-β[0]g];
(* (A^2)(Subscript[R, rθtθ]) *)
Falt=S[A Ddirected[1,β[0]]-β[1]h];
(* (A^2)(- (Subscript[R, rθrθ]) - (1/2)(Subscript[R, θϕθϕ])) *)
R=S[-A Ddirected[1,β[1]]+β[0]h+M];
Ralt=S[(1/(β[1]))(Ddirected[1,A M]-β[0]F)];
(* (A^2)(- (Subscript[R, tθtθ]) + (1/2)(Subscript[R, θϕθϕ])) *)
P=S[-A Ddirected[0,β[0]]+β[1]g-M];
Palt=S[-(1/(β[0]))(Ddirected[0,A M]+β[1]F)];
(* (A^2)(- (Subscript[R, trtr]) - (Subscript[R, tθtθ]) + (Subscript[R, rθrθ])) *)
Pperp=S[(1/2)(A Ddirected[1,g+β[1]]-A Ddirected[0,h+β[0]]+(g^2)-(h^2))];
Pperpalt[1]=S[-(1/(2β[0]))(A Ddirected[0,R]+A Ddirected[1,F]+h(R+P)+2g F)];
Pperpalt[2]=S[(1/(2β[1]))(A Ddirected[1,P]+A Ddirected[0,F]+g(R+P)+2h F)];
];
Module[{temp},Do[temp=SPrint[S[CovariantRicci[i,j]-(1/(A^2))(Switch[{i,j},{0,0},R+P+2Pperp,{0,1},-2F,{1,1},R+P-2Pperp,{2,2}|{3,3},R-P,_,0])]];If[temp=!=0,Print["Difference[",i,",",j,"] = ",temp]],{i,FirstIndex,LastIndex},{j,i,LastIndex}]];
Module[{temp},Do[temp=SPrint[S[CovariantRicci[i,j]-(1/(A^2))(Switch[{i,j},{0,0},Ralt+Palt+2Pperpalt[1],{0,1},-2Falt,{1,1},Ralt+Palt-2Pperpalt[2],{2,2}|{3,3},Ralt-Palt,_,0])]];If[temp=!=0,Print["Differencealt[",i,",",j,"] = ",temp]],{i,FirstIndex,LastIndex},{j,i,LastIndex}]]
,_,Null
]
]
]

(* Explicit solution for Ricci tensor in Arnowitt-Deser-Misner formalism *)
If[VierbeinOption===ArnowittDeserMisner,
(* K *)
Clear[ExtrinsicCurvatureScalar];
Module[{i,j},ExtrinsicCurvatureScalar=S[(∑_(i=1)^LastIndex)(∑_(j=1)^LastIndex)InverseTetradMetric[i,j]ExtrinsicCurvature[i,j]]];
(* (Subscript[R, ij]^((3))) *)
Clear[CovariantRicciSpace];
Module[{k,l},
CovariantRicciSpace[i_,j_]:=0/;i===0||j===0;
CovariantRicciSpace[i_Integer,j_Integer]:=CovariantRicciSpace[j,i]/;i>j;
CovariantRicciSpace[i_Integer,j_Integer]:=CovariantRicciSpace[i,j]=S[(∑_(k=1)^LastIndex)(∑_(l=1)^LastIndex)InverseTetradMetric[k,l]CovariantRiemannSpace[k,i,l,j]]];
Module[{k,l,m,n,temp},
Do[temp=SPrint[S[CovariantRicci[i,j]-Which[
i==0&&j==0,
-Ddirected[0,ExtrinsicCurvatureScalar]-(∑_(k=1)^LastIndex)(∑_(l=1)^LastIndex)InverseTetradMetric[k,l]((∑_(m=1)^LastIndex)(∑_(n=1)^LastIndex)InverseTetradMetric[m,n]ExtrinsicCurvature[k,m]ExtrinsicCurvature[l,n]+DDCovariantSpaceLapse[k,l])
,i==0,
(∑_(k=1)^LastIndex)(∑_(l=1)^LastIndex)InverseTetradMetric[k,l]DCovariantSpaceExtrinsicCurvature[k,l,j]-(∂_(Coordinate[j]))ExtrinsicCurvatureScalar
,True,
-DCovariantTimeExtrinsicCurvature[i,j]-ExtrinsicCurvatureScalar ExtrinsicCurvature[i,j]-DDCovariantSpaceLapse[i,j]+CovariantRicciSpace[i,j]
]]];
If[True,Print["Difference[",i,",",j,"] = ",temp]],{i,FirstIndex,LastIndex},{j,i,LastIndex}]];
]

The Ricci scalar

(* Ricci scalar R *)
Clear[RicciScalar];
RicciScalar=Module[{i,j},S[(∑_(i=FirstIndex)^LastIndex)(∑_(j=FirstIndex)^LastIndex)InverseTetradMetric[i,j]CovariantRicci[i,j]]];
SPrint[RicciScalar]

(* Explicit solution for Ricci scalar for general spherical metric *)
If[VierbeinOption===GeneralSpherical,
If[FirstIndex==0&&((TimeOption===Time&&(SpaceOption===Polar||SpaceOption===ComovingPolar))||(TimeOption===Conformal&&SpaceOption===ComovingPolar)),
Switch[TetradMetricOption,
Cartesian,
Clear[RScalar];
RScalar=2(-R+P+2Pperp);
Module[{temp},temp=SPrint[S[RicciScalar+MetricSign(1/(A^2))RScalar]];If[temp=!=0,Print["Difference = ",temp]]]
,_,Null
]
]
]

(* Explicit solution for Ricci scalar in Arnowitt-Deser-Misner formalism *)
If[VierbeinOption===ArnowittDeserMisner,
(* (R^((3))) *)
Clear[RicciScalarSpace];
Module[{i,j},
RicciScalarSpace:=RicciScalarSpace=S[(∑_(i=1)^LastIndex)(∑_(j=1)^LastIndex)InverseTetradMetric[i,j]CovariantRicciSpace[i,j]]];
Module[{k,l,m,n,temp},
temp=SPrint[S[RicciScalar-
-2Ddirected[0,ExtrinsicCurvatureScalar]-(ExtrinsicCurvatureScalar^2)-(∑_(k=1)^LastIndex)(∑_(l=1)^LastIndex)InverseTetradMetric[k,l]((∑_(m=1)^LastIndex)(∑_(n=1)^LastIndex)InverseTetradMetric[m,n]ExtrinsicCurvature[k,m]ExtrinsicCurvature[l,n]+2DDCovariantSpaceLapse[k,l])
+RicciScalarSpace]];
If[True,Print["Difference] = ",temp]]];
]

Compute the covariant Einstein tensor (both indices down)

(* Einstein tensor (Subscript[G, ij]) *)
Clear[CovariantEinstein];
CovariantEinstein[i_Integer,j_Integer]:=CovariantEinstein[j,i]/;i>j;
CovariantEinstein[i_Integer,j_Integer]:=CovariantEinstein[i,j]=S[CovariantRicci[i,j]-(1/2)TetradMetric[i,j]RicciScalar]

Module[{temp},Do[temp=SPrint[CovariantEinstein[i,j]];If[temp=!=0,Print["CovariantEinstein[",i,",",j,"] = ",temp]],{i,FirstIndex,LastIndex},{j,i,LastIndex}]]

(* Explicit solution for tetrad connections in perturbed FLRW in Newtonian gauge *)
If[Perturb&&VierbeinOption===FriedmannLemaitreRobertsonWalker&&TimeOption===Conformal&&SpaceOption===Cartesian,
Switch[GaugeOption
,Newtonian,
Print[0,0," : ",S[CovariantEinstein[0,0]-(3((((((∂_η)a[η])/((a[η])^2))))^2)+(ϵ/((a[η])^2))(-6(((∂_η)a[η])/(a[η]))Ff+2(∑_(i=1)^3)(∑_(j=1)^3)InverseTetradMetric[i,j](∂_(Coordinate[i],Coordinate[j]))Φf))]];
Do[Print[0,i," : ",S[CovariantEinstein[0,i]-(ϵ/((a[η])^2))(2(∂_(Coordinate[i]))Ff+(1/2)(∑_(j=1)^3)(∑_(k=1)^3)InverseTetradMetric[j,k](∂_(Coordinate[j],Coordinate[k]))Wf[i])]],{i,1,LastIndex}];
Do[Print[i,j," : ",S[CovariantEinstein[i,j]-((-2(((∂_(η,η))a[η])/((a[η])^3))+((((((∂_η)a[η])/((a[η])^2))))^2))KroneckerDelta[i,j]+(ϵ/((a[η])^2))((2((∂_η)Ff+2(((∂_η)a[η])/(a[η]))Ff)+2((((∂_(η,η))a[η])/(a[η]))-2((((((∂_η)a[η])/(a[η]))))^2))Ψf)KroneckerDelta[i,j]-((∂_(Coordinate[i],Coordinate[j]))(Ψf-Φf)-(∑_(k=1)^3)(∑_(l=1)^3)InverseTetradMetric[k,l](∂_(Coordinate[k],Coordinate[l]))(Ψf-Φf)KroneckerDelta[i,j])+(1/2)(∂_η)((∂_(Coordinate[i]))Wf[j]+(∂_(Coordinate[j]))Wf[i])+(((∂_η)a[η])/(a[η]))((∂_(Coordinate[i]))Wf[j]+(∂_(Coordinate[j]))Wf[i])-((∂_(η,η))hijf[i,j]+2(((∂_η)a[η])/(a[η]))(∂_η)hijf[i,j]-(∑_(k=1)^3)(∑_(l=1)^3)InverseTetradMetric[k,l](∂_(Coordinate[k],Coordinate[l]))hijf[i,j])))]],{i,1,LastIndex},{j,1,LastIndex}];
]]

(* Explicit solution for Einstein tensor for general spherical metric *)
If[VierbeinOption===GeneralSpherical,
If[FirstIndex==0&&((TimeOption===Time&&(SpaceOption===Polar||SpaceOption===ComovingPolar))||(TimeOption===Conformal&&SpaceOption===ComovingPolar)),
Switch[TetradMetricOption,
Cartesian,
Module[{temp},Do[temp=SPrint[S[CovariantEinstein[i,j]-(2/(A^2))Switch[{i,j},{0,0},R,{0,1},-F,{1,1},P,{2,2}|{3,3},Pperp,_,0]]];If[temp=!=0,Print["Difference[",i,",",j,"] = ",temp]],{i,FirstIndex,LastIndex},{j,i,LastIndex}]];
Module[{temp},Do[temp=SPrint[S[CovariantEinstein[i,j]-(2/(A^2))Switch[{i,j},{0,0},Ralt,{0,1},-Falt,{1,1},Palt,{2,2},Pperpalt[1],{3,3},Pperpalt[2],_,0]]];If[temp=!=0,Print["Differencealt[",i,",",j,"] = ",temp]],{i,FirstIndex,LastIndex},{j,i,LastIndex}]]
,_,Null
]
]
]

(* Explicit solution for Dupin cyclide *)
If[ValueQ[fDupinCyclide],
If[FirstIndex==0,
Print[Simplify[CovariantEinstein[0,0]-(1/2)RicciScalar]];Print[Simplify[(LastIndex-2)CovariantEinstein[0,0]+(∑_(i=1)^LastIndex)CovariantEinstein[i,i]]];
];
(* the Einstein tensor of a torus is a sum of products of these factors *)
Clear[ftorusx,ftorusy,ftorusz,ftorusw];
If[LastIndex≥1,ftorusx=(1/rx)];
If[LastIndex≥2,ftorusy=((Cos[x])/(ry+rx Cos[x]))];
If[LastIndex≥3,ftorusz=((Cos[x]Cos[y])/(rz+(ry+rx Cos[x])Cos[y]))];
If[LastIndex≥4,ftorusw=((Cos[x]Cos[y]Cos[z])/(rw+(rz+(ry+rx Cos[x])Cos[y])Cos[z]))];
(* the Einstein tensor of a Dupin cyclide is a sum of products of these factors *)
Clear[fDupinx,fDupiny,fDupinz,fDupinw];
If[LastIndex≥1,fDupinx=Simplify[ftorusx(ρf/.Solve[(1/ftorusx)==0,rx][[1]])-(1/ftorusx)]];
If[LastIndex≥2,fDupiny=Simplify[ftorusy(ρf/.Solve[(1/ftorusy)==0,rx][[1]])-(1/ftorusy)]];
If[LastIndex≥3,fDupinz=Simplify[ftorusz(ρf/.Solve[(1/ftorusz)==0,rx][[1]])-(1/ftorusz)]];
If[LastIndex≥4,fDupinw=Simplify[ftorusw(ρf/.Solve[(1/ftorusw)==0,rx][[1]])-(1/ftorusw)]];
Switch[LastIndex
(* 2 spatial dimensions *)
,2,
If[FirstIndex==0,
Simplify[(CovariantEinstein[0,0]/.fDupinCyclide)-fDupinx fDupiny]
]
(* 3 spatial dimensions *)
,3,
Module[{temp},Do[temp=SPrint[(CovariantEinstein[i,j]/.fDupinCyclide)+Switch[{i,j},
{1,1},fDupiny fDupinz,
{2,2},fDupinx fDupinz,
{3,3},fDupinx fDupiny,
_,0]];If[True,Print[Subscript["G",i,j],": ",temp]],{i,1,LastIndex},{j,i,LastIndex}]]
(* 4 spatial dimensions *)
,4,
Module[{temp},Do[temp=SPrint[(CovariantEinstein[i,j]/.fDupinCyclide)+Switch[{i,j},
{1,1},fDupiny fDupinz+fDupiny fDupinw+fDupinz fDupinw,
{2,2},fDupinx fDupinz+fDupinx fDupinw+fDupinz fDupinw,
{3,3},fDupinx fDupiny+fDupinx fDupinw+fDupiny fDupinw,
{4,4},fDupinx fDupiny+fDupinx fDupinz+fDupiny fDupinz,
_,0]];If[True,Print[Subscript["G",i,j],": ",temp]],{i,1,LastIndex},{j,i,LastIndex}]]
]
]

Compute the completely covariant Weyl tensor (all indices down)

(* Minus the Ricci part of the Riemann tensor: add this to Riemann to get Weyl. *)
Clear[CovariantRicciCombo];
CovariantRicciCombo[k_,k_,m_,n_]:=0;
CovariantRicciCombo[k_,l_,m_,m_]:=0;
CovariantRicciCombo[k_Integer,l_Integer,m_,n_]:=-CovariantRicciCombo[l,k,m,n]/;k>l;
CovariantRicciCombo[k_,l_,m_Integer,n_Integer]:=-CovariantRicciCombo[k,l,n,m]/;m>n;
CovariantRicciCombo[k_Integer,l_Integer,m_Integer,n_Integer]:=CovariantRicciCombo[m,n,k,l]/;k<l&&m<n&&(k>m||(k==m&&l>n));
CovariantRicciCombo[k_Integer,l_Integer,m_Integer,n_Integer]:=-CovariantRicciCombo[k,n,l,m]-CovariantRicciCombo[k,m,n,l]/;k<m<n<l;
CovariantRicciCombo[k_Integer,l_Integer,m_Integer,n_Integer]:=CovariantRicciCombo[k,l,m,n]=S[-(1/2)(TetradMetric[k,m]CovariantRicci[l,n]-TetradMetric[k,n]CovariantRicci[l,m]+TetradMetric[l,n]CovariantRicci[k,m]-TetradMetric[l,m]CovariantRicci[k,n])+(1/6)(TetradMetric[k,m]TetradMetric[l,n]-TetradMetric[k,n]TetradMetric[l,m])RicciScalar];

(* Weyl tensor (Subscript[C, klmn]) *)
Clear[CovariantWeyl];
CovariantWeyl[k_,k_,m_,n_]:=0;
CovariantWeyl[k_,l_,m_,m_]:=0;
CovariantWeyl[k_Integer,l_Integer,m_,n_]:=-CovariantWeyl[l,k,m,n]/;k>l;
CovariantWeyl[k_,l_,m_Integer,n_Integer]:=-CovariantWeyl[k,l,n,m]/;m>n;
CovariantWeyl[k_Integer,l_Integer,m_Integer,n_Integer]:=CovariantWeyl[m,n,k,l]/;k<l&&m<n&&(k>m||(k==m&&l>n));
CovariantWeyl[k_Integer,l_Integer,m_Integer,n_Integer]:=-CovariantWeyl[k,n,l,m]-CovariantWeyl[k,m,n,l]/;k<m<n<l;
CovariantWeyl[k_Integer,l_Integer,m_Integer,n_Integer]:=CovariantWeyl[k,l,m,n]=S[CovariantRiemann[k,l,m,n]+CovariantRicciCombo[k,l,m,n]];

Do[If[(i<j&&k<l&&(i<k||(i==k&&j≤l))),temp=SPrint[CovariantWeyl[i,j,k,l]];If[temp=!=0,Print["CovariantWeyl[",i,",",j,",",k,",",l,"] = ",temp]]],{i,FirstIndex,LastIndex},{j,i,LastIndex},{k,i,LastIndex},{l,k,LastIndex}]

(* Explicit solution for Weyl tensor for general spherical metric *)
If[VierbeinOption===GeneralSpherical,
If[FirstIndex==0&&((TimeOption===Time&&(SpaceOption===Polar||SpaceOption===ComovingPolar))||(TimeOption===Conformal&&SpaceOption===ComovingPolar)),
Switch[TetradMetricOption,
Cartesian,
Clear[WeylScalar];
Which[TimeOption===Time&&(SpaceOption===Polar||SpaceOption===ComovingPolar),
WeylScalar=(1/3)(R-P+Pperp-((3M)/(r^3))),
TimeOption===Conformal&&SpaceOption===ComovingPolar,
WeylScalar=(1/3)(R-P+Pperp-3M)];
Module[{temp},Do[If[(i<j&&k<l&&(i<k||(i==k&&j≤l))),temp=SPrint[S[CovariantWeyl[i,j,k,l]+MetricSign(1/(A^2))Switch[{i,j,k,l},{0,1,0,1},-2,{0,2,0,2}|{0,3,0,3},1,{1,2,1,2}|{1,3,1,3},-1,{2,3,2,3},2,_,0]WeylScalar]];If[temp=!=0,Print["Difference[",i,",",j,",",k,",",l,"] = ",temp]]],{i,FirstIndex,LastIndex},{j,i,LastIndex},{k,i,LastIndex},{l,k,LastIndex}]]
,_,Null
]
]
]

Define complexification projection operator

(* (1/2)(((Subscript[δ, i])^k)((Subscript[δ, j])^l) + (I/2)((Subscript[ε, ij])^kl))
The factor of (1/2) makes this a projection operator P, satisfying (P^2)=P. *)
Clear[ComplexificationProjectionOperator];
ComplexificationProjectionOperator[i_Integer,j_Integer,k_Integer,l_Integer]:=ComplexificationProjectionOperator[i,j,k,l]=Module[{m,n},S[(1/2)(KroneckerDelta[i,k]KroneckerDelta[j,l]+(I/2)(∑_(m=FirstIndex)^LastIndex)InverseTetradMetric[k,m](∑_(n=FirstIndex)^LastIndex)InverseTetradMetric[l,n]TotallyAntisymmetric[i,j,m,n])]]

(*Do[temp=SPrint[ComplexificationProjectionOperator[i,j,k,l]];If[temp=!=0,Print["ComplexificationProjectionOperator[",i,",",j,",",k,",",l,"] = ",temp]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex},{k,FirstIndex,LastIndex},{l,FirstIndex,LastIndex}]*)

Complexified Weyl tensor

Clear[ComplexifiedCovariantWeyl];
ComplexifiedCovariantWeyl[i_,i_,k_,l_]:=0;
ComplexifiedCovariantWeyl[i_,j_,k_,k_]:=0;
ComplexifiedCovariantWeyl[i_Integer,j_Integer,k_,l_]:=-ComplexifiedCovariantWeyl[j,i,k,l]/;i>j;
ComplexifiedCovariantWeyl[i_,j_,k_Integer,l_Integer]:=-ComplexifiedCovariantWeyl[i,j,l,k]/;k>l;
ComplexifiedCovariantWeyl[i_Integer,j_Integer,k_Integer,l_Integer]:=ComplexifiedCovariantWeyl[k,l,i,j]/;i<j&&k<l&&(i>k||i==k&&j>l);
ComplexifiedCovariantWeyl[i_Integer,j_Integer,k_Integer,l_Integer]:=-ComplexifiedCovariantWeyl[i,l,j,k]-ComplexifiedCovariantWeyl[i,k,l,j]/;i<k<l<j;
ComplexifiedCovariantWeyl[i_,j_,k_,l_]:=ComplexifiedCovariantWeyl[i,j,k,l]=Module[{m,n,p,q},S[(∑_(m=FirstIndex)^LastIndex)(∑_(n=FirstIndex)^LastIndex)ComplexificationProjectionOperator[i,j,m,n](∑_(p=FirstIndex)^LastIndex)(∑_(q=FirstIndex)^LastIndex)ComplexificationProjectionOperator[k,l,p,q]CovariantWeyl[m,n,p,q]]];

Do[If[i<j&&k<l&&(i<k||(i==k&&j≤l)),temp=SPrint[ComplexifiedCovariantWeyl[i,j,k,l]];If[temp=!=0,Print["ComplexifiedCovariantWeyl[",i,",",j,",",k,",",l,"] = ",temp]]],{i,FirstIndex,LastIndex},{j,i,LastIndex},{k,i,LastIndex},{l,k,LastIndex}]

(* Explicit solution for complexified Weyl tensor for general spherical metric *)
If[VierbeinOption===GeneralSpherical,
If[FirstIndex==0&&((TimeOption===Time&&(SpaceOption===Polar||SpaceOption===ComovingPolar))||(TimeOption===Conformal&&SpaceOption===ComovingPolar)),
Switch[TetradMetricOption,
Cartesian,
Module[{temp},
Do[If[(i<j&&k<l&&(i<k||(i==k&&j≤l))),temp=SPrint[S[ComplexifiedCovariantWeyl[i,j,k,l]+MetricSign(1/(2(A^2)))Switch[{i,j,k,l},{0,1,0,1},-2,{0,2,0,2}|{0,3,0,3},1,{1,2,1,2}|{1,3,1,3},-1,{2,3,2,3},2,{0,1,2,3},-2I,{0,2,1,3},-I,{0,3,1,2},I,_,0]WeylScalar]];If[temp=!=0,Print["Difference[",i,",",j,",",k,",",l,"] = ",temp]]],{i,FirstIndex,LastIndex},{j,i,LastIndex},{k,i,LastIndex},{l,k,LastIndex}]]
,_,Null
]
]
]

Electromagnetic field tensor

Clear[ElectricFieldVector,MagneticFieldVector];
Which[VierbeinOption===GeneralSpherical,
Which[FirstIndex===0&&((TimeOption===Time&&(SpaceOption===Polar||SpaceOption===ComovingPolar))||(TimeOption===Conformal&&SpaceOption===ComovingPolar)),
Switch[TetradMetricOption,
Minkowski,
Which[TimeOption===Time&&SpaceOption===Polar,
ElectricFieldVector={((Q[t,r])/(r^2)),0,0};
,TimeOption===Time&&SpaceOption===ComovingPolar,
ElectricFieldVector={((Q[t,x])/(r^2)),0,0};
,TimeOption===Conformal&&SpaceOption===ComovingPolar,
ElectricFieldVector={((Q[x])/A),0,0};
];
MagneticFieldVector={0,0,0};
,_,Null
]
];
,VierbeinOption===BoyerLindquist,
Which[TimeOption===Time&&SpaceOption===Polar,
(* (Subscript[E, r]) + I (Subscript[B, r]) = (Q/(((r - I a Cos[θ]))^2)) *)
ElectricFieldVector={Q(((r^2)-(a^2)((Cos[θ])^2))/((ρ[r,θ])^4)),0,0};
MagneticFieldVector=Switch[TetradMetricOption,NewmanPenrose,I,_,1]{Q((2a r Cos[θ])/((ρ[r,θ])^4)),0,0};
];
]

If[ValueQ[ElectricFieldVector]&&ValueQ[MagneticFieldVector],
Clear[CovariantElectromagneticFieldMatrix];
CovariantElectromagneticFieldMatrix=(({{0, ElectricFieldVector[[1]], ElectricFieldVector[[2]], ElectricFieldVector[[3]]}, {-ElectricFieldVector[[1]], 0, -MagneticFieldVector[[3]], MagneticFieldVector[[2]]}, {-ElectricFieldVector[[2]], MagneticFieldVector[[3]], 0, -MagneticFieldVector[[1]]}, {-ElectricFieldVector[[3]], -MagneticFieldVector[[2]], MagneticFieldVector[[1]], 0}}));
]

If[ValueQ[CovariantElectromagneticFieldMatrix],
Clear[CovariantElectromagneticField];
CovariantElectromagneticField[i_Integer,j_Integer]:=CovariantElectromagneticFieldMatrix[[i+1-FirstIndex,j+1-FirstIndex]];
Table[SPrint[CovariantElectromagneticField[i,j]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}]//MatrixForm
]

If[ValueQ[CovariantElectromagneticFieldMatrix],
Clear[CovariantElectromagneticEnergyMomentum];
Module[{k,l,m,n},CovariantElectromagneticEnergyMomentum[i_Integer,j_Integer]:=CovariantElectromagneticEnergyMomentum[i,j]=S[(MetricSign/(4π))((∑_(k=FirstIndex)^LastIndex)(∑_(l=FirstIndex)^LastIndex)InverseTetradMetric[k,l]CovariantElectromagneticField[i,k]CovariantElectromagneticField[j,l]-(1/4)InverseTetradMetric[i,j](∑_(k=FirstIndex)^LastIndex)(∑_(l=FirstIndex)^LastIndex)(∑_(m=FirstIndex)^LastIndex)(∑_(n=FirstIndex)^LastIndex)InverseTetradMetric[k,l]InverseTetradMetric[m,n]CovariantElectromagneticField[k,m]CovariantElectromagneticField[l,n])]];
]

If[ValueQ[CovariantElectromagneticFieldMatrix],
Do[temp=SPrint[CovariantElectromagneticEnergyMomentum[i,j]];If[temp=!=0,Print["CovariantElectromagneticEnergyMomentum[",i,",",j,"] = ",temp]],{i,FirstIndex,LastIndex},{j,i,LastIndex}]
]

If[ValueQ[CovariantElectromagneticFieldMatrix],
Clear[CovariantElectricCurrent];
CovariantElectricCurrent[i_Integer]:=Module[{j,k,l,m,n},S[(1/(4π))(∑_(j=FirstIndex)^LastIndex)(∑_(k=FirstIndex)^LastIndex)InverseTetradMetric[j,k](Ddirected[j,CovariantElectromagneticField[k,i]]-(∑_(l=FirstIndex)^LastIndex)(∑_(m=FirstIndex)^LastIndex)InverseTetradMetric[l,m](CovariantConnection[l,k,j]CovariantElectromagneticField[m,i]+CovariantConnection[l,i,j]CovariantElectromagneticField[k,m]))]];
]

If[ValueQ[CovariantElectromagneticFieldMatrix],
Do[temp=SPrint[CovariantElectricCurrent[i]];If[temp=!=0,Print["CovariantElectricCurrent[",i,"] = ",temp]],{i,0,3}]
]

If[ValueQ[CovariantElectromagneticFieldMatrix],
(* FIX: overall sign? *)
Clear[CovariantMagneticCurrent];
Module[{j,k,l,m,n},CovariantMagneticCurrent[i_Integer]:=S[(∑_(j=FirstIndex)^LastIndex)(∑_(k=FirstIndex)^LastIndex)(∑_(l=FirstIndex)^LastIndex)Signature[{i,j,k,l}](Ddirected[j,CovariantElectromagneticField[k,l]]-(∑_(m=FirstIndex)^LastIndex)(∑_(n=FirstIndex)^LastIndex)InverseTetradMetric[m,n](CovariantConnection[m,k,j]CovariantElectromagneticField[n,l]+CovariantConnection[m,l,j]CovariantElectromagneticField[k,n]))]];
]

If[ValueQ[CovariantElectromagneticFieldMatrix],
Do[temp=SPrint[CovariantMagneticCurrent[i]];If[temp=!=0,Print["CovariantElectromagneticEnergyMomentum[",i,",",j,"] = ",temp]],{i,FirstIndex,LastIndex}]
]

(* Explicit solution for electrogmagnetic current vector for general spherical metric *)
If[VierbeinOption===GeneralSpherical,
Which[FirstIndex===0&&((TimeOption===Time&&(SpaceOption===Polar||SpaceOption===ComovingPolar))||(TimeOption===Conformal&&SpaceOption===ComovingPolar)),
Clear[jf,Qf];
Which[TimeOption===Time&&SpaceOption===Polar,
jf=(1/(4π (r^2)));
Qf=Q[t,r];
,TimeOption===Time&&SpaceOption===ComovingPolar,
jf=(1/(4π (r^2)));
Qf=Q[t,x];
,TimeOption===Conformal&&SpaceOption===ComovingPolar,
jf=(1/(4π (A^2)));
Qf=A Q[x];
];
Switch[TetradMetricOption,
Minkowski,
Module[{i},Do[temp=SPrint[S[CovariantElectricCurrent[i]+MetricSign jf Switch[i,0,Ddirected[1,Qf],1,Ddirected[0,Qf],_,0]]];If[temp=!=0,Print["Difference[",i,"] = ",temp]],{i,FirstIndex,LastIndex}]]
,_,Null
]
]
]

Complexified electromagnetic field tensor

Clear[ComplexifiedCovariantElectromagneticField];
ComplexifiedCovariantElectromagneticField[i_,i_]:=0;
ComplexifiedCovariantElectromagneticField[i_Integer,j_Integer]:=-ComplexifiedCovariantElectromagneticField[j,i]/;i>j;
ComplexifiedCovariantElectromagneticField[i_,j_]:=ComplexifiedCovariantElectromagneticField[i,j]=Module[{k,l,m,n},S[(∑_(k=FirstIndex)^LastIndex)(∑_(l=FirstIndex)^LastIndex)ComplexificationProjectionOperator[i,j,k,l]CovariantElectromagneticField[k,l]]];

If[ValueQ[CovariantElectromagneticFieldMatrix],
Table[SPrint[ComplexifiedCovariantElectromagneticField[i,j]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}]//MatrixForm
]

Clear[ComplexifiedCovariantElectromagneticCurrent];
ComplexifiedCovariantElectromagneticCurrent[i_Integer]:=Module[{j,k,l,m,n},S[(1/(4π))(∑_(j=FirstIndex)^LastIndex)(∑_(k=FirstIndex)^LastIndex)InverseTetradMetric[j,k](Ddirected[j,ComplexifiedCovariantElectromagneticField[k,i]]-(∑_(l=FirstIndex)^LastIndex)(∑_(m=FirstIndex)^LastIndex)InverseTetradMetric[l,m](CovariantConnection[l,k,j]ComplexifiedCovariantElectromagneticField[m,i]+CovariantConnection[l,i,j]ComplexifiedCovariantElectromagneticField[k,m]))]];

If[ValueQ[CovariantElectromagneticFieldMatrix],
Do[temp=SPrint[ComplexifiedCovariantElectromagneticCurrent[i]];If[temp=!=0,Print["ComplexifiedCovariantElectromagneticCurrent[",i,"] = ",temp]],{i,0,3}]
]

(* Explicit solution for electrogmagnetic current vector for general spherical metric *)
If[VierbeinOption===GeneralSpherical,
Which[FirstIndex==0&&((TimeOption===Time&&(SpaceOption===Polar||SpaceOption===ComovingPolar))||(TimeOption===Conformal&&SpaceOption===ComovingPolar)),
Switch[TetradMetricOption,
Cartesian,
Module[{i},Do[temp=SPrint[S[ComplexifiedCovariantElectromagneticCurrent[i]+MetricSign(jf/2)Switch[i,0,Ddirected[1,Qf],1,Ddirected[0,Qf],_,0]]];If[temp=!=0,Print["Difference[",i,"] = ",temp]],{i,FirstIndex,LastIndex}]]
,_,Null
]
]
]