Kerr Newman Metric

[Yukterez]

Kerr-Newman, second order differential equations of motion for a charged particle and photons. Animations by Simon Tyran, Vienna (Yukterez)

This is the english version.   Deutschsprachige Version auf kerr.newman.yukterez.net und Yukipedia.

  ↗
Animation Parameter: Bookkeeper time (t) and viewpoint (xy, zy). Charge of the BH: ℧=0.9, spin: a=0.4, charge of the testparticle: q=1, v0=vθ0=c/2  

  ↗
Animation parameter: Spinparameter a and elektric charge ℧. tmax=min(500, t(τplunge-10⁻⁶))

  ↗
Animation parameter: Bookkeeper coordinate time (t) and viewpoint angle (xy, zy). ●: FREFO, ●: ZAMO

  ↗
Shadow of an extremal Kerr-Newman BH with a²+℧²=M², Angle of view: edge on. For other parameters see here. Raw material: Commons.

Line element in Boyer Lindquist coordinates, metric signature (+,-,-,-):



Shorthand terms:



with the spin parameter â=Jc/G/M or in dimensionless units a=â/M, the specific electric charge Ω=⚼·√(K/G) and the dimensionless charge ℧=Ω/M. Here we use the units G=M=c=K=1 with lengths in GM/c² and times in GM/c³. The relation between the mass-equivalent of the total energy and the irreducible mass Mirr is



For testparticles with mass μ=-1, for photons μ=0. The specific charge of the test particle is q. Transformation rule for co- and contravariant indices (superscripted letters are not powers but indices):



Co- and contravariant metric:

Elektromagnetic potential:



Covariant elektromagnetic tensor:



Contravarianter Maxwell-tensor:



Einstein-Tensor Ĝ:



Covariant components:











Contravariant components:

1. Constant of motion: Carter's constant

2. Constant of motion: Total energy

3. Constant of motion: axial angular momentum



The relation between the first proper time derivatives and the local velocity components is







and the total 3-velocity



for testparticles, and v=1 for photons. ω is the Frame-Dragging angular velocity observed at infinity:



Gravitative time dilatation, Lapse function:



Escape velocity for a neutral particle:



First proper time derivatives of the coordinates:



The equations of motion for a neutral particle (with μ=-1, q=0) and a photon (with μ=0, q=0) are therefore:











and those for an electrically charged test particle (with μ=-1, q≠0):











First proper time derivatives of the space coordinates, initial conditions in terms of the local 3-velocity:





Radial coordinates of the horizons and ergospheres:

[Yukterez]

Input: Line element and Maxwell tensor, Output: Einstein tensor, Christoffel symbols and equations of motion

Code: Alles auswählen

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

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 Δ;
                    (* Koordinaten, Dimensionen, magnetisches Monopol, elektrische Ladung *)
x={t, r, θ, φ}; n=4; P=0; Ω=℧;
                                                                         "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, 4}, {k, 1, 4}], Reals];
Subscript["F", μσ] -> MatrixForm[F]

f=FullSimplify[Table[Sum[
inversemetrik[[i, k]] inversemetrik[[j, l]] F[[k, l]],
{k, 1, 4}, {l, 1, 4}], {i, 1, 4}, {j, 1, 4}], Reals];
"F"^μσ -> MatrixForm[f]
                                                                           "Einstein Tensor"
G=FullSimplify[Table[Sum[
-2(inversemetrik[[k,l]] F[[i, k]] F[[j, l]] - metrik[[i, j]] F[[k, l]] f[[k, l]]),
{k, 1, 4}, {l, 1, 4}], {i, 1, 4}, {j, 1, 4}], Reals];
Subscript["G", μσ] -> MatrixForm[G]

H=FullSimplify[Table[Sum[
inversemetrik[[i, k]] inversemetrik[[j, l]] G[[k, l]], {k, 1, 4}, {l, 1, 4}],
{i, 1, 4}, {j, 1, 4}], Reals];
"G"^μσ -> MatrixForm[H]
                                                                       "Christoffel Symbole"
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}]];
christoffelsymbole=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}];

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

list[christoffelsymbole[[1]][[2]][[1]]]
list[christoffelsymbole[[1]][[3]][[1]]]
list[christoffelsymbole[[1]][[4]][[2]]]
list[christoffelsymbole[[1]][[4]][[3]]]
list[christoffelsymbole[[2]][[1]][[1]]]
list[christoffelsymbole[[2]][[2]][[2]]]
list[christoffelsymbole[[2]][[3]][[2]]]
list[christoffelsymbole[[2]][[3]][[3]]]
list[christoffelsymbole[[2]][[4]][[1]]]
list[christoffelsymbole[[2]][[4]][[4]]]
list[christoffelsymbole[[3]][[1]][[1]]]
list[christoffelsymbole[[3]][[2]][[2]]]
list[christoffelsymbole[[3]][[3]][[2]]]
list[christoffelsymbole[[3]][[3]][[3]]]
list[christoffelsymbole[[3]][[4]][[1]]]
list[christoffelsymbole[[3]][[4]][[4]]]
list[christoffelsymbole[[4]][[2]][[1]]]
list[christoffelsymbole[[4]][[3]][[1]]]
list[christoffelsymbole[[4]][[4]][[2]]]
list[christoffelsymbole[[4]][[4]][[3]]]
                                                                      "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]]''[τ]==FullSimplify[rplc[geodäsie[[i]]], Reals]}, {i, 1, n}];

bewegungsgleichung[[1]][[1]]
bewegungsgleichung[[2]][[1]]
bewegungsgleichung[[3]][[1]]
bewegungsgleichung[[4]][[1]]

ClearAll[Σ, Δ, Χ];
                                                                            "Zeitdilatation"
t'[τ]==rple[Simplify[Normal[Solve[
-μ==gtt Т^2+grr r\.b4^2+gθθ θ\.b4^2+gφφ φ\.b4^2 + 2 gtφ Т φ\.b4, Т, Reals]]][[2]][[1]][[2]]]

Simulator code for photons, charged and neutral test particles

Code: Alles auswählen

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||| Mathematica Syntax | KERR-NEWMAN-SIMULATOR | yukterez.net | Version 24.03.2018 ||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]
 
mt1={"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}};
mt2={"EventLocator", "Event"-> (r[τ]-1001/1000 rA)};
mt3={"ImplicitRungeKutta", "DifferenceOrder"-> 20};
mt4={"EquationSimplification"-> "Residual"};
mt0=Automatic;
mta=mt4;
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];              (* Integrationsende *)
 
r0=7;                                                                      (* Startradius *)
θ0=π/2;                                                                    (* Breitengrad *)
φ0=0;                                                                       (* Längengrad *)
a=9/10;                                                                  (* Spinparameter *)
℧=2/5;                                          (* spezifische Ladung des schwarzen Lochs *)
q=-1/10;                                            (* spezifische Ladung des Testkörpers *)
 
v0=4/10;                                                        (* Anfangsgeschwindigkeit *)
α0=0;                                                        (* vertikaler Abschusswinkel *)
ψ0=ArcTan[5/6];                                                 (* Bahninklinationswinkel *)
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 2) GESCHWINDIGKEITS-, ENERGIE UND DREHIMPULSKOMPONENTEN ||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
vr0=v0 Sin[α0];                                     (* radiale Geschwindigkeitskomponente *)
vθ0=v0 Cos[α0] Sin[ψ0];                      (* longitudinale  Geschwindigkeitskomponente *)
vφ0=v0 Cos[α0] Cos[ψ0];                        (* latitudinale Geschwindigkeitskomponente *)
 
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];
 
ε=Sqrt[δ Ξ/χ]/j[v0]+Lz ω0; εζ:=Sqrt[Δ Σ/Χ]/j[ν]+Lz ωζ;   (* Energie und Impulskomponenten *)
Lz=vφ0 Ы/j[v0]; Lζ:=vφ0 я/j[ν];
pθ0=vθ0 Sqrt[Ξ]/j[v0]; pθζ:=θ'[τ] Σ;
pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];
Q=Limit[pθ0^2+(Lz^2 Csc[θ1]^2-a^2 (ε^2+μ)) Cos[θ1]^2, θ1->θ0];        (* Carter Konstante *)
Qζ:=pθζ^2+(Lz^2 Csc[θ[τ]]^2-a^2 (εζ^2+μ)) Cos[θ[τ]]^2;                (* Carter Parameter *)
k=Q+Lz^2+a^2 (ε^2+μ); kζ:=Qζ+Lz^2+a^2 (εζ^2+μ);                               (* Carter k *)

μ=If[v0==1, 0, -1];                                          (* Baryon: μ=-1, Photon: μ=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+μ 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;
Δ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;
 
т[τ_]:=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 *)
vr[τ_]:=R'[τ] Sqrt[Σi[τ]/Δi[τ]];
vθ[τ_]:=Θ'[τ] Sqrt[Σi[τ]];
vφ[τ_]:=Φ'[τ] Σi[τ]/яi[τ]-й[τ] ю[τ];
V[τ_]:=If[μ==0, 1, Re[Sqrt[-ς[τ]^2+ю[τ]^2]/ю[τ]]];        (* lokale Dreiergeschwindigkeit *)
vd[τ_]:=Abs[Sqrt[(Δ Σ-Χ(ε-Lz ωζ)^2)/(μ Χ (ε-Lz ωζ)^2)]];
v[τ_]:=If[μ==0, 1, Evaluate[vlt'[τ]/.sol][[1]]];
ν:=If[μ==0, 1, Re[Sqrt[(Δ Σ-Χ(εζ-Lζ ωζ)^2)/(μ Χ (εζ-Lζ ωζ)^2)]]];
dst[τ_]:=Evaluate[str[τ]/.sol][[1]];                                           (* Strecke *)
 
pΘ[τ_]:=Evaluate[Ξ θ'[τ] /. sol][[1]];
pR[τ_]:=Evaluate[r'[τ] Ξ/δ /. sol][[1]];
sh[τ_]:=Re[Sqrt[ß[τ]^2-Ω[τ]^2]];
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[N[Simplify[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''[τ]==(4 (((a^2+a^2 Cos[2 θ[τ]]+2 (℧^2-r[τ]) r[τ]) (a^2+r[τ]^2) r'[τ] t'[τ])/(a^2+℧^2-2 r[τ]+r[τ]^2)+a^2 (-℧^2+2 r[τ]) Sin[2 θ[τ]] t'[τ] θ'[τ]-(1/(a^2+℧^2-2 r[τ]+r[τ]^2))a (a^4+4 ℧^2 r[τ]^3-6 r[τ]^4-3 a^2 r[τ] (-℧^2+r[τ])+a^2 Cos[2 θ[τ]] (a^2+(℧^2-r[τ]) r[τ])) Sin[θ[τ]]^2 r'[τ] φ'[τ]-2 a^3 Cos[θ[τ]] (-℧^2+2 r[τ]) Sin[θ[τ]]^3 θ'[τ] φ'[τ]))/(a^2+a^2 Cos[2 θ[τ]]+2 r[τ]^2)^2,
 
t'[0]==t10,
t[0]==0,
 
r''[τ]==(1/(8 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3))(-((8 (a^2 Cos[θ[τ]]^2+(a^2+℧^2-a^2 Cos[θ[τ]]^2) r[τ]-r[τ]^2) (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 (r'[τ])^2)/(a^2+℧^2-2 r[τ]+r[τ]^2))+8 (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) (a^2+℧^2-2 r[τ]+r[τ]^2) (t'[τ])^2+16 a^2 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[θ[τ]] r'[τ] θ'[τ]+8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 (a^2+℧^2-2 r[τ]+r[τ]^2) (θ'[τ])^2-16 a (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) (a^2+℧^2-2 r[τ]+r[τ]^2) Sin[θ[τ]]^2 t'[τ] φ'[τ]+(a^2+℧^2-2 r[τ]+r[τ]^2) Sin[θ[τ]]^2 (a^2 (3 a^2+4 ℧^2+4 (a^2-℧^2) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]]) r[τ]+16 a^2 Cos[θ[τ]]^2 r[τ]^3+8 r[τ]^5-8 a^2 r[τ]^2 Sin[θ[τ]]^2+2 a^4 Sin[2 θ[τ]]^2) (φ'[τ])^2),
 
r'[0]==r10,
r[0]==r0,
 
θ''[τ]==(1/(16 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3))(-((16 a^2 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[θ[τ]] (r'[τ])^2)/(a^2+℧^2-2 r[τ]+r[τ]^2))-8 a^2 (℧^2-2 r[τ]) Sin[2 θ[τ]] (t'[τ])^2-32 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 r'[τ] θ'[τ]+16 a^2 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[θ[τ]] (θ'[τ])^2+16 a (℧^2-2 r[τ]) (a^2+r[τ]^2) Sin[2 θ[τ]] t'[τ] φ'[τ]+(a^4 (3 a^2-5 ℧^2+4 (a^2+℧^2) Cos[2 θ[τ]]+(a^2+℧^2) Cos[4 θ[τ]])+a^2 (11 a^2-8 ℧^2+4 (3 a^2+2 ℧^2) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]]) r[τ]^2+8 a^2 (2+Cos[2 θ[τ]]) r[τ]^4+8 r[τ]^6+8 a^4 (3+Cos[2 θ[τ]]) r[τ] Sin[θ[τ]]^2+32 a^2 r[τ]^3 Sin[θ[τ]]^2) Sin[2 θ[τ]] (φ'[τ])^2),
 
θ'[0]==Θ10,
θ[0]==θ0,
 
φ''[τ]==-(1/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2))(-((8 a (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) r'[τ] t'[τ])/(a^2+℧^2-2 r[τ]+r[τ]^2))+8 a Cot[θ[τ]] (℧^2-2 r[τ]) t'[τ] θ'[τ]+((a^2 (3 a^2+8 ℧^2+4 a^2 Cos[2 θ[τ]]+a^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) r'[τ] φ'[τ])/(a^2+℧^2-2 r[τ]+r[τ]^2)+Cot[θ[τ]] (a^2 (3 a^2-4 ℧^2+4 (a^2+℧^2) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]])+16 a^2 Cos[θ[τ]]^2 r[τ]^2+8 r[τ]^4+16 a^2 r[τ] Sin[θ[τ]]^2) θ'[τ] φ'[τ]),
 
φ'[0]==Φ10,
φ[0]==φ0,
 
str'[τ]==vd[τ]/Max[1*^-16, Abs[Sqrt[1-vd[τ]^2]]],
str[0]==0,
vlt'[τ]==vd[τ],
vlt[0]==0
 
};
                                         
DG2={
 
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'[τ]==vd[τ]/Max[1*^-16, Abs[Sqrt[1-vd[τ]^2]]],
str[0]==0,
vlt'[τ]==vd[τ],
vlt[0]==0
 
};
 
DGL=If[q==0, DG1, DG2];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 8) INTEGRATION |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
sol=NDSolve[DGL, {t, r, θ, φ, vlt, str}, {τ, 0, tmax},
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.2r0;                                                                   (* Plot Range *)
VP={r0, r0, r0};                                                     (* Perspektive x,y,z *)
d1=10;                                                                    (* Schweiflänge *)
plp=50;                                                            (* Flächenplot 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[" d 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 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[" 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[Q]], 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[" 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[" ω 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 obsvd"], " = ", s[n0[ß[T]]], s["c"], s[dp]},
{s[" v escpe"], " = ", s[n0[ж[T]]], s["c"], s[dp]},
{s[" v delay"], " = ", s[n0[sh[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]]],
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.003], Gray},
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]]],
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, TMax, TMax, TMax}]]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 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[" d 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 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[" 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[Q]], 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[" 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[" ω 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 obsvd"], " = ", s[n0[ß[tp]]], s["c"], s[dp]},
{s[" v escpe"], " = ", s[n0[ж[tp]]], s["c"], s[dp]},
{s[" v delay"], " = ", s[n0[sh[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, 0, If[tp<0, Min[-1*^-16, tp+d1/3], Max[1*^-16, tp-d1/3]]},
PlotStyle-> {Thickness[0.003], Gray},
PlotPoints-> Automatic,
MaxRecursion-> mrec]]],
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]]],
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, tMax, tMax, tMax}]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 13) EXPORTOPTIONEN |||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
(* Export als HTML Dokument *)
(* Export["dateiname.html", EvaluationNotebook[], "GraphicsOutput" -> "PNG"] *)
(* Export direkt als Bildsequenz *)
(* Do[Export["dateiname" <> ToString[tk] <> ".png", Rasterize[...]   ], {tk, 0, 10, 5}]   *)
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||| http://kerr.newman.yukerez.net ||||| Simon Tyran, Vienna |||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

Converter for initial conditions and constants of motion

Code: Alles auswählen

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE A: Lokale Geschwindigkeit nach Erhaltungsgrößen ε, Lz, Q ||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 7;
θ0 = π/2;
φ0 = 0;
a  = 9/10;
μ  =-1;
℧  = 0;
q  = 0;
 
(* || Erhaltungsgrößen Gesamtenergie, axialer Drehimpuls & Carter Konstante eingeben |||| *)
 
ε  = (72 Sqrt[3/2136769])/7+5 Sqrt[3581/105087];
Lz = (2 Sqrt[105087/61])/35;
Q  = 700/183;
 
(* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *)
 
ε0 = (Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)]+(a vφ0 (2 r0-℧^2) Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))/Sqrt[1-v0^2];
L0 = (vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^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 ((vφ0^2 (-(a^2+r0^2)^2+a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2))/((-1+v0^2) (r0^2+a^2 Cos[θ0]^2))-a^2 (-1+(Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)]+(a vφ0 (2 r0-℧^2) Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))^2/(1-v0^2)));
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code A"
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[%]

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || Syntax: Mathematica ||| http://kerr.newman.yukterez.net ||| Simon Tyran, Vienna  || *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                                                   
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE B: Lokale Geschwindigkeit nach ersten Eigenzeitableitungen  |||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 7;
θ0 = π/2;
φ0 = 0;
a  = 9/10;
μ  =-1;
℧  = 0;
q  = 0;
 
(* || Startwerte für die ersten Eigenzeitableitungen eingeben ||||||||||||||||||||||||||| *)
 
dt = (5 Sqrt[35029/10743])/7;
dr = 0;
dθ = 10/(7 Sqrt[1281]);
dφ = (300 Sqrt[3/125438849])/7+40 Sqrt[3/2136769];
 
(* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *)
 
ε0 = (Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)]+(a vφ0 (2 r0-℧^2) Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))/Sqrt[1-v0^2];
L0 = (vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^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 ((vφ0^2 (-(a^2+r0^2)^2+a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2))/((-1+v0^2) (r0^2+a^2 Cos[θ0]^2))-a^2 (-1+(Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)]+(a vφ0 (2 r0-℧^2) Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))^2/(1-v0^2)));
 
(* || Benötigte Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
Ξ=r0^2+a^2 Cos[θ0]^2;
δ=r0^2-2r0+a^2;
j[v_]:=Sqrt[1+μ v^2];
pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];
pθ0=vθ0 Sqrt[Ξ]/j[v0];
 
dT=(-a^2 (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) (q r0 Sqrt[1-v0^2] ℧+(a^2+r0^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-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2));
dR=-(vr0/((-1+v0^2) Sqrt[-((r0^2+a^2 Cos[θ0]^2)/((-1+v0^2) (a^2-2 r0+r0^2+℧^2)))]));
dΘ=vθ0/(Sqrt[1-v0^2] Sqrt[r0^2+a^2 Cos[θ0]^2]);
dΦ=(Csc[θ0] (a q r0 ℧ Sin[θ0]+(vφ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)])/Sqrt[1-v0^2]+(1/Sqrt[1-v0^2])a (2 r0-℧^2) Sin[θ0] (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 vφ0 (2 r0-℧^2) Sin[θ0])/((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)]))))/((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2));
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code B"
Reduce[dT==dt && dR==dr && dΘ==dθ && dΦ==dφ && v0>0, {v0,vr0,vφ0,vθ0}]
N[%]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || Syntax: Mathematica ||| http://kerr.newman.yukterez.net ||| Simon Tyran, Vienna  || *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                                                   
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE C: Erste Eigenzeitableitungen nach Erhaltungsgrößen ε, Lz, Q ||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 7;
θ0 = π/2;
φ0 = 0;
a  = 9/10;
μ  =-1;
℧  = 0;
q  = 0;
 
(* || Erhaltungsgrößen Gesamtenergie, axialer Drehimpuls & Carter Konstante eingeben |||| *)
 
ε  = (72 Sqrt[3/2136769])/7+5 Sqrt[3581/105087];
Lz = (2 Sqrt[105087/61])/35;
Q  = 700/183;
 
(* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *)
 
ε0 = (Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)]+(a vφ0 (2 r0-℧^2) Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))/Sqrt[1-v0^2];
L0 = (vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^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 ((vφ0^2 (-(a^2+r0^2)^2+a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2))/((-1+v0^2) (r0^2+a^2 Cos[θ0]^2))-a^2 (-1+(Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)]+(a vφ0 (2 r0-℧^2) Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))^2/(1-v0^2)));
 
(* || Benötigte Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
Ξ=r0^2+a^2 Cos[θ0]^2;
δ=r0^2-2r0+a^2;
j[v_]:=Sqrt[1+μ v^2];
pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];
pθ0=vθ0 Sqrt[Ξ]/j[v0];
 
dT=(-a^2 (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) (q r0 Sqrt[1-v0^2] ℧+(a^2+r0^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-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2));
dR=-(vr0/((-1+v0^2) Sqrt[-((r0^2+a^2 Cos[θ0]^2)/((-1+v0^2) (a^2-2 r0+r0^2+℧^2)))]));
dΘ=vθ0/(Sqrt[1-v0^2] Sqrt[r0^2+a^2 Cos[θ0]^2]);
dΦ=(Csc[θ0] (a q r0 ℧ Sin[θ0]+(vφ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)])/Sqrt[1-v0^2]+(1/Sqrt[1-v0^2])a (2 r0-℧^2) Sin[θ0] (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 vφ0 (2 r0-℧^2) Sin[θ0])/((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)]))))/((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2));
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code C"
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[%]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || Syntax: Mathematica ||| http://kerr.newman.yukterez.net ||| Simon Tyran, Vienna  || *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                                                   
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE D: Erste Eigenzeitableitungen nach lokalen Geschwindigkeiten ||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 7;
θ0 = π/2;
φ0 = 0;
a  = 9/10;
μ  =-1;
℧  = 0;
q  = 0;
 
(* || Startwerte für die ersten Eigenzeitableitungen eingeben ||||||||||||||||||||||||||| *)
 
vr0 = 0;
vθ0 = 2/Sqrt[61];
vφ0 = 12/(5 Sqrt[61]);
v0  = Sqrt[vr0^2+vθ0^2+vφ0^2];
 
(* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *)
 
ε0 = (Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)]+(a vφ0 (2 r0-℧^2) Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))/Sqrt[1-v0^2];
L0 = (vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^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 ((vφ0^2 (-(a^2+r0^2)^2+a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2))/((-1+v0^2) (r0^2+a^2 Cos[θ0]^2))-a^2 (-1+(Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)]+(a vφ0 (2 r0-℧^2) Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))^2/(1-v0^2)));
 
(* || Benötigte Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
Ξ=r0^2+a^2 Cos[θ0]^2;
δ=r0^2-2r0+a^2;
j[v_]:=Sqrt[1+μ v^2];
pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];
pθ0=vθ0 Sqrt[Ξ]/j[v0];
 
dT=(-a^2 (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) (q r0 Sqrt[1-v0^2] ℧+(a^2+r0^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-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2));
dR=-(vr0/((-1+v0^2) Sqrt[-((r0^2+a^2 Cos[θ0]^2)/((-1+v0^2) (a^2-2 r0+r0^2+℧^2)))]));
dΘ=vθ0/(Sqrt[1-v0^2] Sqrt[r0^2+a^2 Cos[θ0]^2]);
dΦ=(Csc[θ0] (a q r0 ℧ Sin[θ0]+(vφ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)])/Sqrt[1-v0^2]+(1/Sqrt[1-v0^2])a (2 r0-℧^2) Sin[θ0] (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 vφ0 (2 r0-℧^2) Sin[θ0])/((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)]))))/((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2));
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code D"
Reduce[dT==dt && dR==dr && dΘ==dθ && dΦ==dφ, {dt,dr,dθ,dφ}]
N[%]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || Syntax: Mathematica ||| http://kerr.newman.yukterez.net ||| Simon Tyran, Vienna  || *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                                                   
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE E: Erhaltungsgrößen ε, Lz, Q nach lokalen Geschwindigkeiten |||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 7;
θ0 = π/2;
φ0 = 0;
a  = 9/10;
μ  =-1;
℧  = 0;
q  = 0;
 
(* || Startwerte für die ersten Eigenzeitableitungen eingeben ||||||||||||||||||||||||||| *)
 
vr0 = 0;
vθ0 = 2/Sqrt[61];
vφ0 = 12/(5 Sqrt[61]);
v0  = Sqrt[vr0^2+vθ0^2+vφ0^2];
 
(* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *)
 
ε0 = (Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)]+(a vφ0 (2 r0-℧^2) Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))/Sqrt[1-v0^2];
L0 = (vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^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 ((vφ0^2 (-(a^2+r0^2)^2+a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2))/((-1+v0^2) (r0^2+a^2 Cos[θ0]^2))-a^2 (-1+(Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)]+(a vφ0 (2 r0-℧^2) Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))^2/(1-v0^2)));
 
(* || Benötigte Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
Ξ=r0^2+a^2 Cos[θ0]^2;
δ=r0^2-2r0+a^2;
j[v_]:=Sqrt[1+μ v^2];
pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];
pθ0=vθ0 Sqrt[Ξ]/j[v0];
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code E"
Reduce[ε==ε0 && Lz==L0 && Q==Q0, {ε,Lz,Q}]
N[%]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || Syntax: Mathematica ||| http://kerr.newman.yukterez.net ||| Simon Tyran, Vienna  || *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                                                   
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE F: Erhaltungsgrößen ε, Lz, Q  nach ersten Eigenzeitableitungen  |||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 7;
θ0 = π/2;
φ0 = 0;
a  = 9/10;
μ  =-1;
℧  = 0;
q  = 0;
 
(* || Startwerte für die ersten Eigenzeitableitungen eingeben ||||||||||||||||||||||||||| *)
 
dt = (5 Sqrt[35029/10743])/7;
dr = 0;
dθ = 10/(7 Sqrt[1281]);
dφ = (300 Sqrt[3/125438849])/7+40 Sqrt[3/2136769];
 
(* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *)
 
ε0 = (Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)]+(a vφ0 (2 r0-℧^2) Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))/Sqrt[1-v0^2];
L0 = (vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^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 ((vφ0^2 (-(a^2+r0^2)^2+a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2))/((-1+v0^2) (r0^2+a^2 Cos[θ0]^2))-a^2 (-1+(Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)]+(a vφ0 (2 r0-℧^2) Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))^2/(1-v0^2)));
 
(* || Benötigte Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
Ξ=r0^2+a^2 Cos[θ0]^2;
δ=r0^2-2r0+a^2;
j[v_]:=Sqrt[1+μ v^2];
pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];
pθ0=vθ0 Sqrt[Ξ]/j[v0];
 
dT=(-a^2 (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) (q r0 Sqrt[1-v0^2] ℧+(a^2+r0^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-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2));
dR=-(vr0/((-1+v0^2) Sqrt[-((r0^2+a^2 Cos[θ0]^2)/((-1+v0^2) (a^2-2 r0+r0^2+℧^2)))]));
dΘ=vθ0/(Sqrt[1-v0^2] Sqrt[r0^2+a^2 Cos[θ0]^2]);
dΦ=(Csc[θ0] (a q r0 ℧ Sin[θ0]+(vφ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)])/Sqrt[1-v0^2]+(1/Sqrt[1-v0^2])a (2 r0-℧^2) Sin[θ0] (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 vφ0 (2 r0-℧^2) Sin[θ0])/((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)]))))/((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2));
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code F"
Solve[ε==ε0 && Lz==L0 && Q==Q0 && dT==dt && dR==dr && dΘ==dθ && dΦ==dφ, {ε,Lz,Q}]
N[%]

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

The comments in the codes are in german, but you can translate them with Google Translate or any other translation site if something remains unclear. Most of the code is self explanatory anyway.

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