(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||| Koordinatentransformation von Schwarzschild Droste in isoptrope Koordinaten |||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

d1={dt,dr,dθ,dφ}; (* Droste *)
d2={dt,dř,dθ,dφ}; (* Isotrop *)
n=4; rs=2M;

dr=((dř (Sqrt[r]) (Sqrt[r-rs]))/ř) ; (* Differential *)
я=ř(((1+M/2/ř))^2); (* r in Terms von ř *)

(* Droste Koeffizienten *)
g11=gtt=1-rs/r;
g22=grr=-1/(1-rs/r);
g33=gθθ=-(r^2);
g44=gφφ=-(r^2) ((Sin[θ])^2);
g12=g13=g14=g23=g24=g34=0;

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

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}]];
M2=Simplify[MatrixForm[m2/.r->я]];

M1->M2 (* links: Droste, rechts: isotrope Koordinaten *)

Quit[]

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

smp[y_]:=Simplify[FullSimplify[y,Reals],r>0&&M>0&&0<θ<π/2]; list[y_]:=y[[1]]==y[[2]];
rplc[y_]:=(((((((y/.t->t[τ])/.ř->ř[τ])/.θ->θ[τ])/.φ->φ[τ])/.Derivative[1][t[τ]]->
t'[τ])/.Derivative[1][ř[τ]]->ř'[τ])/.Derivative[1][θ[τ]]->
θ'[τ])/.Derivative[1][φ[τ]]->φ'[τ]

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

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

"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=smp[Inverse[mt]];
"g"^μσ->MatrixForm[it]
mx=ParallelTable[smp[Sum[it[[i,k]] mt[[k,j]],{k,1,n}]],{i,1,n},{j,1,n}];
Subsuperscript["g","μ","σ"]->MatrixForm[mx]

"Maxwell Tensor"
A={0,0,0,0};
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= 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[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]]

"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]]]
př[τ]->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->vř[τ]
rplc[V[3]]/.vx->vθ[τ]
rplc[V[4]]/.vx->vφ[τ]

r->(((1+(M/(2 ř))))^2) ř