(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||| Mathematica Syntax | GEODESIC SOLVER | geodesics.yukterez.net | Version 21.01.2020 |||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

ClearAll["Local`*"]; smp[y_]:=Simplify[y, Reals]; list[y_]:=y[[1]]==y[[2]];
rplc[y_]:=(((((((y/.t->t[τ])/.r->r[τ])/.θ->θ[τ])/.φ->φ[τ])/.Derivative[1][t[τ]]->
t'[τ])/.Derivative[1][r[τ]]->r'[τ])/.Derivative[1][θ[τ]]->
θ'[τ])/.Derivative[1][φ[τ]]->φ'[τ]

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


                           (* Dimensionen, elektrische Ladung, Spin, Vakuumenergie, Masse *)
x={t, r, θ, φ}; n=4; Ω=℧; ℧=℧; a=a; Λ=0; M=1;

                                                                         "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]

"Maxwell Tensor"
A={((r ℧)/((r^2)+(a^2) ((Cos[θ])^2))),-((r ℧ (Sqrt[2 r-(℧^2)]))/((Sqrt[(a^2)+(r^2)]) ((a^2)+(-2+r) r+(℧^2)))),0,-((a r ℧ ((Sin[θ])^2))/((r^2)+(a^2) ((Cos[θ])^2)))};
F=smp[Table[((D[A[[j]], x[[k]]]-D[A[[k]], x[[j]]])), {j, 1, n}, {k, 1, n}]];
Subscript["F", μσ] -> MatrixForm[F]
f=smp[Table[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]

"Christoffelsymbole"
chr=smp[Table[(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=Table[If[UnsameQ[chr[[i, j, k]], 0],
{ToString[Γ[i, j, k]] "->", chr[[i, j, k]]}], {i, 1, n}, {j, 1, n}, {k, 1, j}];
TableForm[Partition[DeleteCases[Flatten[crs], Null], 2]]

"gemischter Riemann Tensor"
rmn=smp[Table[
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=Table[If[UnsameQ[rmn[[i, j, k, l]], 0],
{ToString[R[i, j, k, l]] "->", rmn[[i, j, k, l]]}],
{i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, k - 1}];
TableForm[Partition[DeleteCases[Flatten[rie], Null], 2]]
                                                            (* kovarianter Riemann Tensor *)
rcv=Table[Sum[mt[[i, j]] rmn[[j, k, l, m]], {j, 1, 4}],
{i, 1, n}, {k, 1, n}, {l, 1, n}, {m, 1, n}];
                                                        (* kontravarianter Riemann Tensor *)
rcn=Table[Sum[it[[m, i]] it[[h, j]] it[[o, k]] it[[p, l]] rcv[[i, j, k, l]],
{i, 1, 4}, {j, 1, n}, {k, 1, n}, {l, 1, n}],
{m, 1, 4}, {h, 1, n}, {o, 1, n}, {p, 1, n}];

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

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

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

"Einstein Tensor"
est=smp[ric-Ř mt/2];
Subscript["G", μσ] -> MatrixForm[est]
ein=smp[Table[Sum[
mt[[i, k]] mt[[j, l]] est[[k, l]], {k, 1, n}, {l, 1, n}],
{i, 1, n}, {j, 1, n}]];
"G"^μσ -> MatrixForm[smp[ein]]

"Stress Energie Tensor"
set=smp[est+Λ mt]/8/π;                                                   
Subscript["T", μσ] -> MatrixForm[set]
sei=smp[Table[Sum[
mt[[i, k]] mt[[j, l]] set[[k, l]], {k, 1, n}, {l, 1, n}],
{i, 1, n}, {j, 1, n}]];
"T"^μσ -> MatrixForm[smp[sei]]

"Bewegungsgleichungen"
geo=smp[Table[-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=Table[{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][t][τ]->ṫ[[1, 1, 2]] || ṫ[[2, 1, 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]]]

"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φ[τ]