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

ClearAll["Local`*"]; LaunchKernels[4];
smp[y_]:=FullSimplify[y,Reals];
list[y_]:=y[[1]]==y[[2]];
rplc[y_]:=(((((((y/.t->t[τ])/.X->X[τ])/.Y->Y[τ])/.Z->Z[τ])/.Derivative[1][t[τ]]->
t'[τ])/.Derivative[1][r[τ]]->r'[τ])/.Derivative[1][Y[τ]]->Y'[τ])/.Derivative[1][Z[τ]]->Z'[τ]
rplc[y_]:=y

(* kovariante metrische Komponenten *)
g11=gtt=1/(((1+((E^(-((t (Sqrt[Λ]))/(Sqrt[3])))) M1)/(√((((X-X1))^2)+(((Y-Y1))^2)+(((Z-Z1))^2)))+((E^(-((t (Sqrt[Λ]))/(Sqrt[3])))) M2)/(√((((X-X2))^2)+(((Y-Y2))^2)+(((Z-Z2))^2)))))^2);
g22=gXX=-(E^((2 t (Sqrt[Λ]))/(Sqrt[3]))) (((1+((E^(-((t (Sqrt[Λ]))/(Sqrt[3])))) M1)/(√((((X-X1))^2)+(((Y-Y1))^2)+(((Z-Z1))^2)))+((E^(-((t (Sqrt[Λ]))/(Sqrt[3])))) M2)/(√((((X-X2))^2)+(((Y-Y2))^2)+(((Z-Z2))^2)))))^2);
g33=gYY=-(E^((2 t (Sqrt[Λ]))/(Sqrt[3]))) (((1+((E^(-((t (Sqrt[Λ]))/(Sqrt[3])))) M1)/(√((((X-X1))^2)+(((Y-Y1))^2)+(((Z-Z1))^2)))+((E^(-((t (Sqrt[Λ]))/(Sqrt[3])))) M2)/(√((((X-X2))^2)+(((Y-Y2))^2)+(((Z-Z2))^2)))))^2);
g44=gZZ=-(E^((2 t (Sqrt[Λ]))/(Sqrt[3]))) (((1+((E^(-((t (Sqrt[Λ]))/(Sqrt[3])))) M1)/(√((((X-X1))^2)+(((Y-Y1))^2)+(((Z-Z1))^2)))+((E^(-((t (Sqrt[Λ]))/(Sqrt[3])))) M2)/(√((((X-X2))^2)+(((Y-Y2))^2)+(((Z-Z2))^2)))))^2);
g12=gtX=0;
g13=gtY=0;
g14=gtZ=0;
g23=gXY=0;
g24=grZ=0;
g34=gYZ=0;

(* Dimensionen *)
x={t, X, Y, Z}; n=4; Λ=3(H^2);

"Metrischer Tensor"
mt=smp[{
{g11, g12, g13, g14},
{g12, g22, g23, g24},
{g13, g23, g33, g34},
{g14, g24, g34, g44}
}];
mtc=Subscript["g", μσ] -> MatrixForm[mt]
it=smp[Inverse[mt]];
mtk="g"^μσ -> MatrixForm[it]
mx=ParallelTable[smp[Sum[
it[[i, k]] mt[[k, j]], {k, 1, n}]],
{i, 1, n}, {j, 1, n}];
hp=Ap E^(-I ω (t-Z/c));hx=Ax E^(-I ω (t-Z/c));
mtc
mtk
"|g|"->Det[mt]

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

"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]]
geodesic3=equ[[3]][[1]]
geodesic4=equ[[4]][[1]]

"totale Zeitdilatation"
zd=Sum[mt[[μ, ν]] x[[μ]]' x[[ν]]', {μ, 1, n}, {ν, 1, n}];                       
ṫ=Quiet[rplc[smp[Normal[Solve[
-μ==(zd/.t'->ť), ť]]]]];
((Derivative[1][s])^2) =="ds²/dτ² == -μ" == smp[rplc[zd]]                     
Derivative[1][t]->ṫ[[1, 1, 2]]||ṫ[[2, 1, 2]]||rplc[Sqrt[(c^2)it[[1, 1]]]]/Sqrt[1-(μ^2) (((v/c))^2)]

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

"lokale Geschwindigkeit"
V[x_]:=smp[Normal[Solve[vx Sqrt[-mt[[x, x]]]/Sqrt[1-(μ^2) (((v/c))^2)]-(1-(μ^2) (((v/c))^2)) q A[[x]]==
p[x], vx]][[1, 1]]];                                                   
rplc[V[2]]/.vx->vX
rplc[V[3]]/.vx->vY
rplc[V[4]]/.vx->vZ