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

ClearAll["Local`*"];
smp[y_]:=Simplify[y, Reals];
list[y_]:=y[[1]]==y[[2]];
rplc[y_]:=y

(* kovariante metrische Komponenten *)
g11=guu=(32/r) Exp[-r/2]; 
g22=gww=-guu; g12=0;

r=2 (1+ProductLog[(((u^2)-(w^2))/E)]);

(* Dimensionen *)
x={u, w}; n=2;

"Metrischer Tensor"
mt=smp[{
{g11, g12},
{g12, g22}}];
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]
md=Det[mt]; "|g|" -> smp[md]

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

"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 (* -Λ mt *) )/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]
Ť = smp[Sum[it[[i, j]] set[[i, j]], {i, 1, n}, {j, 1, n}]]; "T" -> Ť

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

"totale Zeitdilatation"
H=Sum[mt[[μ, ν]] x[[μ]]' x[[ν]]', {μ, 1, n}, {ν, 1, n}];                       
ṫ=Quiet[rplc[smp[Normal[Solve[
-μ==(H/.u'->ť), ť]]]]];
Derivative[1][s]^2 == "ds²/dτ² == -μ" == smp[rplc[H]]                     
Derivative[1][u]->ṫ[[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[[μ]]);
pu->rplc[smp[p[1]]]
pw->rplc[smp[p[2]]]

"linearer Impuls"
ṕ[i_]:=rplc[smp[Sqrt[-Sum[mt[[i, k]] x[[i]]' x[[k]]', {k, 1, n}]]]];
ṕu->ṕ[1]
ṕw->ṕ[2]

"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[1]]/.vx->vu
rplc[V[2]]/.vx->vw

"Transformation"
R=r; Clear[r];
Eq[u_,w_,r_,t_]:={u==(Sqrt[(r/2)-1])Exp[(r/4)]Cosh[(t/4)]&&w==(Sqrt[(r/2)-1])Exp[(r/4)]Sinh[(t/4)]};
uw[r_,t_]:=N[FindInstance[Eq[u,w,r,t],{u,w}, Reals,WorkingPrecision->24]];
rt[u_,w_]:=N[FindInstance[Eq[u,w,r,t],{r,t}, Reals,WorkingPrecision->24]];

u=3; w=1;
N[rt[u,w]]

Clear[u, w]

t=1.3862943611199;r=4.0802251644714;
N[uw[r,t]]

Clear[r,t]

v=Simplify[Solve[0==guu (du^2)+gww (dw^2), du],r>0]
û=v[[2,1,2]]/dw;

U=Integrate[1/û, {w,w1,w2}]

p1=Show[ParallelTable[Plot[
Evaluate[{+1,-1}U],{w2,-10,10},
Frame->True, PlotStyle->{{Blue,Thick},{Red,Thick}},ImageSize->480, 
AspectRatio->1,PlotRange->{{-5,5},{-5,5}}],
{w1,-10,10,2}]]

p2=Show[ParallelTable[ContourPlot[
(r/2-1)Exp[r/2]==(u^2)-(w^2),{u,-5,5},{w,-5,5},
ImageSize->476, ContourStyle->{Black,Thickness[0.002]},
AspectRatio->1,PlotRange->{{-5,5},{-5,5}},
MaxRecursion->10, WorkingPrecision->32],
{r,{5,4,3,2,0}}]]

p3=Show[Table[Graphics[{Opacity[1],InfiniteLine[
{0,0},{uw[4,t][[1,1,2]],uw[4,t][[1,2,2]]}]},
PlotRange->{{-5,5},{-5,5}},ImageSize->476,
Frame->True],{t,-4,4}]]

p4=Show[Table[Graphics[{Opacity[1],InfiniteLine[
{0,0},{uw[4,t][[1,2,2]],uw[4,t][[1,1,2]]}]},
PlotRange->{{-5,5},{-5,5}},ImageSize->476,
Frame->True],{t,-4,4}]]

p5=Plot[w/.NSolve[R==0,w][[2]],{u,-5,5}, WorkingPrecision->24,
 PlotStyle->{Opacity[1],Thick,Darker[Red]},
Filling->{Top},FillingStyle->Directive[Opacity[1],Red],
PlotRange->{{-5,5},{-5,5}},ImageSize->476,
Frame->True,AspectRatio->1]

p6=Plot[-w/.NSolve[R==0,w][[2]],{u,-5,5}, WorkingPrecision->24,
 PlotStyle->{Opacity[1],Thick,Darker[Darker[Green]]},
Filling->{Bottom},FillingStyle->Directive[Opacity[1],Darker[Green]],
PlotRange->{{-5,5},{-5,5}},ImageSize->476,
Frame->True,AspectRatio->1]

Show[p1,p2,p3,p4,p5,p6]