Kerr Newman Metrik

[Yukterez]

Das ist die deutschsprachige Version.   English versions are available on en.yukterez.net and yukipedia.

Verwandte Beiträge: Kerr || Reissner Nordström || Schwarzschild || Geodäsie || Gravitationslinsen || Raytracing || Paraboloid  
 

  ←
Schatten eines Kerr-Newman SL mit a²+℧²=M², Betrachtungswinkel: edge on. Für andere Winkel siehe hier. Rohmaterial: ESO/Brunier.   
 

  ←
Zoom auf den Schatten des obigen Kerr-Newman SL mit a²+℧²=M² und eingeblendeten Oberflächen (letztere sind in natura nicht sichtbar)    
 

  ←
Retrograder Orbit eines mit q=1 geladenen Partikels um ein SL mit Spin & Ladung a=√¾ & ℧=⅓. v0 & i0: lokale Startgeschwindigkeit & Inklination    

  ←
Prograder Orbit eines neutralen Testpartikels um ein mit a=0.9 rotierendes und ℧=0.4 elektrisch geladenes Kerr Newman Schwarzloch    

  ←
Selbe Startbedingungen wie oben (R0=x0=7GM/c², v0=0.4c, i0=39.8056°=atan(5/6)rad), aber mit negativ geladenem Testpartikel (q=-¼)      

  ←
Orbit um ein SL mit a=√½ und ℧=√½. Partikel: q=⅓, Lz=0, Startbedingungen: v0=0.57: vr=0, vθ=√(782759/2409750), vφ=-√(6751/96390000)     

  ←
Orbit um ein SL mit a=0.8 und ℧=½. Partikel: q=-½, Lz=0 (man beachte dass die lokale φ-Geschwindigkeit trotz Lz=0 wegen q≠0 hier nicht 0 ist)     

  ←
Orbit eines mit q=½ geladenen Partikels um eine nackte Singularität mit a=¾ und ℧=⅔; farbige Flächen: äußere und innere Ergosphäre      

  ←
Freier Fall eines neutralen Testpartikels auf eine mit mit a=1.5 rotierende und und ℧=0.4 elektrisch geladene nackte Singularität      

  ←
Flucht von neutralen und geladenen Partikeln aus dem Inneren der Ergosphäre, Startposition knapp oberhalb des äußeren Horizonts      

  ←
Orbit eines negativ geladenen Partikels im Inneren eines positiv geladenen Reissner Nordström SL (vergleiche Dokuchaev, Fig. 1)      

[Yukterez]

Linienelement in Boyer-Lindquist-Koordinaten, metrische Signatur (+,-,-,-):



wobei



mit dem Spinparameter â=Jc/G/M bzw. dem dimensionslosen Spinparameter a=â/M und dem spezifischen elektrischen Ladungsparameter Ω=⚼·√(K/G) bzw. dem dimensionlosen Ladungsparameter ℧=Ω/M. Der Artikel verwendet die dimensionslosen natürlichen Einheiten G=M=c=K=1 mit Längen in GM/c² und Zeiten in GM/c³. Der Zusammenhang zwischen dem hier gleich 1 gesetzten Massenäquivalent M und der irreduziblen Masse Mirr ist



Für massebehaftete Testpartikel gilt μ=-1, für Photonen μ=0. Die spezifische Ladung des Testpartikels ist q. Transformationsregel für ko-und kontravariante Indizes (hochgestellte Buchstaben sind hierbei keine Potenzen sondern Indizes):



Ko- und kontravariante Metrik:

Elektromagnetisches Potential:



Kovarianter elektromagnetischer Tensor:



Kontravarianter Maxwelltensor:



Eigenzeitableitungen der Koordinaten:



Damit lauten die Bewegungsgleichungen:









Kanonischer Viererimpuls, lokale Dreiergeschwindigkeit und Koordinatenschnelligkeit:

Insgesamte Zeitdilatation eines neutralen Testpartikels, gewonnen aus dem Linienelement:



Zusätzliche Zeitdilatation eines geladenen Testpartikels:



Insgesamte Zeitdilatation eines geladenen Testpartikels:



Zusammenhang der ersten Eigenzeitableitungen mit den kovarianten Impulskomponenten:

   



Zusammenhang der ersten Eigenzeitableitungen mit den lokalen Dreiergeschwindigkeitskomponenten:

   



Für die lokale Dreiergeschwindigkeit relativ zu einem ZAMO vor Ort wird E nach v aufgelöst:



Radiale Fluchtgeschwindigkeit für ein neutrales Teilchen:



Für die Fluchtgeschwindigkeit eines drehimpulsfreien geladenen Teilchens wird E=1 gesetzt und nach v aufgelöst:



Der dafür benötigte lokale vertikale Abschusswinkel ist π/2 wenn q=0, und wenn q≠0:

1. Erhaltungsgröße: Gesamtenergie E=-pt

2. Erhaltungsgröße: axialer Drehimpuls Lz=+pφ

3. Erhaltungsgröße: Carter-Konstante



mit der koaxialen Drehimpulskomponente (die selber keine Erhaltungsgröße ist)



Das radiale effektive Potential das an seinen Nullstellen die Umkehrpunkte vorgibt ist



mit dem Parameter



Gravitative Zeitdilatation eines ZAMO, unendlich am Horizont:



Zeitdilatation eines neutralen stationären Partikels, unendlich bei der Ergosphäre:



Axialer und koaxialer Gyrationsradius:



Axialer und koaxialer Umfang:



Frame-Dragging Winkelgeschwindigkeit im System den Koordinatenbuchhalters:



Lokale Frame-Dragging Geschwindigkeit relativ zu den Fixsternen (bei der Ergosphäre c):



Radialkoordinaten der Horizonte und Ergosphären:



Kartesische Projektion:

[Yukterez]

Input: Linienelement und Maxwell-Tensor, Output: Bewegungsgleichungen

Code: Alles auswählen

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

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 Δ;
щ=(q ℧ r (a^2+r^2))/(Δ Σ);

              (* Koordinaten, Dimensionen, magnetisches Monopol, elektrische Ladung, Spin *)
x={t, r, θ, φ}; n=4; P=0; Ω=℧; ℧=℧; a=a;

                                                                         "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, n}, {k, 1, n}], Reals];
Subscript["F", μσ] -> MatrixForm[F]
f=Simplify[Table[Sum[
inversemetrik[[i, k]] inversemetrik[[j, l]] F[[k, l]],
{k, 1, n}, {l, 1, n}], {i, 1, n}, {j, 1, n}], Reals];
"F"^μσ -> MatrixForm[f]

                                                                        "Christoffelsymbole"
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}]];
Christoffel=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}];
TableForm[Partition[DeleteCases[Flatten[Christoffel], Null], 2]]
                                               
                                                                            "Riemann Tensor"
riemann=Simplify[Table[
D[christoffel[[i, j, l]], x[[k]]] - D[christoffel[[i, j, k]], x[[l]]] +
Sum[christoffel[[s, j, l]] christoffel[[i, k, s]] -
christoffel[[s, j, k]] christoffel[[i, l, s]],
{s, 1, n}], {i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, n}]];
Riemann=Table[If[UnsameQ[riemann[[i, j, k, l]], 0],
{ToString[R[i, j, k, l]], riemann[[i, j, k, l]]}],
{i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, k - 1}];
TableForm[Partition[DeleteCases[Flatten[Riemann], Null], 2]]
 
                                                                              "Ricci Tensor"
ricci=Simplify[Table[
Sum[riemann[[i, j, i, l]], {i, 1, n}], {j, 1, n}, {l, 1, n}]];
Subscript["Ř", μσ] -> MatrixForm[ricci]
Ricci=Simplify[Table[Sum[
inversemetrik[[i, k]] inversemetrik[[j, l]] ricci[[k, l]], {k, 1, n}, {l, 1, n}],
{i, 1, n}, {j, 1, n}], Reals];
"Ř"^μσ -> MatrixForm[Ricci]
                                                                              "Ricci Skalar"
Ř=Simplify[Sum[inversemetrik[[i, j]] ricci[[i, j]], {i, 1, n}, {j, 1, n}]]; "Ř"->Ř

"Einstein Tensor"
einstein=Simplify[Ricci-Ř metrik/2];
Subscript["G", μσ] -> MatrixForm[einstein]
Einstein=Simplify[Table[Sum[
metrik[[i, k]] metrik[[j, l]] einstein[[k, l]], {k, 1, n}, {l, 1, n}],
{i, 1, n}, {j, 1, n}], Reals];
"G"^μσ -> MatrixForm[Simplify[Einstein]]

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

                                                                      "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]]''[τ]==Simplify[rplc[geodäsie[[i]]], Reals]}, {i, 1, n}];

geodesic1=bewegungsgleichung[[1]][[1]]
geodesic2=bewegungsgleichung[[2]][[1]]
geodesic3=bewegungsgleichung[[3]][[1]]
geodesic4=bewegungsgleichung[[4]][[1]]

                                                                            "Zeitdilatation"
ṫ=rple[Simplify[Normal[Solve[
-μ==gtt ť^2+grr r\.b4^2+gθθ θ\.b4^2+gφφ φ\.b4^2 + 2 gtφ ť φ\.b4, ť, Reals]]][[2, 1, 2]]-щ];
Derivative[1][t][τ]==ṫ











Simulator-Code für Photonen, geladene und neutrale Teilchen

Code: Alles auswählen

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||| Mathematica | kerr.newman.yukterez.net | 06.08.2017 - 18.04.2019, Version 14 |||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]; ClearAll["Local`*"];
Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
 
mt1={"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}};
mt2={"EventLocator", "Event"-> (r[τ]-1001/1000 rA)};
mt3={"ImplicitRungeKutta", "DifferenceOrder"-> 20};
mt4={"EquationSimplification"-> "Residual"};
mt0=Automatic;
mta=mt0;
wp=MachinePrecision;
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 1) STARTBEDINGUNGEN EINGEBEN |||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
A=a;                                   (* pseudosphärisch [BL]: A=0, kartesisch [KS]: A=a *)
 
tmax=300;                                                                    (* Eigenzeit *)
Tmax=300;                                                              (* Koordinatenzeit *)
TMax=Min[Tmax, т[plunge-1*^-3]]; tMax=Min[tmax, plunge-1/1000];       (* Integrationsende *)
 
r0 = Sqrt[7^2-a^2];                                                        (* Startradius *)
r1 = r0+1;                                                   (* Endradius wenn v0=vr0=vr1 *)
θ0 = π/2;                                                                  (* Breitengrad *)
φ0 = 0;                                                                     (* Längengrad *)
a  = 9/10;                                                               (* Spinparameter *)
℧  = 2/5;                                       (* spezifische Ladung des schwarzen Lochs *)
q  = 0;                                             (* spezifische Ladung des Testkörpers *)
 
v0 = 2/5;                                                       (* Anfangsgeschwindigkeit *)
α0 = 0;                                                      (* vertikaler Abschusswinkel *)
i0 = ArcTan[5/6];                                               (* Bahninklinationswinkel *)
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 2) GESCHWINDIGKEITS-, ENERGIE UND DREHIMPULSKOMPONENTEN ||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
vr0=v0 Sin[α0];                                     (* radiale Geschwindigkeitskomponente *)
vθ0=v0 Cos[α0] Sin[i0];                      (* longitudinale  Geschwindigkeitskomponente *)
vφ0=v0 Cos[α0] Cos[i0];                        (* latitudinale Geschwindigkeitskomponente *)

vr[τ_]:=R'[τ]/Sqrt[Δi[τ]] Sqrt[Σi[τ] (1+μ v[τ]^2)];
vθ[τ_]:=Θ'[τ] Sqrt[Σi[τ] (1+μ v[τ]^2)];
vφ[τ_]:=-((Sin[Θ[τ]] Sqrt[1+μ v[τ]^2] (-a^5 ε Cos[Θ[τ]]^2+a q ℧ R[τ]^3-a ε R[τ]^4+
a q ℧ R[τ] (a^2-Δi[τ])+a^2 Δi[τ] (xJ[τ] ε Cot[Θ[τ]]^2+Φ'[τ] Cos[Θ[τ]]^2 Σi[τ])+
R[τ]^2 (-a^3 ε (1+Cos[Θ[τ]]^2)+Δi[τ] (xJ[τ] ε Csc[Θ[τ]]^2+Φ'[τ] Σi[τ]))))/((a^2 Cos[Θ[τ]]^2+
R[τ]^2) (a^2 Sin[Θ[τ]]^2-Δi[τ]) Sqrt[Χi[τ]/Σi[τ]]));
vΦ[τ_]:=Sqrt[v[τ]^2-vθ[τ]^2-vr[τ]^2];
 
x0[A_]:=Sqrt[r0^2+A^2] Sin[θ0] Cos[φ0];                             (* Anfangskoordinaten *)
y0[A_]:=Sqrt[r0^2+A^2] Sin[θ0] Sin[φ0];
z0[A_]:=r0 Cos[θ0];
 
ε0=Sqrt[δ Ξ/χ]/j[v0]+Lz ω0;
ε=ε0+((q r0 ℧)/(r0^2+a^2 Cos[θ0]^2));
εζ:=Sqrt[Δ Σ/Χ]/j[ν]+Lz ωζ+((q r[τ] ℧)/(r[τ]^2+a^2 Cos[θ[τ]]^2));
LZ=vφ0 Ы/j[v0];
Lz=LZ+((q a r0 ℧ Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2));
Lζ:=vφ0 я/j[ν]+((q a r[τ] ℧ Sin[θ[τ]]^2)/(r[τ]^2+a^2 Cos[θ[τ]]^2));
pθ0=vθ0 Sqrt[Ξ]/j[v0]; pθζ:=θ'[τ] Σ;
pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];
Qk=Limit[pθ0^2+(Lz^2 Csc[θ1]^2-a^2 (ε^2+μ)) Cos[θ1]^2, θ1->θ0];       (* Carter Konstante *)
Q=Limit[pθ0^2+(Lz^2 Csc[θ1]^2-a^2 (ε^2+μ)) Cos[θ1]^2, θ1->θ0];
Qζ:=pθζ^2+(Lz^2 Csc[θ[τ]]^2-a^2 (εζ^2+μ)) Cos[θ[τ]]^2;
k=Q+Lz^2+a^2 (ε^2+μ); kζ:=Qζ+Lz^2+a^2 (εζ^2+μ);
 
μ=If[Abs[v0]==1, 0, -1];                                     (* Baryon: μ=-1, Photon: μ=0 *)

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 3) FLUCHTGESCHWINDIGKEIT UND DAFÜR BENÖTIGTER VERTIKALER ABSCHUSSWINKEL ||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

vEsc=If[q==0, ж0, Abs[(\[Sqrt](r0^2 (r0^2 (δ Ξ-χ)+2 q r0 χ ℧-q^2 χ ℧^2)+
2 a^2 r0 (r0 δ Ξ-r0 χ+q χ ℧) Cos[θ0]^2+a^4 (δ Ξ-
χ) Cos[θ0]^4))/(Sqrt[χ] (r0 (r0-q ℧)+a^2 Cos[θ0]^2))]];

αEsc=If[q==0, π/2, ArcCos[-((a q r0 ℧ Sqrt[1-Abs[Sqrt[-a^4 (a^2+r0^2) (2 r0-℧^2) Cos[θ0]^4+
2 a^2 r0 Cos[θ0]^2 (-r0 (a^2+r0^2)^2+q (a^2+r0^2)^2 ℧+r0 (a^2-2 r0+r0^2+℧^2) (r0^2+
a^2 Cos[θ0]^2)+a^2 r0 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2-a^2 q ℧ (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)+
r0^2 (-r0^2 (a^2+r0^2) (2 r0-℧^2)+2 q r0 ℧ ((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)-
q^2 ℧^2 ((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))]/((r0 (r0-q ℧)+
a^2 Cos[θ0]^2) Sqrt[(a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+
℧^2) Sin[θ0]^2])]^2] Sin[θ0])/(Abs[Sqrt[-a^4 (a^2+r0^2) (2 r0-℧^2) Cos[θ0]^4+
2 a^2 r0 Cos[θ0]^2 (-r0 (a^2+r0^2)^2+q (a^2+r0^2)^2 ℧+r0 (a^2-2 r0+r0^2+℧^2) (r0^2+
a^2 Cos[θ0]^2)+a^2 r0 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2-a^2 q ℧ (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)+
r0^2 (-r0^2 (a^2+r0^2) (2 r0-℧^2)+2 q r0 ℧ ((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+
℧^2) Sin[θ0]^2)-q^2 ℧^2 ((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))]/((r0 (r0-q ℧)+
a^2 Cos[θ0]^2) Sqrt[(a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2])] (r0^2+
a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+
a^2 Cos[θ0]^2)]))]];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 4) HORIZONTE UND ERGOSPHÄREN RADIEN ||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
rE=1+Sqrt[1-a^2 Cos[θ]^2-℧^2];                                       (* äußere Ergosphäre *)
RE[A_, w1_, w2_]:=Xyz[xyZ[
{Sqrt[rE^2+A^2] Sin[θ] Cos[φ], Sqrt[rE^2+A^2] Sin[θ] Sin[φ], rE Cos[θ]}, w1], w2];
rG=1-Sqrt[1-a^2 Cos[θ]^2-℧^2];                                       (* innere Ergosphäre *)
RG[A_, w1_, w2_]:=Xyz[xyZ[
{Sqrt[rG^2+A^2] Sin[θ] Cos[φ], Sqrt[rG^2+A^2] Sin[θ] Sin[φ], rG Cos[θ]}, w1], w2];
rA=1+Sqrt[1-a^2-℧^2];                                                 (* äußerer Horizont *)
RA[A_, w1_, w2_]:=Xyz[xyZ[
{Sqrt[rA^2+A^2] Sin[θ] Cos[φ], Sqrt[rA^2+A^2] Sin[θ] Sin[φ], rA Cos[θ]}, w1], w2];
rI=1-Sqrt[1-a^2-℧^2];                                                 (* innerer Horizont *)
RI[A_, w1_, w2_]:=Xyz[xyZ[
{Sqrt[rI^2+A^2] Sin[θ] Cos[φ], Sqrt[rI^2+A^2] Sin[θ] Sin[φ], rI Cos[θ]}, w1], w2];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 5) HORIZONTE UND ERGOSPHÄREN PLOT ||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
horizons[A_, mesh_, w1_, w2_]:=Show[
ParametricPlot3D[RE[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
Mesh -> mesh, PlotPoints -> plp, PlotStyle -> Directive[Blue, Opacity[0.10]]],
ParametricPlot3D[RA[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Cyan, Opacity[0.15]]],
ParametricPlot3D[RI[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Red, Opacity[0.25]]],
ParametricPlot3D[RG[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Red, Opacity[0.35]]]];
BLKS:=Grid[{{horizons[a, 35, 0, 0], horizons[0, 35, 0, 0]}}];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 6) FUNKTIONEN ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
j[v_]:=Sqrt[1+μ v^2];                                                    (* Lorentzfaktor *)
mirr=Sqrt[2-℧^2+2 Sqrt[1-a^2-℧^2]]/2;                                (* irreduzible Masse *)
я=Sqrt[Χ/Σ]Sin[θ[τ]];                                            (* axialer Umfangsradius *)
яi[τ_]:=Sqrt[Χi[τ]/Σi[τ]]Sin[Θ[τ]];
Ы=Sqrt[χ/Ξ]Sin[θ0];

Σ=r[τ]^2+a^2 Cos[θ[τ]]^2;                                    (* poloidialer Umfangsradius *)
Σi[τ_]:=R[τ]^2+a^2 Cos[Θ[τ]]^2;
Ξ=r0^2+a^2 Cos[θ0]^2;

Δ=r[τ]^2-2r[τ]+a^2+℧^2;
Δr[r_]:=r^2-2r+a^2+℧^2;
Δi[τ_]:=R[τ]^2-2R[τ]+a^2+℧^2;
δ=r0^2-2r0+a^2+℧^2;

Χ=(r[τ]^2+a^2)^2-a^2 Sin[θ[τ]]^2 Δ;
Χi[τ_]:=(R[τ]^2+a^2)^2-a^2 Sin[Θ[τ]]^2 Δi[τ];
χ=(r0^2+a^2)^2-a^2 Sin[θ0]^2 δ;

Xj=a Sin[θ0]^2;
xJ[τ_]:=a Sin[Θ[τ]]^2;
XJ=a Sin[θ[τ]]^2;

Pr[r_]:=Abs[ε](r^2+a^2)+℧ q r-a Lz;
Pt[τ_]:=ε(R[τ]^2+a^2)+℧ q R[τ]-a Lz;
Pτ=ε(r[τ]^2+a^2)+℧ q r[τ]-a Lz;

Vr[r_]:=Pr[r]^2-Δr[r](μ^2 r^2+(Lz-a ε)^2+Q);             (* effektives radiales Potential *)
Vθ[θ_]:=Q-Cos[θ]^2(a^2 (μ^2-ε^2)+Lz^2 Sin[θ]^(-2)); (* effektives latitudinales Potential *)

т[τ_]:=Evaluate[t[τ]/.sol][[1]];                        (* Koordinatenzeit nach Eigenzeit *)
д[ξ_]:=Quiet[zt /.FindRoot[т[zt]-ξ, {zt, 0}]];          (* Eigenzeit nach Koordinatenzeit *)
T :=Quiet[д[tk]];                       
 
ю[τ_]:=Evaluate[t'[τ]/.sol][[1]];
γ[τ_]:=If[μ==0, "Infinity", ю[τ]];                                           (* totale ZD *)
R[τ_]:=Evaluate[r[τ]/.sol][[1]];                                (* Boyer-Lindquist Radius *)
Φ[τ_]:=Evaluate[φ[τ]/.sol][[1]];                           
Θ[τ_]:=Evaluate[θ[τ]/.sol][[1]];
ß[τ_]:=Sqrt[X'[τ]^2+Y'[τ]^2+Z'[τ]^2 ]/ю[τ];
 
ς[τ_]:=Sqrt[Χi[τ]/Δi[τ]/Σi[τ]]; ς0=Sqrt[χ/δ/Ξ];                         (* gravitative ZD *)
ω[τ_]:=(a(2R[τ]-℧^2))/Χi[τ]; ω0=(a(2r0-℧^2))/χ; ωζ=(a(2r[τ]-℧^2))/Χ;    (* F-Drag Winkelg *)
Ω[τ_]:=ω[τ] Sqrt[X[τ]^2+Y[τ]^2];            (* Frame Dragging beobachtete Geschwindigkeit *)
й[τ_]:=ω[τ] яi[τ] ς[τ]; й0=ω0 Ы ς0;              (* Frame Dragging lokale Geschwindigkeit *)
 
ж[τ_]:=Sqrt[ς[τ]^2-1]/ς[τ]; ж0=Sqrt[ς0^2-1]/ς0;                  (* Fluchtgeschwindigkeit *)
V[τ_]:=If[μ==0, 1, Re[Sqrt[-ς[τ]^2+ю[τ]^2]/ю[τ]]];
                                                  (* Fluchtgeschwindigkeit von r0 nach r1 *)
vd1:=v1/.NSolve[Sqrt[δ Ξ/χ]/Sqrt[1-v1^2]+((q r0 ℧)/(r0^2+a^2 Cos[θ0]^2))==Sqrt[((a^2+
(-2+r1) r1+℧^2) (r1^2+a^2 Cos[θ0]^2))/((a^2+r1^2)^2-a^2 (a^2+(-2+r1) r1+℧^2) Sin[θ0]^2)]+
(a^2 q r1 ℧ (2 r1-℧^2) Sin[θ0]^2)/((r1^2+a^2 Cos[θ0]^2) ((a^2+r1^2)^2-a^2 (a^2+(-2+r1) r1+
℧^2) Sin[θ0]^2))+((q r1 ℧)/(r1^2+a^2 Cos[θ0]^2))&&v1>0,v1][[1]];
                                                          (* lokale Dreiergeschwindigkeit *)
vd[τ_]:=Abs[(Sqrt[Δ Σ^3 Χ-ε^2 Σ^2 Χ^2-2 a Lz ε Σ^2 Χ ℧^2-a^2 Lz^2 Σ^2 ℧^4+
4 a Lz ε Σ^2 Χ r[τ]+2 q ε Σ Χ^2 ℧ r[τ]+4 a^2 Lz^2 Σ^2 ℧^2 r[τ]+2 a Lz q Σ Χ ℧^3 r[τ]-
4 a^2 Lz^2 Σ^2 r[τ]^2-4 a Lz q Σ Χ ℧ r[τ]^2-q^2 Χ^2 ℧^2 r[τ]^2])/(ε Σ Χ+
a Lz Σ ℧^2-2 a Lz Σ r[τ]-q Χ ℧ r[τ])];

v[τ_]:=If[μ==0, 1, Evaluate[vlt'[τ]/.sol][[1]]];

ν:=If[μ==0, 1, Re[Sqrt[(Δ Σ-Χ(εζ-Lζ ωζ)^2)/(μ Χ (εζ-Lζ ωζ)^2)]]];
vesc[τ_]:=Abs[(\[Sqrt](R[τ]^2 (R[τ]^2 (Δi[τ] Σi[τ]-Χi[τ])+2 q R[τ] Χi[τ] ℧-q^2 Χi[τ] ℧^2)+
2 a^2 R[τ] (R[τ] Δi[τ] Σi[τ]-R[τ] Χi[τ]+q Χi[τ] ℧) Cos[Θ[τ]]^2+a^4 (Δi[τ] Σi[τ]-
Χi[τ]) Cos[Θ[τ]]^4))/(Sqrt[Χi[τ]] (R[τ] (R[τ]-q ℧)+a^2 Cos[Θ[τ]]^2))];

dst[τ_]:=Evaluate[str[τ]/.sol][[1]];                                           (* Strecke *)
 
pΘ[τ_]:=Evaluate[Ξ θ'[τ] /. sol][[1]];
pR[τ_]:=Evaluate[r'[τ] Ξ/δ /. sol][[1]];

epot[τ_]:=ε+μ-ekin[τ];                                             (* potentielle Energie *)
ekin[τ_]:=If[μ==0, ς[τ], 1/Sqrt[1-v[τ]^2]-1];                       (* kinetische Energie *)
 
                                                               (* beobachtete Inklination *)
ink0:=б/. Solve[Z'[0]/ю[0] Cos[б]==-Y'[0]/ю[0] Sin[б]&&б>0&&б<2π&&б<δp[r0, a], б][[1]];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 7) DIFFERENTIALGLEICHUNG |||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
dp= \!\(\*SuperscriptBox[\(Y\),\(Y\)]\); n0[z_]:=Chop[Re[N[Chop[Simplify[Chop[z]]]]]];
 
Φ10=(a(2r0-℧^2) ε+(Ξ-2r0)Lz Csc[θ0]^2)/(Ξ δ);
Θ10=pθ0/Ξ;
r10=pr0 δ/Ξ;
t10=-((a (-℧^2+2 r0) Sin[θ0]^2 Φ10)/(Ξ+℧^2-2 r0))+\[Sqrt]((1/(δ (Ξ+℧^2-2 r0)^2))(Ξ (Ξ+℧^2-
2 r0) (-δ μ+Ξ r10^2+δ Ξ Θ10^2)-δ χ (-Ξ-℧^2+2 r0) Sin[θ0]^2 Φ10^2+a^2 δ (℧^2-
2 r0)^2 Sin[θ0]^4 Φ10^2));
                                       
DGL={
 
t''[τ]==-(((r'[τ] ((a^2+r[τ]^2) (a^2 Cos[θ[τ]]^2 (q ℧-2 t'[τ])+r[τ] (-q ℧ r[τ]+
2 (-℧^2+r[τ]) t'[τ]))+a (2 a^4 Cos[θ[τ]]^2+a^2 ℧^2 (3+Cos[2 θ[τ]]) r[τ]-
a^2 (3+Cos[2 θ[τ]]) r[τ]^2+4 ℧^2 r[τ]^3-6 r[τ]^4) Sin[θ[τ]]^2 φ'[τ]))/(a^2+℧^2+(-2+
r[τ]) r[τ])+a^2 θ'[τ] (Sin[2 θ[τ]] (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])-2 a Cos[θ[τ]] (℧^2-
2 r[τ]) Sin[θ[τ]]^3 φ'[τ]))/(a^2 Cos[θ[τ]]^2+r[τ]^2)^2),
 
t'[0]==1/(δ Ξ Sin[θ0]^2) (Lz (δ Xj-a Sin[θ0]^2 (r0^2+a^2))+ε (-δ Xj^2+
Sin[θ0]^2 (r0^2+a^2)^2)-q ℧ r0 Sin[θ0]^2 (r0^2+a^2)),
t[0]==0,
 
r''[τ]==((-1+r[τ])/(a^2+℧^2+(-2+r[τ]) r[τ])-r[τ]/(a^2 Cos[θ[τ]]^2+r[τ]^2)) r'[τ]^2+
(a^2 Sin[2 θ[τ]] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(8 (a^2 Cos[θ[τ]]^2+
r[τ]^2)^3))(a^2+℧^2+(-2+r[τ]) r[τ]) (8 t'[τ] (a^2 Cos[θ[τ]]^2 (-q ℧+t'[τ])+
r[τ] (q ℧ r[τ]+(℧^2-r[τ]) t'[τ]))+8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2+
8 a Sin[θ[τ]]^2 (a^2 Cos[θ[τ]]^2 (q ℧-2 t'[τ])+r[τ] (-q ℧ r[τ]+2 (-℧^2+r[τ]) t'[τ])) φ'[τ]+
Sin[θ[τ]]^2 (r[τ] (a^2 (3 a^2+4 ℧^2+4 (a-℧) (a+℧) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]])+
8 r[τ] (2 a^2 Cos[θ[τ]]^2 r[τ]+r[τ]^3-a^2 Sin[θ[τ]]^2))+2 a^4 Sin[2 θ[τ]]^2) φ'[τ]^2),
 
r'[0]==(pr0 δ)/Ξ,
r[0]==r0,
 
θ''[τ]==-((a^2 Cos[θ[τ]] Sin[θ[τ]] r'[τ]^2)/((a^2+℧^2+(-2+r[τ]) r[τ]) (a^2 Cos[θ[τ]]^2+
r[τ]^2)))-(2 r[τ] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(16 (a^2 Cos[θ[τ]]^2+
r[τ]^2)^3))Sin[2 θ[τ]] (a^2 (-8 t'[τ] (2 q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])+8 (a^2 Cos[θ[τ]]^2+
r[τ]^2)^2 θ'[τ]^2)+16 a (a^2+r[τ]^2) (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ]) φ'[τ]+(3 a^6-5 a^4 ℧^2+
10 a^4 r[τ]+11 a^4 r[τ]^2-8 a^2 ℧^2 r[τ]^2+16 a^2 r[τ]^3+16 a^2 r[τ]^4+8 r[τ]^6+
a^4 Cos[4 θ[τ]] (a^2+℧^2+(-2+r[τ]) r[τ])+4 a^2 Cos[2 θ[τ]] (a^2+℧^2+(-2+
r[τ]) r[τ]) (a^2+2 r[τ]^2)) φ'[τ]^2),
 
θ'[0]==pθ0/Ξ,
θ[0]==θ0,
 
φ''[τ]==-(1/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2))((r'[τ] (4 a q ℧ (a^2 Cos[θ[τ]]^2-r[τ]^2)-
8 a (a^2 Cos[θ[τ]]^2+(℧^2-r[τ]) r[τ]) t'[τ]+(a^2 (3 a^2+8 ℧^2+a^2 (4 Cos[2 θ[τ]]+
Cos[4 θ[τ]])) r[τ]-4 a^2 (3+Cos[2 θ[τ]]) r[τ]^2+8 (a^2+℧^2+a^2 Cos[2 θ[τ]]) r[τ]^3-
16 r[τ]^4+8 r[τ]^5+2 a^4 Sin[2 θ[τ]]^2) φ'[τ]))/(a^2+℧^2+(-2+r[τ]) r[τ])+
θ'[τ] (8 a Cot[θ[τ]] (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])+(8 Cot[θ[τ]] (a^2+r[τ]^2)^2-
2 a^2 (3 a^2+2 ℧^2+4 (-1+r[τ]) r[τ]) Sin[2 θ[τ]]-a^4 Sin[4 θ[τ]]) φ'[τ])),
 
φ'[0]==1/(δ Ξ Sin[θ0]^2) (ε (-δ Xj+a Sin[θ0]^2 (r0^2+a^2))+Lz (δ-a^2 Sin[θ0]^2)-
q ℧ r0 a Sin[θ0]^2),
φ[0]==φ0,
 
str'[τ]==If[μ==0, 1, vd[τ]/Max[1*^-16, Abs[Sqrt[1-vd[τ]^2]]]],
str[0]==0,
vlt'[τ]==If[μ==0, 1, vd[τ]],
vlt[0]==0
 
};
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 8) INTEGRATION |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
sol=NDSolve[DGL, {t, r, θ, φ, vlt, str}, {τ, 0, tmax+1/1000},
WorkingPrecision-> wp,
MaxSteps-> Infinity,
Method-> mta,
InterpolationOrder-> All,
StepMonitor :> (laststep=plunge; plunge=τ;
stepsize=plunge-laststep;), Method->{"EventLocator",
"Event" :> (If[stepsize<1*^-4, 0, 1])}];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 9) KOORDINATEN |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
X[τ_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]];            (* kartesisch *)
Y[τ_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]];
Z[τ_]:=Evaluate[r[τ] Cos[θ[τ]]/.sol][[1]];
 
x[τ_]:=Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]];       (* Plotkoordinaten *)
y[τ_]:=Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]];
z[τ_]:=Z[τ];
 
XYZ[τ_]:=Sqrt[X[τ]^2+Y[τ]^2+Z[τ]^2]; XY[τ_]:=Sqrt[X[τ]^2+Y[τ]^2];  (* kartesischer Radius *)
 
Xyz[{x_, y_, z_}, α_]:={x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z};      (* Rotationsmatrix *)
xYz[{x_, y_, z_}, β_]:={x Cos[β]+z Sin[β], y, z Cos[β]-x Sin[β]};
xyZ[{x_, y_, z_}, ψ_]:={x, y Cos[ψ]-z Sin[ψ], y Sin[ψ]+z Cos[ψ]};
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 10) PLOT EINSTELLUNGEN |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
PR=1.2 r1;                                                                  (* Plot Range *)
VP={r0, r0, r0};                                                     (* Perspektive x,y,z *)
d1=10;                                                                    (* Schweiflänge *)
plp=50;                                                            (* Flächenplot Details *)
Plp=Automatic;                                                          (* Kurven Details *)
 
w1l=0; w2l=0; w1r=0; w2r=0;                                          (* Startperspektiven *)
Mrec=100; mrec=10;                                       (* Parametric Plot Subdivisionen *)
imgsize=380;                                                                 (* Bildgröße *)
 
s[text_]:=Style[text, FontSize->font]; font=11;                            (* Anzeigestil *)
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 11) PLOT NACH KOORDINATENZEIT ||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
display[T_]:=Grid[{
{s[" t coord"], " = ", s[n0[tk]], s["GM/c³"], s[dp]},
{If[μ==0, s[" affineP"], s[" τ propr"]], " = ", s[n0[T]], s["GM/c³"], s[dp]},
{s[" ṫ total"], " = ", s[n0[γ[T]]], s["dt/dτ"], s[dp]},
{s[" ς gravt"], " = ", s[n0[ς[T]]], s["dt/dτ"], s[dp]},
{s[" γ kinet"], " = ", s[n0[1/Sqrt[1-v[T]^2]]], s["dt/dτ"], s[dp]},
{s[" R carts"], " = ", s[n0[XYZ[T]]], s["GM/c²"], s[dp]},
{s[" x carts"], " = ", s[n0[X[T]]], s["GM/c²"], s[dp]},
{s[" y carts"], " = ", s[n0[Y[T]]], s["GM/c²"], s[dp]},
{s[" z carts"], " = ", s[n0[Z[T]]], s["GM/c²"], s[dp]},
{s[" s dstnc"], " = ", s[n0[dst[T]]], s["GM/c²"], s[dp]},
{s[" r coord"], " = ", s[n0[R[T]]], s["GM/c²"], s[dp]},
{s[" φ longd"], " = ", s[n0[Φ[T] 180/π]], s["deg"], s[dp]},
{s[" θ lattd"], " = ", s[n0[Θ[T] 180/π]], s["deg"], s[dp]},
{s[" d¹r/dτ¹"], " = ", s[n0[R'[T]]], s["c"], s[dp]},
{s[" d¹φ/dτ¹"], " = ", s[n0[Φ'[T]]], s["c\.b3/G/M"], s[dp]},
{s[" d¹θ/dτ¹"], " = ", s[n0[Θ'[T]]], s["c\.b3/G/M"], s[dp]},
{s[" d\.b2r/dτ\.b2"], " = ", s[n0[R''[T]]], s["c⁴/G/M"], s[dp]},
{s[" d\.b2φ/dτ\.b2"], " = ", s[n0[Φ''[T]]], s["c⁶/G\.b2/M\.b2"], s[dp]},
{s[" d\.b2θ/dτ\.b2"], " = ", s[n0[Θ''[T]]], s["c⁶/G\.b2/M\.b2"], s[dp]},
{s[" a SpinP"], " = ", s[n0[a]], s["GM²/c"], s[dp]},
{s[" ℧ cntrl"], " = ", s[n0[℧]], s["Q/M"], s[dp]},
{s[" q prtcl"], " = ", s[n0[q]], s["q/m"], s[dp]},
{s[" M irred"], " = ", s[n0[mirr]], s["M"], s[dp]},
{s[" E kinet"], " = ", s[n0[ekin[T]]], s["mc²"], s[dp]},
{s[" E poten"], " = ", s[n0[epot[T]]], s["mc²"], s[dp]},
{s[" E total"], " = ", s[n0[ε]], s["mc²"], s[dp]},
{s[" CarterQ"], " = ", s[n0[Qk]], s["GMm/c"], s[dp]},
{s[" L axial"], " = ", s[n0[Lz]], s["GMm/c"], s[dp]},
{s[" L polar"], " = ", s[n0[pΘ[T]]], s["GMm/c"], s[dp]},
{s[" p r.mom"], " = ", s[n0[pR[T]]], s["mc"], s[dp]},
{s[" ω fdrag"], " = ", s[n0[Abs[ω[T]]]], s["c³/G/M"], s[dp]},
{s[" v fdrag"], " = ", s[n0[Abs[й[T]]]], s["c"], s[dp]},
{s[" Ω fdrag"], " = ", s[n0[Abs[Ω[T]]]], s["c"], s[dp]},
{s[" v propr"], " = ", s[n0[v[T]/Sqrt[1-v[T]^2]]], s["c"], s[dp]},
{s[" v escpe"], " = ", s[n0[vesc[T]]], s["c"], s[dp]},
{s[" v obsvd"], " = ", s[n0[ß[T]]], s["c"], s[dp]},
{s[" v r,loc"], " = ", s[n0[vr[T]]], s["c"], s[dp]},
{s[" v θ,loc"], " = ", s[n0[vθ[T]]], s["c"], s[dp]},
{s[" v φ,loc"], " = ", s[n0[vφ[T]]], s["c"], s[dp]},
{s[" v local"], " = ", s[n0[v[T]]], s["c"], s[dp]},
{s[" "], s[" "], s["                   "], s["         "]}},
Alignment-> Left, Spacings-> {0, 0}];
 
plot1a[{xx_, yy_, zz_, tk_, w1_, w2_}]:=                                     (* Animation *)
Show[Graphics3D[{
{PointSize[0.009], Red, Point[
Xyz[xyZ[{x[T], y[T], z[T]}, w1], w2]]}},
ImageSize-> imgsize,
PlotRange-> PR,
SphericalRegion->False,
ImagePadding-> 1],
horizons[A, None, w1, w2],
If[a==0, {},
Graphics3D[{{PointSize[0.009], Purple, Point[
Xyz[xyZ[{
Sin[-φ0-ω0 tk+π/2] Sqrt[x0[A]^2+y0[A]^2],
Cos[-φ0-ω0 tk+π/2] Sqrt[x0[A]^2+y0[A]^2],
z0[A]}, w1], w2]]}}]],
If[tk==0, {}, If[a==0, {},
ParametricPlot3D[
Xyz[xyZ[{
Sin[-φ0-ω0 tt+π/2] Sqrt[x0[A]^2+y0[A]^2],
Cos[-φ0-ω0 tt+π/2] Sqrt[x0[A]^2+y0[A]^2],
z0[A]}, w1], w2],
{tt, Max[0, tk-199/100 π/ω0], tk},
PlotStyle -> {Thickness[0.001], Dashed, Purple},
PlotPoints-> Automatic,
MaxRecursion-> mrec]]],
Block[{$RecursionLimit = Mrec},
If[tk==0, {},
ParametricPlot3D[
Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2], {tt, If[TMax<0, Min[0, T+d1], Max[0, T-d1]], T},
PlotStyle-> {Thickness[0.004]},
ColorFunction-> Function[{x, y, z, t},
Hue[0, 1, 0.5, If[TMax<0, Max[Min[(+T+(-t+d1))/d1, 1], 0], Max[Min[(-T+(t+d1))/d1, 1], 0]]]],
ColorFunctionScaling-> False,
PlotPoints-> Automatic,
MaxRecursion-> mrec]]],
If[tk==0, {},
Block[{$RecursionLimit = Mrec},
ParametricPlot3D[
Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2],
{tt, 0, If[Tmax<0, Min[-1*^-16, T+d1/3], Max[1*^-16, T-d1/3]]},
PlotStyle-> {Thickness[0.004], Opacity[0.6], Darker[Gray]},
PlotPoints-> Plp,
MaxRecursion-> mrec]]],
ViewPoint-> {xx, yy, zz}];
 
Quiet[Do[
Print[Rasterize[Grid[{{
plot1a[{0, -Infinity, 0, tk, w1l, w2l}],
plot1a[{0, 0, Infinity, tk, w1r, w2r}],
display[Quiet[д[tk]]]
}}, Alignment->Left]]],
{tk, 0, TMax, TMax/2}]]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 12) PLOT NACH EIGENZEIT ||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
display[T_]:=Grid[{
{If[μ==0, s[" affineP"], s[" τ propr"]], " = ", s[n0[tp]], s["GM/c³"], s[dp]},
{s[" t coord"], " = ", s[n0[т[tp]]], s["GM/c³"], s[dp]},
{s[" ṫ total"], " = ", s[n0[γ[tp]]], s["dt/dτ"], s[dp]},
{s[" ς gravt"], " = ", s[n0[ς[tp]]], s["dt/dτ"], s[dp]},
{s[" γ kinet"], " = ", s[n0[1/Sqrt[1-v[tp]^2]]], s["dt/dτ"], s[dp]},
{s[" R carts"], " = ", s[n0[XYZ[tp]]], s["GM/c²"], s[dp]},
{s[" x carts"], " = ", s[n0[X[tp]]], s["GM/c²"], s[dp]},
{s[" y carts"], " = ", s[n0[Y[tp]]], s["GM/c²"], s[dp]},
{s[" z carts"], " = ", s[n0[Z[tp]]], s["GM/c²"], s[dp]},
{s[" s dstnc"], " = ", s[n0[dst[tp]]], s["GM/c²"], s[dp]},
{s[" r coord"], " = ", s[n0[R[tp]]], s["GM/c²"], s[dp]},
{s[" φ longd"], " = ", s[n0[Φ[tp] 180/π]], s["deg"], s[dp]},
{s[" θ lattd"], " = ", s[n0[Θ[tp] 180/π]], s["deg"], s[dp]},
{s[" d¹r/dτ¹"], " = ", s[n0[R'[tp]]], s["c"], s[dp]},
{s[" d¹φ/dτ¹"], " = ", s[n0[Φ'[tp]]], s["c\.b3/G/M"], s[dp]},
{s[" d¹θ/dτ¹"], " = ", s[n0[Θ'[tp]]], s["c\.b3/G/M"], s[dp]},
{s[" d\.b2r/dτ\.b2"], " = ", s[n0[R''[tp]]], s["c⁴/G/M"], s[dp]},
{s[" d\.b2φ/dτ\.b2"], " = ", s[n0[Φ''[tp]]], s["c⁶/G\.b2/M\.b2"], s[dp]},
{s[" d\.b2θ/dτ\.b2"], " = ", s[n0[Θ''[tp]]], s["c⁶/G\.b2/M\.b2"], s[dp]},
{s[" a SpinP"], " = ", s[n0[a]], s["GM²/c"], s[dp]},
{s[" ℧ cntrl"], " = ", s[n0[℧]], s["Q/M"], s[dp]},
{s[" q prtcl"], " = ", s[n0[q]], s["q/m"], s[dp]},
{s[" M irred"], " = ", s[n0[mirr]], s["M"], s[dp]},
{s[" E kinet"], " = ", s[n0[ekin[tp]]], s["mc²"], s[dp]},
{s[" E poten"], " = ", s[n0[epot[tp]]], s["mc²"], s[dp]},
{s[" E total"], " = ", s[n0[ε]], s["mc²"], s[dp]},
{s[" CarterQ"], " = ", s[n0[Qk]], s["GMm/c"], s[dp]},
{s[" L axial"], " = ", s[n0[Lz]], s["GMm/c"], s[dp]},
{s[" L polar"], " = ", s[n0[pΘ[tp]]], s["GMm/c"], s[dp]},
{s[" p r.mom"], " = ", s[n0[pR[tp]]], s["mc"], s[dp]},
{s[" ω fdrag"], " = ", s[n0[ω[tp]]], s["c³/G/M"], s[dp]},
{s[" v fdrag"], " = ", s[n0[й[tp]]], s["c"], s[dp]},
{s[" Ω fdrag"], " = ", s[n0[Ω[tp]]], s["c"], s[dp]},
{s[" v propr"], " = ", s[n0[v[tp]/Sqrt[1-v[tp]^2]]], s["c"], s[dp]},
{s[" v escpe"], " = ", s[n0[vesc[tp]]], s["c"], s[dp]},
{s[" v obsvd"], " = ", s[n0[ß[tp]]], s["c"], s[dp]},
{s[" v r,loc"], " = ", s[n0[vr[tp]]], s["c"], s[dp]},
{s[" v θ,loc"], " = ", s[n0[vθ[tp]]], s["c"], s[dp]},
{s[" v φ,loc"], " = ", s[n0[vφ[tp]]], s["c"], s[dp]},
{s[" v local"], " = ", s[n0[v[tp]]], s["c"], s[dp]},
{s[" "], s[" "], s["                   "], s["         "]}},
Alignment-> Left, Spacings-> {0, 0}];
 
plot1b[{xx_, yy_, zz_, tk_, w1_, w2_}]:=                                     (* Animation *)
Show[Graphics3D[{
{PointSize[0.009], Red, Point[
Xyz[xyZ[{x[tp], y[tp], z[tp]}, w1], w2]]}},
ImageSize-> imgsize,
PlotRange-> PR,
SphericalRegion->False,
ImagePadding-> 1],
horizons[A, None, w1, w2],
If[a==0, {},
Graphics3D[{{PointSize[0.009], Purple, Point[
Xyz[xyZ[{
Sin[-φ0-ω0 т[tp]+π/2] Sqrt[x0[A]^2+y0[A]^2],
Cos[-φ0-ω0 т[tp]+π/2] Sqrt[x0[A]^2+y0[A]^2],
z0[A]}, w1], w2]]}}]],
If[tk==0, {}, If[a==0, {},
ParametricPlot3D[
Xyz[xyZ[{
Sin[-φ0-ω0 т[tt]+π/2] Sqrt[x0[A]^2+y0[A]^2],
Cos[-φ0-ω0 т[tt]+π/2] Sqrt[x0[A]^2+y0[A]^2],
z0[A]}, w1], w2],
{tt, Max[0, д[т[tp]-199/100 π/ω0]], tp},
PlotStyle -> {Thickness[0.001], Dashed, Purple},
PlotPoints-> Automatic,
MaxRecursion-> 12]]],
If[tk==0, {},
Block[{$RecursionLimit = Mrec},
ParametricPlot3D[
Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2], {tt, If[tp<0, Min[0, tp+d1], Max[0, tp-d1]], tp},
PlotStyle-> {Thickness[0.004]},
ColorFunction-> Function[{x, y, z, t},
Hue[0, 1, 0.5, If[tp<0,
Max[Min[(+tp+(-t+d1))/d1, 1], 0], Max[Min[(-tp+(t+d1))/d1, 1], 0]]]],
ColorFunctionScaling-> False,
PlotPoints-> Automatic,
MaxRecursion-> mrec]]],
If[tk==0, {},
Block[{$RecursionLimit = Mrec},
ParametricPlot3D[
Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2],
{tt, 0, If[tp<0, Min[-1*^-16, tp+d1/3], Max[1*^-16, tp-d1/3]]},
PlotStyle-> {Thickness[0.004], Opacity[0.6], Darker[Gray]},
PlotPoints-> Plp,
MaxRecursion-> mrec]]],
ViewPoint-> {xx, yy, zz}];
 
Do[
Print[Rasterize[Grid[{{
plot1b[{0, -Infinity, 0, tp, w1l, w2l}],
plot1b[{0, 0, +Infinity, tp, w1r, w2r}],
display[tp]
}}, Alignment->Left]]],
{tp, 0, tMax, tMax/2}]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 13) EXPORTOPTIONEN |||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
(* Export als HTML Dokument *)
(* Export["Y:\\export\\dateiname.html", EvaluationNotebook[], "GraphicsOutput" -> "PNG"]  *)
(* Export direkt als Bildsequenz *)
(* Do[Export["Y:\\export\\dateiname" <> ToString[tk] <> ".png", Rasterize[...]            *)
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||| http://kerr.newman.yukerez.net ||||| Simon Tyran, Vienna |||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)











Erhaltungsgrößen Testmonitor zur Probe der numerischen Stabilität

Code: Alles auswählen

N@ε

ε0=With[{τ=0},Evaluate[+((q r[τ] ℧)/(r[τ]^2+a^2 Cos[θ[τ]]^2))+t'[τ] (1-(2 r[τ]-℧^2)/(r[τ]^2+a^2 Cos[θ[τ]]^2))+(a φ'[τ] (2 r[τ]-℧^2) Sin[θ[τ]]^2)/(r[τ]^2+a^2 Cos[θ[τ]]^2)/.sol][[1]]]

ε1=With[{τ=tMax 99/100},Evaluate[+((q r[τ] ℧)/(r[τ]^2+a^2 Cos[θ[τ]]^2))+t'[τ] (1-(2 r[τ]-℧^2)/(r[τ]^2+a^2 Cos[θ[τ]]^2))+(a φ'[τ] (2 r[τ]-℧^2) Sin[θ[τ]]^2)/(r[τ]^2+a^2 Cos[θ[τ]]^2)/.sol][[1]]]

N@Lz

Lz0=With[{τ=0},Evaluate[(q a r [τ]℧ Sin[θ[τ]]^2)/(r[τ]^2+a^2 Cos[θ[τ]]^2)-(a t'[τ]  (2 r[τ]-℧^2) Sin[θ[τ]]^2)/(r[τ]^2+a^2 Cos[θ[τ]]^2)+(φ'[τ] Sin[θ[τ]]^2 ((a^2+r[τ]^2)^2-a^2 (a^2-2 r[τ]+r[τ]^2+℧^2) Sin[θ[τ]]^2))/(r[τ]^2+a^2 Cos[θ[τ]]^2)/.sol][[1]]]

Lz1=With[{τ=tMax 99/100},Evaluate[(q a r [τ]℧ Sin[θ[τ]]^2)/(r[τ]^2+a^2 Cos[θ[τ]]^2)-(a t'[τ]  (2 r[τ]-℧^2) Sin[θ[τ]]^2)/(r[τ]^2+a^2 Cos[θ[τ]]^2)+(φ'[τ] Sin[θ[τ]]^2 ((a^2+r[τ]^2)^2-a^2 (a^2-2 r[τ]+r[τ]^2+℧^2) Sin[θ[τ]]^2))/(r[τ]^2+a^2 Cos[θ[τ]]^2)/.sol][[1]]]











Startbedingungssolver: Umrechung von lokaler Geschwindigkeit, Eigenzeitableitungen und Erhaltungsgrößen

Code: Alles auswählen

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE 1: Lokale Geschwindigkeit nach Erhaltungsgrößen ε, Lz, Q ||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]; ClearAll["Local`*"];
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 5;
θ0 = π/4;
φ0 = 0;
a  = 2/3;
℧  = 1/2;
q  = -1/3;
μ  =-1;
 
(* || Erhaltungsgrößen Gesamtenergie, axialer Drehimpuls & Carter Konstante eingeben |||| *)
 
ε  = -(785010/23679959)+(27 Sqrt[13/23679959])/2+2 Sqrt[256510/1356121];
Lz = (-195+Sqrt[307839467])/17706;
Q  = (390497414585450885-15093984745944 Sqrt[1130]-162128393295 Sqrt[307839467]+2851156320 Sqrt[347858597710])/131213267228553354;
 
(* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *)
 
ε0 = (q r0 ℧)/(r0^2+a^2 Cos[θ0]^2)+Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]/Sqrt[1+v0^2 μ]+(a (2 r0-℧^2) ((a q r0 ℧ Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ]))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2);

L0 = (a q r0 ℧ Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ];

Q0 = (vθ0^2 (r0^2+a^2 Cos[θ0]^2))/(1+v0^2 μ)+Cos[θ0]^2 ((a q r0 ℧ Sin[θ0])/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ])^2-a^2 Cos[θ0]^2 (μ+((q r0 ℧)/(r0^2+a^2 Cos[θ0]^2)+Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]/Sqrt[1+v0^2 μ]+(a (2 r0-℧^2) Sin[θ0] ((a q r0 ℧ Sin[θ0])/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ]))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))^2);
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code 1"
Reduce[ε0==ε && L0==Lz && Q0==Q && vr0^2==v0^2-vφ0^2-vθ0^2 && v0>0, {v0,vr0,vφ0,vθ0}]
N[%]

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax |||| kerr.newman.yukterez.net |||| Simon Tyran, Vienna  ||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                                                   
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE 2: Lokale Geschwindigkeit nach ersten Eigenzeitableitungen  |||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]; ClearAll["Local`*"];
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 5;
θ0 = π/4;
φ0 = 0;
a  = 2/3;
℧  = 1/2;
q  = -1/3;
μ  =-1;
 
(* || Startwerte für die ersten Eigenzeitableitungen eingeben ||||||||||||||||||||||||||| *)
 
dt = -(1053/5822777)+4 Sqrt[208634/1667315];
dr = Sqrt[565/2951]/2;
dθ = 3/Sqrt[2951];
dφ = -(1108809/607414628309)+324 Sqrt[26/13379176835]+6 Sqrt[227/1356121];
 
(* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *)
 
ε0 = (q r0 ℧)/(r0^2+a^2 Cos[θ0]^2)+Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]/Sqrt[1+v0^2 μ]+(a (2 r0-℧^2) ((a q r0 ℧ Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ]))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2);

L0 = (a q r0 ℧ Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ];

Q0 = (vθ0^2 (r0^2+a^2 Cos[θ0]^2))/(1+v0^2 μ)+Cos[θ0]^2 ((a q r0 ℧ Sin[θ0])/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ])^2-a^2 Cos[θ0]^2 (μ+((q r0 ℧)/(r0^2+a^2 Cos[θ0]^2)+Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]/Sqrt[1+v0^2 μ]+(a (2 r0-℧^2) Sin[θ0] ((a q r0 ℧ Sin[θ0])/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ]))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))^2);
 
(* || Benötigte Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
Ξ=r0^2+a^2 Cos[θ0]^2;
δ=r0^2-2r0+a^2+℧^2;
χ=(r0^2+a^2)^2-a^2 Sin[θ0]^2 δ;
Xj=a Sin[θ0]^2;
j[v_]:=Sqrt[1+μ v^2];
pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];
pθ0=vθ0 Sqrt[Ξ]/j[v0];
 
dT=1/(δ Ξ Sin[θ0]^2) (L0 (δ Xj-a Sin[θ0]^2 (r0^2+a^2))+ε0 (-δ Xj^2+Sin[θ0]^2 (r0^2+a^2)^2)-q ℧ r0 Sin[θ0]^2 (r0^2+a^2));
dR=(pr0 δ)/Ξ;
dΘ=pθ0/Ξ;
dΦ=1/(δ Ξ Sin[θ0]^2) (ε0 (-δ Xj+a Sin[θ0]^2 (r0^2+a^2))+L0 (δ-a^2 Sin[θ0]^2)-q ℧ r0 a Sin[θ0]^2);
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code 2"
Reduce[dT==dt && dR==dr && dΘ==dθ && dΦ==dφ && v0>0, {v0,vr0,vφ0,vθ0}]
N[%]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax |||| kerr.newman.yukterez.net |||| Simon Tyran, Vienna  ||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                                                   
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE 3: Erste Eigenzeitableitungen nach Erhaltungsgrößen ε, Lz, Q ||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]; ClearAll["Local`*"];
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 5;
θ0 = π/4;
φ0 = 0;
a  = 2/3;
℧  = 1/2;
q  = -1/3;
μ  =-1;
 
(* || Erhaltungsgrößen Gesamtenergie, axialer Drehimpuls & Carter Konstante eingeben |||| *)
 
ε  = -(785010/23679959)+(27 Sqrt[13/23679959])/2+2 Sqrt[256510/1356121];
Lz = (-195+Sqrt[307839467])/17706;
Q  = (390497414585450885-15093984745944 Sqrt[1130]-162128393295 Sqrt[307839467]+2851156320 Sqrt[347858597710])/131213267228553354;
 
(* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *)
 
ε0 = (q r0 ℧)/(r0^2+a^2 Cos[θ0]^2)+Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]/Sqrt[1+v0^2 μ]+(a (2 r0-℧^2) ((a q r0 ℧ Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ]))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2);

L0 = (a q r0 ℧ Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ];

Q0 = (vθ0^2 (r0^2+a^2 Cos[θ0]^2))/(1+v0^2 μ)+Cos[θ0]^2 ((a q r0 ℧ Sin[θ0])/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ])^2-a^2 Cos[θ0]^2 (μ+((q r0 ℧)/(r0^2+a^2 Cos[θ0]^2)+Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]/Sqrt[1+v0^2 μ]+(a (2 r0-℧^2) Sin[θ0] ((a q r0 ℧ Sin[θ0])/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ]))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))^2);
 
(* || Benötigte Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
Ξ=r0^2+a^2 Cos[θ0]^2;
δ=r0^2-2r0+a^2+℧^2;
χ=(r0^2+a^2)^2-a^2 Sin[θ0]^2 δ;
Xj=a Sin[θ0]^2;
j[v_]:=Sqrt[1+μ v^2];
pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];
pθ0=vθ0 Sqrt[Ξ]/j[v0];
 
dT=1/(δ Ξ Sin[θ0]^2) (L0 (δ Xj-a Sin[θ0]^2 (r0^2+a^2))+ε0 (-δ Xj^2+Sin[θ0]^2 (r0^2+a^2)^2)-q ℧ r0 Sin[θ0]^2 (r0^2+a^2));
dR=(pr0 δ)/Ξ;
dΘ=pθ0/Ξ;
dΦ=1/(δ Ξ Sin[θ0]^2) (ε0 (-δ Xj+a Sin[θ0]^2 (r0^2+a^2))+L0 (δ-a^2 Sin[θ0]^2)-q ℧ r0 a Sin[θ0]^2);
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code 3"
Solve[ε==ε0 && Lz==L0 && Q==Q0 && dt==dT && dr==dR && dθ==dΘ && dφ==dΦ && v0^2==vr0^2+vθ0^2+vφ0^2, {dt,dr,dθ,dφ}]
N[%]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax |||| kerr.newman.yukterez.net |||| Simon Tyran, Vienna  ||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                                                   
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE 4: Erste Eigenzeitableitungen nach lokalen Geschwindigkeiten ||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]; ClearAll["Local`*"];
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 5;
θ0 = π/4;
φ0 = 0;
a  = 2/3;
℧  = 1/2;
q  = -1/3;
μ  =-1;
 
(* || Startwerte für die ersten Eigenzeitableitungen eingeben ||||||||||||||||||||||||||| *)
 
vr0 = 1/4;
vθ0 = 1/4;
vφ0 = 1/4;
v0  = Sqrt[vr0^2+vθ0^2+vφ0^2];
 
(* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *)
 
ε0 = (q r0 ℧)/(r0^2+a^2 Cos[θ0]^2)+Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]/Sqrt[1+v0^2 μ]+(a (2 r0-℧^2) ((a q r0 ℧ Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ]))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2);

L0 = (a q r0 ℧ Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ];

Q0 = (vθ0^2 (r0^2+a^2 Cos[θ0]^2))/(1+v0^2 μ)+Cos[θ0]^2 ((a q r0 ℧ Sin[θ0])/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ])^2-a^2 Cos[θ0]^2 (μ+((q r0 ℧)/(r0^2+a^2 Cos[θ0]^2)+Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]/Sqrt[1+v0^2 μ]+(a (2 r0-℧^2) Sin[θ0] ((a q r0 ℧ Sin[θ0])/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ]))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))^2);
 
(* || Benötigte Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
Ξ=r0^2+a^2 Cos[θ0]^2;
δ=r0^2-2r0+a^2+℧^2;
χ=(r0^2+a^2)^2-a^2 Sin[θ0]^2 δ;
Xj=a Sin[θ0]^2;
j[v_]:=Sqrt[1+μ v^2];
pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];
pθ0=vθ0 Sqrt[Ξ]/j[v0];
 
dT=1/(δ Ξ Sin[θ0]^2) (L0 (δ Xj-a Sin[θ0]^2 (r0^2+a^2))+ε0 (-δ Xj^2+Sin[θ0]^2 (r0^2+a^2)^2)-q ℧ r0 Sin[θ0]^2 (r0^2+a^2));
dR=(pr0 δ)/Ξ;
dΘ=pθ0/Ξ;
dΦ=1/(δ Ξ Sin[θ0]^2) (ε0 (-δ Xj+a Sin[θ0]^2 (r0^2+a^2))+L0 (δ-a^2 Sin[θ0]^2)-q ℧ r0 a Sin[θ0]^2);
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code 4"
Reduce[dT==dt && dR==dr && dΘ==dθ && dΦ==dφ, {dt,dr,dθ,dφ}]
N[%]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax |||| kerr.newman.yukterez.net |||| Simon Tyran, Vienna  ||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                                                   
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE 5: Erhaltungsgrößen ε, Lz, Q nach lokalen Geschwindigkeiten |||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]; ClearAll["Local`*"];
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 5;
θ0 = π/4;
φ0 = 0;
a  = 2/3;
℧  = 1/2;
q  = -1/3;
μ  =-1;
 
(* || Startwerte für die ersten Eigenzeitableitungen eingeben ||||||||||||||||||||||||||| *)
 
vr0 = 1/4;
vθ0 = 1/4;
vφ0 = 1/4;
v0  = Sqrt[vr0^2+vθ0^2+vφ0^2];
 
(* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *)
 
ε0 = (q r0 ℧)/(r0^2+a^2 Cos[θ0]^2)+Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]/Sqrt[1+v0^2 μ]+(a (2 r0-℧^2) ((a q r0 ℧ Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ]))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2);

L0 = (a q r0 ℧ Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ];

Q0 = (vθ0^2 (r0^2+a^2 Cos[θ0]^2))/(1+v0^2 μ)+Cos[θ0]^2 ((a q r0 ℧ Sin[θ0])/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ])^2-a^2 Cos[θ0]^2 (μ+((q r0 ℧)/(r0^2+a^2 Cos[θ0]^2)+Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]/Sqrt[1+v0^2 μ]+(a (2 r0-℧^2) Sin[θ0] ((a q r0 ℧ Sin[θ0])/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ]))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))^2);
 
(* || Benötigte Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
Ξ=r0^2+a^2 Cos[θ0]^2;
δ=r0^2-2r0+a^2+℧^2;
χ=(r0^2+a^2)^2-a^2 Sin[θ0]^2 δ;
Xj=a Sin[θ0]^2;
j[v_]:=Sqrt[1+μ v^2];
pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];
pθ0=vθ0 Sqrt[Ξ]/j[v0];
 
dT=1/(δ Ξ Sin[θ0]^2) (L0 (δ Xj-a Sin[θ0]^2 (r0^2+a^2))+ε0 (-δ Xj^2+Sin[θ0]^2 (r0^2+a^2)^2)-q ℧ r0 Sin[θ0]^2 (r0^2+a^2));
dR=(pr0 δ)/Ξ;
dΘ=pθ0/Ξ;
dΦ=1/(δ Ξ Sin[θ0]^2) (ε0 (-δ Xj+a Sin[θ0]^2 (r0^2+a^2))+L0 (δ-a^2 Sin[θ0]^2)-q ℧ r0 a Sin[θ0]^2);
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code 5"
Reduce[ε==ε0 && Lz==L0 && Q==Q0, {ε,Lz,Q}]
N[%]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax |||| kerr.newman.yukterez.net |||| Simon Tyran, Vienna  ||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                                                   
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE 6: Erhaltungsgrößen ε, Lz, Q  nach ersten Eigenzeitableitungen  |||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]; ClearAll["Local`*"];
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 5;
θ0 = π/4;
φ0 = 0;
a  = 2/3;
℧  = 1/2;
q  = -1/3;
μ  =-1;
 
(* || Startwerte für die ersten Eigenzeitableitungen eingeben ||||||||||||||||||||||||||| *)
 
dt = -(1053/5822777)+4 Sqrt[208634/1667315];
dr = Sqrt[565/2951]/2;
dθ = 3/Sqrt[2951];
dφ = -(1108809/607414628309)+324 Sqrt[26/13379176835]+6 Sqrt[227/1356121];
 
(* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *)
 
ε0 = (q r0 ℧)/(r0^2+a^2 Cos[θ0]^2)+Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]/Sqrt[1+v0^2 μ]+(a (2 r0-℧^2) ((a q r0 ℧ Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ]))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2);

L0 = (a q r0 ℧ Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ];

Q0 = (vθ0^2 (r0^2+a^2 Cos[θ0]^2))/(1+v0^2 μ)+Cos[θ0]^2 ((a q r0 ℧ Sin[θ0])/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ])^2-a^2 Cos[θ0]^2 (μ+((q r0 ℧)/(r0^2+a^2 Cos[θ0]^2)+Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]/Sqrt[1+v0^2 μ]+(a (2 r0-℧^2) Sin[θ0] ((a q r0 ℧ Sin[θ0])/(r0^2+a^2 Cos[θ0]^2)+(vφ0 Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+v0^2 μ]))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))^2);
 
(* || Benötigte Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
Ξ=r0^2+a^2 Cos[θ0]^2;
δ=r0^2-2r0+a^2+℧^2;
χ=(r0^2+a^2)^2-a^2 Sin[θ0]^2 δ;
Xj=a Sin[θ0]^2;
j[v_]:=Sqrt[1+μ v^2];
pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];
pθ0=vθ0 Sqrt[Ξ]/j[v0];
 
dT=1/(δ Ξ Sin[θ0]^2) (L0 (δ Xj-a Sin[θ0]^2 (r0^2+a^2))+ε0 (-δ Xj^2+Sin[θ0]^2 (r0^2+a^2)^2)-q ℧ r0 Sin[θ0]^2 (r0^2+a^2));
dR=(pr0 δ)/Ξ;
dΘ=pθ0/Ξ;
dΦ=1/(δ Ξ Sin[θ0]^2) (ε0 (-δ Xj+a Sin[θ0]^2 (r0^2+a^2))+L0 (δ-a^2 Sin[θ0]^2)-q ℧ r0 a Sin[θ0]^2);
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code 6"
Solve[ε==ε0 && Lz==L0 && Q==Q0 && dT==dt && dR==dr && dΘ==dθ && dΦ==dφ, {ε,Lz,Q}]
N[%]

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax |||| kerr.newman.yukterez.net |||| Simon Tyran, Vienna  ||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)











[Yukterez]

Referenzen und vertiefende Literatur:

• Bičak, Stuchlík & Balek: Charged particle motion in the field of rotating charged black holes and naked singularities
• Stuchlík, Schee & Blaschke: Particle collisions and optical effects in the mining Kerr–Newman spacetimes
• Gal’tsov & Kobialko: Completing characterization of photon orbits in Kerr-Newman metrics
• Cebeci, Özdemir & Şentorun: Motion of charged particles in Kerr-Newman spacetime
• Calvani & Felice: Off-equatorial circular orbits in the Kerr-Newman spacetime
• Hackmann & Xu: Charged particle motion in Kerr-Newman spacetimes
• Semerák & Bičák: Interplay of forces in the Kerr–Newman field
• ICRANet: Generalizations of the Kerr-Newman solution
• Vyacheslav Dokuchaev: Voyage inside a black hole
• Newman & Adamo: Kerr-Newman metric
• Misner, Thorne & Wheeler: MTW

Animations by Simon Tyran, Vienna (Yukterez) - reuse permitted under the Creative Commons License CC BY-SA 4.0