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

smp[y_]:=Simplify[y,Reals&&r>0&&M>0];
list[y_]:=y[[1]]==y[[2]];
rplc[y_]:=y;

(* kovariante metrische Komponenten *)
g11=guu=1-((2 r-(℧^2))/((r^2)+(a^2) ((Cos[θ])^2)));
g22=grr=0;
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)+(℧^2)) ((Sin[θ])^4))/((r^2)+(a^2) ((Cos[θ])^2)));
g12=gur=-1;
g13=guθ=0;
g14=guφ=((a (2 r-(℧^2)) ((Sin[θ])^2))/((r^2)+(a^2) ((Cos[θ])^2)));
g23=grθ=0;
g24=grφ=a ((Sin[θ])^2);
g34=gθφ=0;

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

"Metrischer Tensor"
mt=smp[{
{g11, g12, g13, g14},
{g12, g22, g23, g24},
{g13, g23, g33, g34},
{g14, g24, g34, g44}
}];
Subscript["g", μσ] -> MatrixForm[mt]
it=Simplify[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={((r ℧)/((r^2)+(a^2) ((Cos[θ])^2))),-((r ℧)/((a^2)-2 r+(r^2)+(℧^2))),0,-((a r ℧ ((Sin[θ])^2))/((r^2)+(a^2) ((Cos[θ])^2)))}

"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, 4}, {j, 1, n}, {k, 1, n}, {l, 1, n}]];
"K"->krn

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

"Stress Energie Impuls Tensor"
set=smp[est/8/π];
Subscript["T",μσ]->MatrixForm[set]
sei=ParallelTable[smp[Sum[it[[i,k]] it[[j,l]] set[[k,l]],{k,1,n},{l,1,n}]],{i,1,n},{j,1,n}];
"T"^μσ->MatrixForm[sei]
sem=ParallelTable[smp[Sum[it[[i,k]] set[[k,j]],{k,1,n}]],{i,1,n},{j,1,n}];
Subsuperscript["T","μ","σ"]->MatrixForm[sem]

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

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

Quit[]

(* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
(* ||| KOORDINATEN TRANSFORMATOR | geodesics.yukterez.net | Boyer Lindquist -> Nullkoordinaten |||*)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

"Kerr Newman, Boyer Lindquist"

M=1;

g11=gtt=1-(2 M r-℧^2)/ (r^2+a^2 Cos[θ]^2);
g22=grr=-(r^2+a^2 Cos[θ]^2)/(a^2-2M r+r^2+℧^2);
g33=gθθ=-r^2-a^2 Cos[θ]^2;
g44=gφφ=-((a^2+r^2)^2 Sin[θ]^2-a^2 (a^2-2M r+r^2+℧^2) Sin[θ]^4)/(r^2+a^2 Cos[θ]^2);
g14=gtφ=-(a (℧^2-2 M r) Sin[θ]^2)/(r^2+a^2 Cos[θ]^2);
g12=g13=g23=g24=g34=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={r ℧/(r^2+a^2 Cos[θ]^2), 0, 0, -a r ℧ Sin[θ]^2/(r^2+a^2 Cos[θ]^2)};
Subscript["A", μ] -> A1

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

(* Transformationsregel *)

Σ=r^2+a^2 Cos[θ]^2;
Δ=(r^2+a^2)(1-Λ/3 r^2)-2 M r+℧^2;
Χ=(r^2+a^2)^2-a^2 Sin[θ]^2 Δ;
щ=(q ℧ r (a^2+r^2))/(Δ Σ);
χ=1+Λ/3 a^2;
ж=1+Λ/3 a^2 Cos[θ]^2;
σ=a^2+r^2;
Λ=0;

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

"Kerr Newman, Nullkoordinaten"

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&&r>0];
Subscript["A", μ] -> A2