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

ClearAll["Global`*"]
                                               (* Input: kovariante metrische Komponenten *)
gtt=1-(2r-℧^2)/Σ;
grr=-Σ/Δ;
gθθ=-Σ;
gφφ=-Χ/Σ Sin[θ]^2;
gtφ=a (2r-℧^2) Sin[θ]^2/Σ;
                                                                           (* Abkürzungen *)
Σ=r^2+a^2 Cos[θ]^2;
Δ=r^2-2 r+a^2+℧^2;
Χ=(r^2+a^2)^2-a^2 Sin[θ]^2 Δ;
                    (* Koordinaten, Dimensionen, magnetisches Monopol, elektrische Ladung *)
x={t, r, θ, φ}; n=4; P=0; Ω=℧;
                                                                         "Metrischer Tensor"
metrik={{gtt, 0, 0, gtφ}, {0, grr, 0, 0}, {0, 0, gθθ, 0}, {gtφ, 0, 0, gφφ}};
Subscript["g", μσ] -> MatrixForm[metrik]

inversemetrik=Simplify[Inverse[metrik]];
"g"^μσ -> MatrixForm[inversemetrik]
                                                                            "Maxwell Tensor"
A={(Ω r)/Σ+P/Σ Cos[θ] a, 0, 0, -(Ω r/Σ) Sin[θ]^2 a-P/Σ Cos[θ](r^2+a^2)};
F=Simplify[Table[((D[A[[j]], x[[k]]]-D[A[[k]], x[[j]]])), {j, 1, 4}, {k, 1, 4}], Reals];
Subscript["F", μσ] -> MatrixForm[F]

f=FullSimplify[Table[Sum[
inversemetrik[[i, k]] inversemetrik[[j, l]] F[[k, l]], 
{k, 1, 4}, {l, 1, 4}], {i, 1, 4}, {j, 1, 4}], Reals];
"F"^μσ -> MatrixForm[f]
                                                                           "Einstein Tensor"
G=FullSimplify[Table[Sum[
-2(inversemetrik[[k,l]] F[[i, k]] F[[j, l]] - metrik[[i, j]] F[[k, l]] f[[k, l]]), 
{k, 1, 4}, {l, 1, 4}], {i, 1, 4}, {j, 1, 4}], Reals];
Subscript["G", μσ] -> MatrixForm[G]

H=FullSimplify[Table[Sum[
inversemetrik[[i, k]] inversemetrik[[j, l]] G[[k, l]], {k, 1, 4}, {l, 1, 4}], 
{i, 1, 4}, {j, 1, 4}], Reals];
"G"^μσ -> MatrixForm[H]
                                                                       "Christoffel Symbole"
christoffel:=Simplify[Table[(1/2)Sum[(inversemetrik[[i, s]])
(D[metrik[[s, j]], x[[k]]]+D[metrik[[s, k]], x[[j]]] -D[metrik[[j, k]], x[[s]]]), {s, 1, n}], 
{i, 1, n}, {j, 1, n}, {k, 1, n}]];
christoffelsymbole=Table[If[UnsameQ[christoffel[[i, j, k]], 0], 
{ToString[Γ[i, j, k]], christoffel[[i, j, k]]}], {i, 1, n}, {j, 1, n}, {k, 1, j}];

rplc[y_]:=(((((((y/.t->t[τ])/.r->r[τ])/.θ->θ[τ])/.φ->φ[τ])/.Derivative[1][t[τ]]-> 
t'[τ])/.Derivative[1][r[τ]]->r'[τ])/.Derivative[1][θ[τ]]->θ'[τ])/.Derivative[1][φ[τ]]->φ'[τ];
rple[y_]:=(((((((y/.t->t[τ])/.r->r[τ])/.θ->θ[τ])/.φ->φ[τ])/.Derivative[1][t[τ]]->  
t'[τ])/. r\.b4->r'[τ])/. θ\.b4->θ'[τ])/. φ\.b4->φ'[τ];
list[y_]:=y[[1]]==y[[2]];

list[christoffelsymbole[[1]][[2]][[1]]]
list[christoffelsymbole[[1]][[3]][[1]]]
list[christoffelsymbole[[1]][[4]][[2]]]
list[christoffelsymbole[[1]][[4]][[3]]]
list[christoffelsymbole[[2]][[1]][[1]]]
list[christoffelsymbole[[2]][[2]][[2]]]
list[christoffelsymbole[[2]][[3]][[2]]]
list[christoffelsymbole[[2]][[3]][[3]]]
list[christoffelsymbole[[2]][[4]][[1]]]
list[christoffelsymbole[[2]][[4]][[4]]]
list[christoffelsymbole[[3]][[1]][[1]]]
list[christoffelsymbole[[3]][[2]][[2]]]
list[christoffelsymbole[[3]][[3]][[2]]]
list[christoffelsymbole[[3]][[3]][[3]]]
list[christoffelsymbole[[3]][[4]][[1]]]
list[christoffelsymbole[[3]][[4]][[4]]]
list[christoffelsymbole[[4]][[2]][[1]]]
list[christoffelsymbole[[4]][[3]][[1]]]
list[christoffelsymbole[[4]][[4]][[2]]]
list[christoffelsymbole[[4]][[4]][[3]]]
                                                                      "Bewegungsgleichungen"
geodäsie=Simplify[Table[-Sum[
christoffel[[i, j, k]] x[[j]]' x[[k]]'+q f[[i, k]] x[[j]]' metrik[[j, k]], 
{j, 1, n}, {k, 1, n}], {i, 1, n}]];

bewegungsgleichung=Table[{x[[i]]''[τ]==FullSimplify[rplc[geodäsie[[i]]], Reals]}, {i, 1, n}];

bewegungsgleichung[[1]][[1]]
bewegungsgleichung[[2]][[1]]
bewegungsgleichung[[3]][[1]]
bewegungsgleichung[[4]][[1]]

ClearAll[Σ, Δ, Χ];
                                                                            "Zeitdilatation"
t'[τ]==rple[Simplify[Normal[Solve[
-μ==gtt Т^2+grr r\.b4^2+gθθ θ\.b4^2+gφφ φ\.b4^2 + 2 gtφ Т φ\.b4, Т, Reals]]][[2]][[1]][[2]]]