(* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||| Mathematica Syntax | EINSTEIN-MAXWELL TENSOR+GEODESIC SOLVER | geodesics.yukterez.net |||| *)
(* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

ClearAll["Local`*"];
smp[y_]:=Simplify[y,r>0&&-M<a<M&&M>0];
list[y_]:=y[[1]]==y[[2]];
rplc[y_]:=y;

(* kovariante metrische Komponenten *)
g11=gtt=1-((2 M r)/((r^2)+(a^2) ((Cos[θ])^2)));
g22=grr=-(((a^2)+2 (r^2)+(a^2) Cos[2 θ])/(2 ((a^2)+r (-2 M+r))));
g33=gθθ=-(r^2)-(a^2) ((Cos[θ])^2);
g44=gφφ=-(((a^2)+(r^2)) ((Sin[θ])^2))-((2 (a^2) M r ((Sin[θ])^4))/((r^2)+(a^2) ((Cos[θ])^2)));
g12=gtr=0;
g13=gtθ=0;
g14=gtφ=-((2 a M r ((Sin[θ])^2))/((r^2)+(a^2) ((Cos[θ])^2)));
g23=grθ=0;
g24=grφ=0;
g34=gθφ=0;

(* Dimensionen *)
x={t, r, θ, φ}; n=4; 

"Metrischer Tensor"
mt={
{g11, g12, g13, g14},
{g12, g22, g23, g24},
{g13, g23, g33, g34},
{g14, g24, g34, g44}
};
Subscript["g", μσ] -> MatrixForm[mt]
it=smp[Inverse[mt]];
"g"^μσ -> MatrixForm[it]
mx=Table[smp[Sum[
it[[i, k]] mt[[k, j]], {k, 1, n}]],
{i, 1, n}, {j, 1, n}];
Subsuperscript["g", "μ", "σ"] -> MatrixForm[mx]
md=Det[mt]; "|g|" -> smp[md]

"Elektromagnetisches Vektorpotential"
A={0,0,0,0}; "A"^μ->smp[A]

"Maxwell Tensor"
F=ParallelTable[smp[((D[A[[j]], x[[k]]]-D[A[[k]], x[[j]]]))], {j, 1, n}, {k, 1, n}];
Subscript["F", μσ] -> MatrixForm[F]
f=smp[ParallelTable[Sum[
it[[i, k]] it[[j, l]] F[[k, l]],
{k, 1, n}, {l, 1, n}], {i, 1, n}, {j, 1, n}]];
"F"^μσ -> MatrixForm[f]
џ=ParallelTable[smp[Sum[
it[[i, k]] F[[k, j]], {k, 1, n}]],
{i, 1, n}, {j, 1, n}];
Subsuperscript["F", "μ", "σ"] -> MatrixForm[џ]

"Christoffelsymbole"
chr=ParallelTable[smp[(1/2)Sum[(it[[i, s]])
(D[mt[[s, j]], x[[k]]]+D[mt[[s, k]], x[[j]]] -D[mt[[j, k]], x[[s]]]), {s, 1, n}]],
{i, 1, n}, {j, 1, n}, {k, 1, n}];
crs=ParallelTable[If[UnsameQ[chr[[i, j, k]], 0],
{Subsuperscript["Γ", ToString[j] <> ToString[k], i] -> chr[[i, j, k]]}],
{i, 1, n}, {j, 1, n}, {k, 1, j}];
TableForm[DeleteCases[Flatten[crs], Null]]

"gemischter Riemann Tensor"
rmn=ParallelTable[smp[
D[chr[[i, j, l]], x[[k]]] - D[chr[[i, j, k]], x[[l]]] +
Sum[chr[[s, j, l]] chr[[i, k, s]] -
chr[[s, j, k]] chr[[i, l, s]],
{s, 1, n}]], {i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, n}];
rie=ParallelTable[If[UnsameQ[rmn[[i, j, k, l]], 0],
{Subsuperscript["R", ToString[j] <> ToString[k] <> ToString[l], i] -> rmn[[i, j, k, l]]}],
{i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, k - 1}];
TableForm[DeleteCases[Flatten[rie], Null]]
(* kovarianter Riemann Tensor *)
rcv=ParallelTable[Sum[mt[[i, j]] rmn[[j, k, l, m]], {j, 1, n}],
{i, 1, n}, {k, 1, n}, {l, 1, n}, {m, 1, n}];
(* kontravarianter Riemann Tensor *)
rcn=ParallelTable[Sum[it[[m, i]] it[[h, j]] it[[o, k]] it[[p, l]] rcv[[i, j, k, l]],
{i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, n}],
{m, 1, n}, {h, 1, n}, {o, 1, n}, {p, 1, n}];

"Ricci Tensor"
rcc=ParallelTable[smp[
Sum[rmn[[i, j, i, l]], {i, 1, n}]], {j, 1, n}, {l, 1, n}];
Subscript["Ř", μσ] -> MatrixForm[rcc]
ric=ParallelTable[smp[Sum[
it[[i, k]] it[[j, l]] rcc[[k, l]], {k, 1, n}, {l, 1, n}]],
{i, 1, n}, {j, 1, n}];
"Ř"^μσ -> MatrixForm[ric]
rck=ParallelTable[smp[Sum[
it[[i, k]] rcc[[k, j]], {k, 1, n}]],
{i, 1, n}, {j, 1, n}];
Subsuperscript["Ř", "μ", "σ"] -> MatrixForm[rck]

"Ricci Skalar"
Ř=smp[Sum[it[[i, j]] rcc[[i, j]], {i, 1, n}, {j, 1, n}]]; "Ř"->Ř

"Kretschmann Skalar"
krn= smp[Sum[rcv[[i, j, k, l]] rcn[[i, j, k, l]],
{i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, n}]];
"K"->krn

"Bewegungsgleichungen"
geo=ParallelTable[smp[-Sum[
chr[[i, j, k]] x[[j]]' x[[k]]'+q f[[i, k]] x[[j]]' mt[[j, k]],
{j, 1, n}, {k, 1, n}]], {i, 1, n}];

equ=ParallelTable[{x[[i]]''==smp[rplc[geo[[i]]]]}, {i, 1, n}];

geodesic1=equ[[1]][[1]]
geodesic2=equ[[2]][[1]]
geodesic3=equ[[3]][[1]]
geodesic4=equ[[4]][[1]]

"totale Zeitdilatation"
H=Sum[mt[[μ, ν]] x[[μ]]' x[[ν]]', {μ, 1, n}, {ν, 1, n}];                       
ṫ=Quiet[rplc[smp[Normal[Solve[
-μ==(H/.t'->ť), ť]]]]];
Derivative[1][s]^2 == "ds²/dτ² == -μ" == smp[rplc[H]]                     
Derivative[1][t]->ṫ[[1, 1, 2]]||ṫ[[2, 1, 2]]||rplc[Sqrt[it[[1, 1]]]]/Sqrt[1-μ^2 v^2]

"kovarianter Viererimpuls"
p[μ_]:=-(Sum[mt[[μ, ν]]*x[[ν]]', {ν, 1, n}]+q A[[μ]]);
pt->rplc[smp[p[1]]]
pr->rplc[smp[p[2]]]
pθ->rplc[smp[p[3]]]
pφ->rplc[smp[p[4]]]

"linearer 3er-Impuls"
ṕ[i_]:=rplc[smp[Sqrt[-Sum[mt[[i, k]] x[[i]]' x[[k]]',
{k, 1, n}]]]]-Sqrt[1-μ^2 v^2] q A[[i]];
ṕr->ṕ[2]
ṕθ->ṕ[3]
ṕφ->ṕ[4]

"lokale Geschwindigkeit"
V[x_]:=smp[Normal[Solve[vx Sqrt[-mt[[x, x]]]/Sqrt[1-μ^2 v^2]-(1-μ^2 v^2) q A[[x]]==
p[x], vx]][[1, 1]]];
rplc[V[2]]/.vx->vr
rplc[V[3]]/.vx->vθ
rplc[V[4]]/.vx->vφ

"Frame Dragging Winkelgeschwindigkeit"
"ω"->-smp[mt[[1,4]]/mt[[4,4]]]

"Frame Dragging lokale Geschwindigkeit"
"ν"->smp[(Sqrt[mt[[1,4]]it[[1,4]]])]

"Horizonte & Ergosphären"
rh=smp[Solve[1/grr==0,r]]/.r->rH
re=smp[Solve[gtt==0,r]]/.r->rE

Quit[]

(* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
(*  KOORDINATEN TRANSFORMATOR|geodesics.yukterez.net|Kerr Schild (Finkelstein) -> Boyer Lindquist *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

d1={dT,dr,dθ,dΦ};
d2={dt,dr,dθ,dφ};
n=4;

"Kerr, Original Finkelstein Style"

g11=gtt=1-((2 M r)/((r^2)+(a^2) ((Cos[θ])^2)));
g22=grr=-1-((2 M r)/((r^2)+(a^2) ((Cos[θ])^2)));
g33=gθθ=-(r^2)-(a^2) ((Cos[θ])^2);
g44=gφφ=-(((a^2)+(r^2)) ((Sin[θ])^2))-((2 (a^2) M r ((Sin[θ])^4))/((r^2)+(a^2) ((Cos[θ])^2)));
g12=gtr=-((2 M r)/((r^2)+(a^2) ((Cos[θ])^2)));
g13=gtθ=0;
g14=gtφ=-((2 a M r ((Sin[θ])^2))/((r^2)+(a^2) ((Cos[θ])^2)));
g23=grθ=0;
g24=grφ=a (-1-((2 M r)/((r^2)+(a^2) ((Cos[θ])^2)))) ((Sin[θ])^2);
g34=gθφ=0;

m1={
{g11,g12,g13,g14},
{g12,g22,g23,g24},
{g13,g23,g33,g34},
{g14,g24,g34,g44}};

M1=MatrixForm[m1];
Subscript["g",μσ]->M1

A1={0,0,0,0};
Subscript["A",μ]->A1

"Kerr, Boyer Lindquist"

dT=dt+2M r dr/((r^2)-2M r+(a^2));
dΦ=dφ-a  dr/((r^2)-2M r+(a^2));

m2=FullSimplify[Table[Sum[D[d1[[k]],d2[[i]]] D[d1[[s]],d2[[j]]] m1[[k,s]],{k,1,n},{s,1,n}],{i,1,n},{j,1,n}],Reals];
M2=MatrixForm[m2];
Subscript["g",μσ]->M2

A2=FullSimplify[Table[Sum[D[d1[[k]],d2[[i]]] A1[[k]],{k,1,n}],{i,1,n}],Reals];
Subscript["A",μ]->A2