Kerr Newman Metrik

Deutschsprachige Version
Benutzeravatar
Yukterez
Administrator
Beiträge: 207
Registriert: Mi 21. Okt 2015, 02:16

Kerr Newman Metrik

Beitragvon Yukterez » So 6. Aug 2017, 23:28

Bild

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

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

Bild  ←
Schatten eines ungefütterten Kerr-Newman SL mit a²+℧²=M², Betrachtungswinkel: edge on. Andere Winkel siehe hier. Hintergrund: ESO/Brunier.   
Bild

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

Bild  ←
Akkretionsscheibe um ein geladenes und rotierendes SL mit a=0.95, ℧=0.3, Innenradius ri=isco, Außenradius ra=10, Blickwinkel=89°    
Bild

Bild  ←
Verzerrte Oberfläche eines geladenen und rotierenden Körpers mit Erdoberfläche. a=0.95, ℧=0.3, Ausdehnung=1.01r+, Blickwinkel=89°    
Bild

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

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

Bild  ←
Prograder Orbit eines negativ geladenen Testpartikels (q=-¼) um ein wie oben rotierendes und positiv geladenes schwarzes Loch  Bild

Bild  ←
Nichtäquatorialer und retrograder Photonenorbit um ein mit a=½ und ℧=½ geladenes schwarzes Loch, konstanter Boyer Lindquist Radius    Bild

Bild  ←
Polarer Photonenorbit (Lz=0) um eine mit a=½ rotierende und mit ℧=¾ geladene nackte Singularität, konstanter Boyer Lindquist Radius    Bild

Bild  ←
Polarer Orbit (Lz=0) eines positiv geladenen Testpartikels (q=⅓) um ein positiv geladenes und rotierendes schwarzes Loch (℧=a=0.7)   Bild

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

Bild  ←
Geodätischer Orbit um eine nackte Singularität mit den gleichen Spin- und Ladungsparametern wie ein Beispiel weiter oben    
Bild

Bild  ←
Nichtäquatorialer und retrograder Photonenorbit um eine mit a=0.9 rotierende und ℧=0.9 geladene nackte Singularität    
Bild

Bild  ←
Retrograder Photonenorbit um eine nackte Singularität (a=0.99, ℧=0.99). Lokale äquatoriale Inklination: -2.5rad=-143.23945°    
Bild

Bild  ←
Stationärer Photonenorbit (E=0) um eine Ringsingularität (a=½, ℧=1). Außer bei r=1, θ=90° ist v framedrag überall <c, daher keine Ergosphären.    Bild

Bild  ←
Äquatorialer retrograder Photonenorbit, Singularität bei r=0→R=√(r²+a²)=a=½. Ergoring (grün) bei r=1, Umkehrpunkte bei r=0.8 und r=1.3484    Bild

Bild  ←
Retrograder Photonenorbit der 3. Art innerhalb eines mit a=¾ rotierenden und ℧=½ geladenen schwarzes Lochs, konstanter BL Radius    Bild


Bild  ←
Orbit eines negativ geladenen Partikels im Inneren eines positiv geladenen Reissner Nordström SL (vergleiche Dokuchaev, Fig. 1)    Bild
Bild
Simon Tyran aka Симон Тыран @ minds || vk || wikipedia || stackexchange || wolframBild

Benutzeravatar
Yukterez
Administrator
Beiträge: 207
Registriert: Mi 21. Okt 2015, 02:16

Kerr Newman Metrik

Beitragvon Yukterez » Di 22. Aug 2017, 21:46

Bild

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

Bild

wobei

Bild

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

Bild

Effektive Masse:

Bild

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):

Bild

Ko- und kontravariante Metrik:

Bild

Elektromagnetisches Potential:

Bild

Kovarianter elektromagnetischer Tensor:

Bild

Kontravarianter Maxwelltensor:

Bild

Eigenzeitableitungen der Koordinaten:

Bild

Damit lauten die Bewegungsgleichungen:

Bild

Bild

Bild

Bild

Kanonischer Viererimpuls:

Bild

Insgesamte Zeitdilatation eines neutralen Testpartikels, gewonnen aus dem Linienelement:

Bild

Zusätzliche Zeitdilatation eines geladenen Testpartikels:

Bild

Insgesamte Zeitdilatation eines geladenen Testpartikels:

Bild

Zusammenhang der ersten Eigenzeitableitungen mit den kovarianten Impulskomponenten:

Bild

Bild

Zusammenhang der ersten Eigenzeitableitungen mit den lokalen Dreiergeschwindigkeitskomponenten:

Bild     Bild

Bild

mit dem kontrahierten transversalen elektromagnetischen Potential

Bild

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

Bild

und das latitudinale Potential

Bild

mit dem Parameter

Bild

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

Bild

Radiale Fluchtgeschwindigkeit für ein neutrales Teilchen:

Bild

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

Bild

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

Bild

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

Bild

3. Erhaltungsgröße: Carter-Konstante

Bild

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

Bild

Die azimutalen und latitudinalen Stoßparameter sind

Bild

Gravitative Zeitdilatation eines ZAMO, unendlich am Horizont:

Bild

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

Bild

Frame-Dragging Winkelgeschwindigkeit im System den Koordinatenbuchhalters:

Bild

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

Bild

mit dem Zusammenhang

Bild

Axialer und koaxialer Gyrationsradius:

Bild

Axialer und koaxialer Umfang:

Bild

Die Radien der äquatorialen Photonenkreisorbits sind implizit gegeben durch:

Bild

Der letzte stabile Orbit (ISCO) eines neutralen Partikels ergibt sich aus dem Polynom

Bild

Radialkoordinaten der Horizonte und Ergosphären:

Bild

Kartesische Projektion:

Bild

Bild
Bild
Simon Tyran aka Симон Тыран @ minds || vk || wikipedia || stackexchange || wolframBild

Benutzeravatar
Yukterez
Administrator
Beiträge: 207
Registriert: Mi 21. Okt 2015, 02:16

Kerr Newman, Code

Beitragvon Yukterez » Do 28. Sep 2017, 07:50

Bild

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 - 12.04.2019, Version 22 |||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]; ClearAll["Local`*"];
Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
 
wp=MachinePrecision;
mt1=Automatic;
mt2={"EquationSimplification"-> "Residual"};
mt3={"ImplicitRungeKutta", "DifferenceOrder"-> 20};
mt4={"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}};
mt5={"EventLocator", "Event"-> (r[τ]-1001/1000 rA)};
mta=mt1;
dgl=1;
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 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+2;                                                   (* 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 *)

vrj[τ_]:=R'[τ]/Sqrt[Δi[τ]] Sqrt[Σi[τ] (1+μ v[τ]^2)];
vθj[τ_]:=Θ'[τ] Sqrt[Σi[τ] (1+μ v[τ]^2)];
vφj[τ_]:=Evaluate[(-(((a^2 Cos[θ[τ]]^2+r[τ]^2) (a^2+℧^2-2 r[τ]+r[τ]^2) Sin[θ[τ]] Sqrt[1-
μ^2 v[τ]^2] (-φ'[τ]-(a q ℧ r[τ])/((a^2 Cos[θ[τ]]^2+r[τ]^2) (a^2+℧^2-2 r[τ]+r[τ]^2))+
(ε Csc[θ[τ]]^2 (a (-a^2-℧^2+2 r[τ]-r[τ]^2) Sin[θ[τ]]^2+a (a^2+
r[τ]^2) Sin[θ[τ]]^2))/((a^2 Cos[θ[τ]]^2+r[τ]^2) (a^2+℧^2-2 r[τ]+r[τ]^2))+(a q ℧ r[τ] (a^2+
℧^2-2 r[τ]+r[τ]^2-a^2 Sin[θ[τ]]^2))/((a^2 Cos[θ[τ]]^2+r[τ]^2)^2 (a^2+℧^2-2 r[τ]+
r[τ]^2) (1-μ^2 v[τ]^2))))/((a^2+℧^2-2 r[τ]+r[τ]^2-a^2 Sin[θ[τ]]^2) Sqrt[((a^2+r[τ]^2)^2-
a^2 (a^2+℧^2-2 r[τ]+r[τ]^2) Sin[θ[τ]]^2)/(a^2 Cos[θ[τ]]^2+r[τ]^2)]))) /. sol][[1]]
vtj[τ_]:=Sqrt[vrj[τ]^2+vθj[τ]^2+vφj[τ]^2];
vr[τ_]:=vrj[τ]/vtj[τ]*v[τ];
vθ[τ_]:=vθj[τ]/vtj[τ]*v[τ];
vφ[τ_]:=vφj[τ]/vtj[τ]*v[τ];
VΦ[τ_]:=Sqrt[v[τ]^2-vθ[τ]^2-vr[τ]^2];
Vφ[τ_]:=If[q==0, Vφ[τ], VΦ[τ]];
 
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)) j[v0]^2;
Lζ:=vφ0 я/j[ν]+0((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+μ);
                                                                                  (* ISCO *)
isco = rISCO/.Solve[0 == rISCO (6 rISCO-rISCO^2-9 ℧^2+3 a^2)+4 ℧^2 (℧^2-a^2)-8 a (rISCO-℧^2)^(3/2) && rISCO>=rA, rISCO][[1]];
μ=If[Abs[v0]==1, 0, If[Abs[v0]<1, -1, 1]];                   (* Baryon: μ=-1, Photon: μ=0 *)

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 3) FLUCHTGESCHWINDIGKEIT UND BENÖTIGTER 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))]];
                                                  (* horizontaler Photonenkreiswinkel, i0 *)
iP[r0_, a_]:=Normal[iPh/.NSolve[1/(8 (r0^2+a^2 Cos[θ0]^2)^3) (a^2+(-2+r0) r0+
℧^2) (8 r0 (r0^2+a^2 Cos[θ0]^2) Sin[iPh]^2+1/((a^2-2 r0+r0^2+℧^2) (r0^2+
a^2 Cos[θ0]^2)) 8 a (Cos[iPh] Sin[θ0] (a^2-2 r0+r0^2+℧^2-a^2 Sin[θ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)]+(a (a^2+r0^2) Sin[θ0]^2+
a (-a^2+2 r0-r0^2-℧^2) Sin[θ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)]+(a (2 r0-℧^2) Cos[iPh] 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)])/((a^2+r0^2)^2-a^2 (a^2-
2 r0+r0^2+℧^2) Sin[θ0]^2))) (-(1/((a^2-2 r0+r0^2+℧^2) (r0^2+
a^2 Cos[θ0]^2)))2 a^2 Cot[θ0]^2 (Cos[iPh] Sin[θ0] (-a (a^2+r0^2) Sin[θ0]^2+a (a^2-2 r0+
r0^2+℧^2) Sin[θ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)]+((a^2+r0^2)^2 Sin[θ0]^2+a^2 (-a^2+2 r0-r0^2-
℧^2) Sin[θ0]^4) (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)]+(a (2 r0-℧^2) Cos[iPh] 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)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+
℧^2) Sin[θ0]^2)))+1/((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2)) 2 r0 (r0-
℧^2) Csc[θ0]^2 (Cos[iPh] Sin[θ0] (-a (a^2+r0^2) Sin[θ0]^2+a (a^2-2 r0+r0^2+
℧^2) Sin[θ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)]+((a^2+r0^2)^2 Sin[θ0]^2+a^2 (-a^2+2 r0-r0^2-℧^2) Sin[θ0]^4) (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)]+
(a (2 r0-℧^2) Cos[iPh] 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)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))))+1/((a^2-2 r0+r0^2+
℧^2) (r0^2+a^2 Cos[θ0]^2)) 8 Csc[θ0]^2 (Cos[iPh] Sin[θ0] (-a (a^2+r0^2) Sin[θ0]^2+a (a^2-
2 r0+r0^2+℧^2) Sin[θ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)]+((a^2+r0^2)^2 Sin[θ0]^2+a^2 (-a^2+2 r0-r0^2-℧^2) Sin[θ0]^4) (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)]+
(a (2 r0-℧^2) Cos[iPh] 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)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))) (1/((a^2-2 r0+r0^2+
℧^2) (r0^2+a^2 Cos[θ0]^2)) a^2 Cot[θ0]^2 (Cos[iPh] Sin[θ0] (-a (a^2+r0^2) Sin[θ0]^2+
a (a^2-2 r0+r0^2+℧^2) Sin[θ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)]+((a^2+r0^2)^2 Sin[θ0]^2+a^2 (-a^2+2 r0-r0^2-
℧^2) Sin[θ0]^4) (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)]+(a (2 r0-℧^2) Cos[iPh] 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)])/((a^2+r0^2)^2-a^2 (a^2-
2 r0+r0^2+℧^2) Sin[θ0]^2)))+1/((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2)) r0 (-r0+
℧^2) Csc[θ0]^2 (Cos[iPh] Sin[θ0] (-a (a^2+r0^2) Sin[θ0]^2+a (a^2-2 r0+r0^2+
℧^2) Sin[θ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)]+
((a^2+r0^2)^2 Sin[θ0]^2+a^2 (-a^2+2 r0-r0^2-℧^2) Sin[θ0]^4) (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)]+(a (2 r0-
℧^2) Cos[iPh] 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)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))))+1/((a^2-2 r0+r0^2+
℧^2)^2 (r0^2+a^2 Cos[θ0]^2)^2) Csc[θ0]^2 (Cos[iPh] Sin[θ0] (a^2-2 r0+r0^2+℧^2-
a^2 Sin[θ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)]+
(a (a^2+r0^2) Sin[θ0]^2+a (-a^2+2 r0-r0^2-℧^2) Sin[θ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)]+(a (2 r0-
℧^2) Cos[iPh] 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)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)))^2 (r0 (a^2 (3 a^2+
4 ℧^2+4 (a-℧) (a+℧) Cos[2 θ0]+a^2 Cos[4 θ0])+8 r0 (r0^3+2 a^2 r0 Cos[θ0]^2-
a^2 Sin[θ0]^2))+2 a^4 Sin[2 θ0]^2))==0,iPh,Reals]][[1]]/.C[1]->0
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 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-μ^2 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_]:=ε(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;
pτ=ε(r0^2+a^2)+℧ q r0-a Lz;

Vr[r_]:=Pr[r]^2-Δr[r](μ^2 r^2+(Lz-a ε)^2+Q);             (* effektives radiales Potential *)
Vτ=Pτ^2-Δ(μ^2 r[τ]^2+(Lz-a ε)^2+Q);
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));
                                       
DG1={
 
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
 
};

DG2={
 
t'[τ]==1/(Δ Σ Sin[θ[τ]]^2) (Lz (Δ XJ-a Sin[θ[τ]]^2 (r[τ]^2+a^2))+ε (-Δ XJ^2+
Sin[θ[τ]]^2 (r[τ]^2+a^2)^2)-q ℧ r[τ] Sin[θ[τ]]^2 (r[τ]^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/(Δ Σ Sin[θ[τ]]^2) (ε (-Δ XJ+a Sin[θ[τ]]^2 (r[τ]^2+a^2))+Lz (Δ-a^2 Sin[θ[τ]]^2)-
q ℧ r[τ] a Sin[θ[τ]]^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
 
};

DG3={
 
t'[τ]==1/(Δ Σ Sin[θ[τ]]^2) (Lz (Δ XJ-a Sin[θ[τ]]^2 (r[τ]^2+a^2))+ε (-Δ XJ^2+
Sin[θ[τ]]^2 (r[τ]^2+a^2)^2)-q ℧ r[τ] Sin[θ[τ]]^2 (r[τ]^2+a^2)),
t[0]==0,
 
r'[τ]==Sign[vr0] Sqrt[Vτ]/Σ,
r[0]==r0,

θ'[τ]==Sign[vθ0] Sqrt[Q-Cos[θ[τ]]^2 (a^2 (μ^2-ε^2)+Lz^2/Sin[θ[τ]]^2)]/Σ,
θ[0]==θ0,

φ'[τ]==1/(Δ Σ Sin[θ[τ]]^2) (ε (-Δ XJ+a Sin[θ[τ]]^2 (r[τ]^2+a^2))+Lz (Δ-a^2 Sin[θ[τ]]^2)-
q ℧ r[τ] a Sin[θ[τ]]^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
 
};

DGL = If[dgl==1, DG1, If[dgl==2, DG2, DG3]];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 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=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, FontFamily->"Consolas", FontSize->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"], " = ", If[μ==0, s[n0[ς[T]]], 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[N[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.011], Red, Point[
Xyz[xyZ[{x[T], y[T], z[T]}, w1], w2]]}},
ImageSize-> imgsize,
PlotRange-> {
{-(9 Sign[Abs[xx]]+1) PR, +(9 Sign[Abs[xx]]+1) PR},
{-(9 Sign[Abs[yy]]+1) PR, +(9 Sign[Abs[yy]]+1) PR},
{-(9 Sign[Abs[zz]]+1) PR, +(9 Sign[Abs[zz]]+1) PR}
},
SphericalRegion->False,
ImagePadding-> 1],

horizons[A, None, w1, w2],

If[a==0, {}, ParametricPlot3D[
Xyz[xyZ[{
Sin[prm] a,
Cos[prm] a,
0}, w1], w2],
{prm, 0, 2π},
PlotStyle -> {Thickness[0.005], Orange}]],

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-1/2 π/ω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/1}]]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 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"], " = ", If[μ==0, s[n0[ς[tp]]], 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[N[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.011], Red, Point[
Xyz[xyZ[{x[tp], y[tp], z[tp]}, w1], w2]]}},
ImageSize-> imgsize,
PlotRange-> {
{-(9 Sign[Abs[xx]]+1) PR, +(9 Sign[Abs[xx]]+1) PR},
{-(9 Sign[Abs[yy]]+1) PR, +(9 Sign[Abs[yy]]+1) PR},
{-(9 Sign[Abs[zz]]+1) PR, +(9 Sign[Abs[zz]]+1) PR}
},
SphericalRegion->False,
ImagePadding-> 1],

horizons[A, None, w1, w2],

If[a==0, {}, ParametricPlot3D[
Xyz[xyZ[{
Sin[prm] a,
Cos[prm] a,
0}, w1], w2],
{prm, 0, 2π},
PlotStyle -> {Thickness[0.005], Orange}]],

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]-1/2 π/ω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/1}]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 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 |||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)










Interpolationsfunktion für den schnelleren Output vieler Einzelbilder bei dafür etwas längerer Anlaufzeit:

Code: Alles auswählen

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 11) PLOT NACH KOORDINATENZEIT ||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

Plp=Round[Max[4, 2tk]];

rint=Evaluate[Interpolation[Table[{tint,R[д[tint]]},{tint,0,TMax,1/2}]]];
uint=Evaluate[Interpolation[Table[{tint,Θ[д[tint]]},{tint,0,TMax,1/2}]]];
pint=Evaluate[Interpolation[Table[{tint,Φ[д[tint]]},{tint,0,TMax,1/2}]]];
zint[τ_]:=rint[τ] Cos[uint[τ]];
xint[τ_]:=Sqrt[rint[τ]^2+a^2] Sin[uint[τ]] Cos[pint[τ]];
yint[τ_]:=Sqrt[rint[τ]^2+a^2] Sin[uint[τ]] Sin[pint[τ]];

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"], " = ", If[μ==0, s[n0[ς[T]]], s[n0[1/Sqrt[1-v[T]^2]]]], s["dt/dτ"], s[dp]},
{s[" R carts"], " = ", s[n0[Sqrt[xint[tk]^2+yint[tk]^2+zint[tk]^2]]], s["GM/c²"], s[dp]},
{s[" x carts"], " = ", s[n0[xint[tk]]], s["GM/c²"], s[dp]},
{s[" y carts"], " = ", s[n0[yint[tk]]], s["GM/c²"], s[dp]},
{s[" z carts"], " = ", s[n0[zint[tk]]], s["GM/c²"], s[dp]},
{s[" s dstnc"], " = ", s[n0[dst[T]]], s["GM/c²"], s[dp]},

{s[" r coord"], " = ", s[n0[rint[tk]]], s["GM/c²"], s[dp]},
{s[" φ longd"], " = ", s[n0[pint[tk] 180/π]], s["deg"], s[dp]},
{s[" θ lattd"], " = ", s[n0[uint[tk] 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[N[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}];

hz=horizons[A, None, 0, 0];
 
plot1a[{xx_, yy_, zz_, tk_, w1_, w2_}]:=                                     (* Animation *)
Show[

Graphics3D[{
{PointSize[0.011], Red, Point[
Xyz[xyZ[{xint[tk], yint[tk], zint[tk]}, w1], w2]]}},
ImageSize-> imgsize,
PlotRange-> {
{-(9 Sign[Abs[xx]]+1) PR, +(9 Sign[Abs[xx]]+1) PR},
{-(9 Sign[Abs[yy]]+1) PR, +(9 Sign[Abs[yy]]+1) PR},
{-(9 Sign[Abs[zz]]+1) PR, +(9 Sign[Abs[zz]]+1) PR}
},
SphericalRegion->False,
ImagePadding-> 1],

horizons[A, None, w1, w2],

If[a==0, {}, ParametricPlot3D[
Xyz[xyZ[{
Sin[prm] a,
Cos[prm] a,
0}, w1], w2],
{prm, 0, 2π},
PlotStyle -> {Thickness[0.005], Orange}]],

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[{xint[tt], yint[tt], zint[tt]}, w1], w2], {tt, If[TMax<0, Min[0, tk+d1], Max[0, tk-d1]], tk},
PlotStyle-> {Thickness[0.004]},
ColorFunction-> Function[{x, y, z, t},
Hue[0, 1, 0.5, If[TMax<0, Max[Min[(+tk+(-t+d1))/d1, 1], 0], Max[Min[(-tk+(t+d1))/d1, 1], 0]]]],
ColorFunctionScaling-> False,
PlotPoints-> Automatic,
MaxRecursion-> mrec]]],

If[tk==0, {},
Block[{$RecursionLimit = Mrec},
ParametricPlot3D[
Xyz[xyZ[{xint[tt], yint[tt], zint[tt]}, w1], w2],
{tt, 0, If[Tmax<0, Min[-1*^-16, tk+d1/3], Max[1*^-16, tk-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}]]










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

N@Q

Q0=With[{τ=0},Evaluate[((θ'[τ] (r[τ]^2+a^2 Cos[θ[τ]]^2))^2+(((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))^2 Csc[θ[τ]]^2-a^2 ((((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))^2+μ)) Cos[θ[τ]]^2)/.sol][[1]]]

Q1=With[{τ=tMax 99/100},Evaluate[((θ'[τ] (r[τ]^2+a^2 Cos[θ[τ]]^2))^2+(((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))^2 Csc[θ[τ]]^2-a^2 ((((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))^2+μ)) Cos[θ[τ]]^2)/.sol][[1]]]

Plot[{
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]],
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]],
Evaluate[((θ'[τ] (r[τ]^2+a^2 Cos[θ[τ]]^2))^2+(((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))^2 Csc[θ[τ]]^2-a^2 ((((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))^2+μ)) Cos[θ[τ]]^2)/.sol][[1]]
}, {τ, 0, tMax 99/100}, PlotPoints->100, ImageSize->300, Frame->True]










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

Code: Alles auswählen

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE 1: Erste Eigenzeitableitungen nach lokalen Geschwindigkeiten ||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]; ClearAll["Local`*"];
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = Sqrt[7^2-a^2];
θ0 = π/2;
φ0 = 0;
a  = 9/10;
℧  = 2/5;
q  = -1/2;
μ  =-1;
 
(* || Startwerte für die lokalen Geschwindigkeitskomponenten eingeben ||||||||||||||||||||||||||| *)
 
vr0 = 0;
vθ0 = 2/5 Sin[(2 π)/9];
vφ0 = 2/5 Cos[(2 π)/9];
v0  = Sqrt[vr0^2+vθ0^2+vφ0^2];
 
(* || Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ε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 μ^2]+(a (2 r0-℧^2) ((a q r0 (1-v0^2 μ^2) ℧ 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 μ^2]))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2);
 
L0 = (a q r0 (1-v0^2 μ^2) ℧ 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 μ^2];
 
Q0 = (vθ0^2 (r0^2+a^2 Cos[θ0]^2))/(1-v0^2 μ^2)+Cos[θ0]^2 (Csc[θ0]^2 ((a q r0 (1-v0^2 μ^2) ℧ 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 μ^2])^2-a^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 μ^2]+(a (2 r0-℧^2) ((a q r0 (1-v0^2 μ^2) ℧ 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 μ^2]))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))^2));
 
Ξ=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 1"
FindInstance[dT==dt && dR==dr && dΘ==dθ && dΦ==dφ, {dt,dr,dθ,dφ}, Reals]
N[%]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax |||| kerr.newman.yukterez.net |||| Simon Tyran, Vienna  ||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                                                   
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE 2: Erhaltungsgrößen ε, Lz, Q nach lokalen Geschwindigkeiten |||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]; ClearAll["Local`*"];
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = Sqrt[7^2-a^2];
θ0 = π/2;
φ0 = 0;
a  = 9/10;
℧  = 2/5;
q  = -1/2;
μ  =-1;
 
(* || Startwerte für die lokalen Geschwindigkeitskomponenten eingeben ||||||||||||||||||||||||||| *)
 
vr0 = 0;
vθ0 = 2/5 Sin[(2 π)/9];
vφ0 = 2/5 Cos[(2 π)/9];
v0  = Sqrt[vr0^2+vθ0^2+vφ0^2];
 
(* || Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ε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 μ^2]+(a (2 r0-℧^2) ((a q r0 (1-v0^2 μ^2) ℧ 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 μ^2]))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2);
 
L0 = (a q r0 (1-v0^2 μ^2) ℧ 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 μ^2];
 
Q0 = (vθ0^2 (r0^2+a^2 Cos[θ0]^2))/(1-v0^2 μ^2)+Cos[θ0]^2 (Csc[θ0]^2 ((a q r0 (1-v0^2 μ^2) ℧ 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 μ^2])^2-a^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 μ^2]+(a (2 r0-℧^2) ((a q r0 (1-v0^2 μ^2) ℧ 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 μ^2]))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))^2));
 
Ξ=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"
FindInstance[ε==ε0 && Lz==L0 && Q==Q0, {ε,Lz,Q}, Reals]
N[%]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax |||| kerr.newman.yukterez.net |||| Simon Tyran, Vienna  ||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE 3: Lokale Geschwindigkeit nach ersten Eigenzeitableitungen  |||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]; ClearAll["Local`*"];
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = Sqrt[7^2-a^2];
θ0 = π/2;
φ0 = 0;
a  = 9/10;
℧  = 2/5;
q  = -1/2;
μ  =-1;
 
(* || Startwerte für die ersten Eigenzeitableitungen eingeben ||||||||||||||||||||||||||| *)
 
dt = -((1701 (-24095+4 Sqrt[4819])+47026091025 Sqrt[(21 (1229-5 Sqrt[4819]))/(5902951+405 Sqrt[4819])]+142231604345 Sqrt[(101199 (1229-5 Sqrt[4819]))/(5902951+405 Sqrt[4819])])/(487677981 (-1229+5 Sqrt[4819])));
dr = 0;
dθ = (20 Sin[(2 π)/9])/Sqrt[101199];
dφ = (10064917571310-509342021892 Sqrt[4819]+3576385309875 Sqrt[(101199 (1229-5 Sqrt[4819]))/(5902951+405 Sqrt[4819])]-257016180174650625 Sqrt[(3687-15 Sqrt[4819])/(41320657+2835 Sqrt[4819])]+13988810657375 Sqrt[1/21 (5902951+405 Sqrt[4819])] Cos[(2 π)/9]-713519331725 Sqrt[4819/21 (5902951+405 Sqrt[4819])] Cos[(2 π)/9])/(116113805 (-3622484152+14508505 Sqrt[4819]));
 
(* || Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ε0 = ((q r0 ℧)/(r0^2+a^2 Cos[θ0]^2))+dt (1-(2 r0-℧^2)/(r0^2+a^2 Cos[θ0]^2))+(a dφ (2 r0-℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2);
L0 = (q a r0 ℧ Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)-(a dt  (2 r0-℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)+(dφ Sin[θ0]^2 ((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))/(r0^2+a^2 Cos[θ0]^2);
Q0 = ((dθ (r0^2+a^2 Cos[θ0]^2))^2+(((q a r0 ℧ Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)-(a dt  (2 r0-℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)+(dφ Sin[θ0]^2 ((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))/(r0^2+a^2 Cos[θ0]^2))^2 Csc[θ0]^2-a^2 ((((q r0 ℧)/(r0^2+a^2 Cos[θ0]^2))+dt (1-(2 r0-℧^2)/(r0^2+a^2 Cos[θ0]^2))+(a dφ (2 r0-℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2))^2+μ)) Cos[θ0]^2);
 
Ξ=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);

v0j = Abs[(Sqrt[δ Ξ^3 χ-ε0^2 Ξ^2 χ^2-2 a L0 ε0 Ξ^2 χ ℧^2-a^2 L0^2 Ξ^2 ℧^4+4 a L0 ε0 Ξ^2 χ r0+2 q ε0 Ξ χ^2 ℧ r0+4 a^2 L0^2 Ξ^2 ℧^2 r0+2 a L0 q Ξ χ ℧^3 r0-4 a^2 L0^2 Ξ^2 r0^2-4 a L0 q Ξ χ ℧ r0^2-q^2 χ^2 ℧^2 r0^2])/(ε0 Ξ χ+a L0 Ξ ℧^2-2 a L0 Ξ r0-q χ ℧ r0)];
vrj = dr/Sqrt[δ] Sqrt[Ξ (1+μ v0j^2)];
vθj = dθ Sqrt[Ξ (1+μ v0j^2)];
vφj = -(((a^2 Cos[θ0]^2+r0^2) (a^2+℧^2-2 r0+r0^2) Sin[θ0] Sqrt[1-μ^2 v0j^2] (-dφ-(a q ℧ r0)/((a^2 Cos[θ0]^2+r0^2) (a^2+℧^2-2 r0+r0^2))+(ε0 Csc[θ0]^2 (a (-a^2-℧^2+2 r0-r0^2) Sin[θ0]^2+a (a^2+r0^2) Sin[θ0]^2))/((a^2 Cos[θ0]^2+r0^2) (a^2+℧^2-2 r0+r0^2))+(a q ℧ r0 (a^2+℧^2-2 r0+r0^2-a^2 Sin[θ0]^2))/((a^2 Cos[θ0]^2+r0^2)^2 (a^2+℧^2-2 r0+r0^2) (1-μ^2 v0j^2))))/((a^2+℧^2-2 r0+r0^2-a^2 Sin[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+℧^2-2 r0+r0^2) Sin[θ0]^2)/(a^2 Cos[θ0]^2+r0^2)]));

(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code 3"
Simplify[Solve[v0==v0j && vr0==vrj && vθ0==vθj && vφ0==vφj, {v0,vr0,vφ0,vθ0}, Reals]]
N[%]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax |||| kerr.newman.yukterez.net |||| Simon Tyran, Vienna  ||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                               
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE 4: Lokale Geschwindigkeit nach Erhaltungsgrößen ε, Lz, Q ||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]; ClearAll["Local`*"];
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = Sqrt[7^2-a^2];
θ0 = π/2;
φ0 = 0;
a  = 9/10;
℧  = 2/5;
q  = -1/2;
μ  =-1;
 
(* || Erhaltungsgrößen Gesamtenergie, axialer Drehimpuls & Carter Konstante eingeben |||| *)
 
ε = 0.90688763;
Lz = 2.3240259;
Q = 3.7925614;
 
(* || 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];
P=ε(r0^2+a^2)+℧ q r0-a Lz;
Vr=P^2-δ(μ^2 r0^2+(Lz-a ε)^2+Q);
Vθ=Q-Cos[θ0]^2(a^2(μ^2-ε^2)+Lz^2/Sin[θ0]^2);

dT=Abs[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))];
dR=Sqrt[Abs[Vr]]/Ξ;
dΘ=Sqrt[Abs[Vθ]]/Ξ;
dΦ=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);

v0j=Abs[(Sqrt[δ Ξ^3 χ-ε^2 Ξ^2 χ^2-2 a Lz ε Ξ^2 χ ℧^2-a^2 Lz^2 Ξ^2 ℧^4+4 a Lz ε Ξ^2 χ r0+2 q ε Ξ χ^2 ℧ r0+4 a^2 Lz^2 Ξ^2 ℧^2 r0+2 a Lz q Ξ χ ℧^3 r0-4 a^2 Lz^2 Ξ^2 r0^2-4 a Lz q Ξ χ ℧ r0^2-q^2 χ^2 ℧^2 r0^2])/(ε Ξ χ+a Lz Ξ ℧^2-2 a Lz Ξ r0-q χ ℧ r0)];
vrj=Abs[dR Sqrt[Ξ (1+μ v0j^2)]/Sqrt[δ]];
vθj=Abs[dΘ Sqrt[Ξ (1+μ v0j^2)]];
vφj=-(((a^2 Cos[θ0]^2+r0^2) (a^2+℧^2-2 r0+r0^2) Sin[θ0] Sqrt[1-μ^2 v0j^2] (-dΦ-(a q ℧ r0)/((a^2 Cos[θ0]^2+r0^2) (a^2+℧^2-2 r0+r0^2))+(ε Csc[θ0]^2 (a (-a^2-℧^2+2 r0-r0^2) Sin[θ0]^2+a (a^2+r0^2) Sin[θ0]^2))/((a^2 Cos[θ0]^2+r0^2) (a^2+℧^2-2 r0+r0^2))+(a q ℧ r0 (a^2+℧^2-2 r0+r0^2-a^2 Sin[θ0]^2))/((a^2 Cos[θ0]^2+r0^2)^2 (a^2+℧^2-2 r0+r0^2) (1-μ^2 v0j^2))))/((a^2+℧^2-2 r0+r0^2-a^2 Sin[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+℧^2-2 r0+r0^2) Sin[θ0]^2)/(a^2 Cos[θ0]^2+r0^2)]));
vtj=Sqrt[vrj^2+vθj^2+vφj^2];
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code 4"
FullSimplify[FindInstance[{v0==Re@v0j && vθ0==Re@vθj && vφ0==Re@vφj && vr0==Re@Sqrt[v0^2-vφ0^2-vθ0^2]}, {v0, vr0, vθ0, vφ0}]]
N[%]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax |||| kerr.newman.yukterez.net |||| Simon Tyran, Vienna  ||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                               
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE 5: Erste Eigenzeitableitungen nach Erhaltungsgrößen ε, Lz, Q ||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]; ClearAll["Local`*"];
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = Sqrt[7^2-a^2];
θ0 = π/2;
φ0 = 0;
a  = 9/10;
℧  = 2/5;
q  = -1/2;
μ  =-1;
 
(* || Erhaltungsgrößen Gesamtenergie, axialer Drehimpuls & Carter Konstante eingeben |||| *)
 
ε = 0.90688763;
Lz = 2.3240259;
Q = 3.7925614;
 
(* || 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;
щ=(q ℧ r0 (a^2+r0^2))/(δ Ξ);

gtt=1-(2 r0-℧^2)/Ξ;
grr=-Ξ/δ;
gθθ=-Ξ;
gφφ=-χ/Ξ Sin[θ0]^2;
gtφ=a (2r0-℧^2) Sin[θ0]^2/Ξ;

P=ε(r0^2+a^2)+℧ q r0-a Lz;
Vr=P^2-δ(μ^2 r0^2+(Lz-a ε)^2+Q);
Vθ=Q-Cos[θ0]^2(a^2(μ^2-ε^2)+Lz^2/Sin[θ0]^2);

dT=Abs[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))];
dR=Sqrt[Abs[Vr]]/Ξ;
dΘ=Sqrt[Abs[Vθ]]/Ξ;
dΦ=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);

(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code 5"
FullSimplify[FindInstance[{dt==dT && dθ==dΘ && dφ==dΦ && gtt dt^2+grr dr^2+gθθ dθ^2+gφφ dφ^2+gtφ dt dφ==-μ}, {dt, dr, dθ, dφ}]]
N[%]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax |||| kerr.newman.yukterez.net |||| Simon Tyran, Vienna  ||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)









◬ Zur Überprüfung der numerischen Stabilität wird nach dem Plotten der Bahn der Erhaltungsgrößen Testmonitor (Codefenster 4) ausgeführt. dg1: alle Bewegungsgleichung über die zweite Zeitableitung (empfohlene Integrationsmethode mt0 oder mt3), dg2: r' und θ' über zweite, t' und φ' über die erste Zeitableitung (empfohlene Integrationsmethode mt2), dg3: alle Bewegungsgleichgungen über die erste Zeitableitung (diese Methode funktioniert wegen dem ± vor r' und θ' nur bis zu den radialen und poloidialen Umkehrpunkten und wird außer zu Testzwecken nicht empfohlen).

▣ Raytracing Code für schwarze Löcher, stellare Objekte und Akkretionsscheiben: raytracing.yukterez.net
Bild
Simon Tyran aka Симон Тыран @ minds || vk || wikipedia || stackexchange || wolframBild

Benutzeravatar
Yukterez
Administrator
Beiträge: 207
Registriert: Mi 21. Okt 2015, 02:16

Kerr Newman Metrik

Beitragvon Yukterez » Di 27. Mär 2018, 21:46

Bild
Simon Tyran aka Симон Тыран @ minds || vk || wikipedia || stackexchange || wolframBild


Zurück zu „Yukterez Notizblock“

Wer ist online?

Mitglieder in diesem Forum: 0 Mitglieder und 0 Gäste