Seite 1 von 2

Kerr Metrik

Verfasst: Mi 22. Jun 2016, 05:12
von Yukterez
Update: ENGLISH VERSION Bild Bild
Bild

Verwandte Beiträge: Kerr Newman Metrik || Schwarzschild Metrik || Geodätengleichung || Gravitationslinseneffekt
Bild

Alle Formeln sind in natürlichen Einheiten:

G=M=c=1, d.h. alle Längen haben die Einheit GM/c² und Zeiten GM/c³. Die metrische Signatur ist time-positive (+,-,-,-). Dabei ist a der Spinparameter (für schwarze Löcher 0≤a≤M), M das gesamte Massenäquivalent des Zentralkörpers, und Mirr seine unreduzierbare Masse:

Bild
Bild

Metrik und Koordinaten:

Die der Kürze wegen zusammengefassten Terme sind:

Bild

Kovariante metrische Koeffizienten:

Bild

Kontravariante Metrik-Komponenten:

Bild

Dabei steht a für den dimensionslosen Spinparameter Jc/G/M². Mit der Transformationsregel in kartesische Koordinaten:

Bild

lautet das Linienelement in Boyer-Lindquist-Koordinaten:

Bild

Metrischer Tensor (t,r,θ,Ф):

Bild

Bild

Μit der Transformation:

Bild

mit der Koordinatenzeit T und dem Azimuthalwinkel ψ:

Bild

lautet die Metrik in Kerr-Schild-Koordinaten (T,r,θ,ψ):

Bild
Bild

Bewegungsgleichungen in Boyer-Lindquist-Koordinaten:

Koordinatenzeitableitung nach der Eigenzeit (dt/dτ):

Bild

Radialkoordinatenableitung (dr/dτ):

Bild

Radiale Impulskomponentenableitung:

Bild

Radialimpulskomponente:

Bild

Breitengradableitung (dθ/dτ):

Bild

Drehimpulsableitung auf der θ-Achse (dpθ/dτ):

Bild

Latitudinaldrehimpulskomponente:

Bild

Längengradableitung (dФ/dτ):

Bild

Drehimpulsableitung auf der Ф-Achse (pФ/dτ):

Bild

Longitudinaldrehimpulskomponente:

Bild

Erhaltungsgröße Carter-Konstante:

Bild

Erhaltungsgröße Carter k:

Bild

Erhaltungsgröße Gesamtenergie:

Bild

Erhaltungsgröße Drehimpuls entlang Ф:

Bild

Lokale 3er-Geschwindigkeit auf der r-Achse:

Bild

Lokale 3er-Geschwindigkeit auf der θ-Achse:

Bild

Lokale 3er-Geschwindigkeit auf der Ф-Achse:

Bild

Lokale 3er-Geschwindigkeit, total:

Bild

Für massebehaftete Testteilchen gilt μ=-1 und für Photonen μ=-0. a ist der Spinparameter und δ der Bahninklinationswinkel. Mit α als dem vertikalen Abschusswinkel ergeben sich die Komponenten der Geschwindigkeit (relativ zum ZAMO)

Bild

Aus der Unendlichkeit beobachtete Geschwindigkeit:

Bild

Radiale Fluchtgeschwindigkeit:

Bild

Bild

Frame-Dragging Winkelgeschwindigkeit (dФ/dt):

Bild

Gravitative Zeitdilatationskomponente (dt/dτ):

Bild

Axialer und koaxialer Gyrationsradius:

Bild

Axialer und koaxialer Umfang:

Bild
Bild

Links: radialer Einfall in Boyer-Lindquist-Koordinaten (die darstellbare Bahn endet an der Koordinatensingularität am äußeren Ereignishorizont, wo der Partikel mit stehengebliebener Uhr am rotierenden Ereignishorizont einfriert und im System des weit entfernten Beobachters unendlich viele Korotationen ausführt), rechts: die selbe Situation in Kerr-Schild-Koordinaten (die Bahn kann damit bis zur Ringsingularität hin fortgesetzt werden, da die Korotation heraustransformiert und die unendlich vielen Umdrehungen in infinitesimaler Eigenzeit am Horizont dadurch vermieden werden). Der Testpartikel startet relativ zum stationären Beobachter at infinity aus der Ruhelage (die lokale Initialgeschwindigkeit ist somit die negative Frame-Dragging-Geschwindigkeit). Startbedingungen: r=5, θ=π/2, v=vφ=-vZAMO=-√(3/527)·6/5 (lokal), vObserved=0. Perspektive: polare Ansicht, Projektion: kartesisch, Animationsparameter: Eigenzeit τ:

Bild

Bild

Simulator Code, Boyer-Lindquist & Kerr-Schild, Partikel & Photonen:

Code: Alles auswählen

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| Mathematica Syntax | http://kerr.yukterez.net | Version: 30.10.2017  |||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

ClearAll["Global`*"]

mt1={"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}};
mt2={"EventLocator", "Event"-> (r[t]-1000001/1000000 rA)};
mt3={"ImplicitRungeKutta", "DifferenceOrder"-> 20};
mt4={"EquationSimplification"-> "Residual"};
mt0=Automatic;
mta=mt0;
wp=MachinePrecision;

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 1) STARTBEDINGUNGEN EINGEBEN |||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

A=a;                                             (* pseudosphärisch: A=0, kartesisch: A=a *)
crd=1;                                       (* Boyer-Lindquist: crd=1, Kerr-Schild crd=2 *)
dsp=1;                                                                   (* Display Modus *)

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 *)
μ=If[v0==1, 0, -1];                                          (* Baryon: μ=-1, Photon: μ=0 *)

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

x0BL[A_]:=Sqrt[r0^2+A^2] Sin[θ0] Cos[φ0];                           (* Anfangskoordinaten *)
y0BL[A_]:=Sqrt[r0^2+A^2] Sin[θ0] Sin[φ0];
z0[A_]:=r0 Cos[θ0];

x0KS[A_]:=((r0 Cos[φ0]+A Sin[φ0]) Sin[θ0]);
y0KS[A_]:=((r0 Sin[φ0]-A Cos[φ0]) Sin[θ0]);

x0[A_]:=If[crd==1, x0BL[A], x0KS[A]];
y0[A_]:=If[crd==1, y0BL[A], y0KS[A]];

ε=Sqrt[δ Ξ/χ]/j[v0]+Lz ω0;                           (* Energie und Drehimpulskomponenten *)
Lz=vφ0 Ы/j[v0];
pθ0=vθ0 Sqrt[Ξ]/j[v0];
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 *)
k=Q+Lz^2+a^2 (ε^2+μ);                                                         (* Carter k *)

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 3) RADIUS NACH GESCHWINDIGKEIT UND VICE VERSA ||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

rPro=2 (1+Cos[2/3 ArcCos[-a]]);                          (* prograder Photonenorbitradius *)
rRet=2 (1+Cos[2/3 ArcCos[+a]]);                        (* retrograder Photonenorbitradius *)
rTeo=1+2 Sqrt[1-a^3/3] Cos[ArcCos[(1-a^2)/(1-a^2/3)^(3/2)]/3];

δp[r_,a_]:=Quiet[δi/.NSolve[(a^4(-1+r)+2(-3+r)r^4+a^2r(6+r(-5+3 r))+   (* Eq. Ink. Winkel *)
4a Sqrt[a^2+(-2+r)r](a^2+3 r^2)Cos[δi]-a^2(3+r)(a^2+(-2+r)r)Cos[2δi])/(2r(a^2+
(-2+r)r)(r^3+a^2(2+r)))==0&&δi<=π&&δi>=0,δi][[1]]];

vPro=(a^2-2a Sqrt[r0]+r0^2)/(Sqrt[a^2+(-2+r0)r0](a+r0^(3/2)));  (* Kreisgeschwindigkeit + *)
vRet=(a^2+2a Sqrt[r0]+r0^2)/(Sqrt[a^2+(-2+r0)r0](a-r0^(3/2)));  (* Kreisgeschwindigkeit - *)

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 4) HORIZONTE UND ERGOSPHÄREN RADIEN ||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

rE=1+Sqrt[1-a^2 Cos[θ]^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];                                           (* 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];                                                     (* ä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];                                                     (* 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[(Sqrt[1-a^2]+1)/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;
Δi[τ_]:=R[τ]^2-2R[τ]+a^2;
δ=r0^2-2r0+a^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 δ;
ц=2r[τ]/Σ; ц0=2r0/Ξ;

т[τ_]:=Evaluate[t[τ]/.sol][[1]];                        (* Koordinatenzeit nach Eigenzeit *)
д[ξ_] :=Quiet[Ξ /.FindRoot[т[Ξ]-ξ, {Ξ, 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]];
ß[τ_]:=Re[Sqrt[X'[τ]^2+Y'[τ]^2+Z'[τ]^2 ]/ю[τ]];

ς[τ_]:=Sqrt[Χi[τ]/Δi[τ]/Σi[τ]]; ς0=Sqrt[χ/δ/Ξ];                         (* gravitative ZD *)
ω[τ_]:=2R[τ] a/Χi[τ]; ω0=2r0 a/χ; ωd=2r[τ] a/Χ;   (* Frame Dragging Winkelgeschwindigkeit *)
Ω[τ_]:=ω[τ] 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 *)
      
vd[τ_]:=Abs[-((\[Sqrt](-a^4(ε-Lz ωd)^2-2 a^2r[τ]^2 (ε-Lz ωd)^2-
        r[τ]^4(ε-Lz ωd)^2+Δ(Σ+a^2 Sin[θ[τ]]^2 (ε-
        Lz ωd)^2)))/(Sqrt[-(a^2+r[τ]^2)^2+
        a^2 Sin[θ[τ]]^2 Δ](ε - Lz ωd)))];         
    
v[τ_]:=If[μ==0, 1, Evaluate[vlt'[τ]/.sol][[1]]];          (* lokale Dreiergeschwindigkeit *)
dst[τ_]:=Evaluate[str[τ]/.sol][[1]];                                           (* Strecke *)
      
pΘ[τ_]:=Evaluate[pθ[τ] /. sol][[1]]; pΘks[τ_]:=Σi[τ] Θ'[τ];                     (* Impuls *)
pR[τ_]:=Evaluate[pr[τ] /. sol][[1]]; pRks[τ_]:=Σi[τ]/Δi[τ] R'[τ];
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[z]];

(* Boyer-Lindquist-Koordinaten *)

pr2τ[τ_]:=1/(Σ Δ) (((r[τ]^2+a^2)μ-k)(r[τ]-1)+r[τ] Δ μ+2r[τ](r[τ]^2+a^2) ε^2-2 a ε Lz)-(2pr[τ]^2 (r[τ]-1))/Σ;
pθ2τ[τ_]:=(Sin[θ[τ]]Cos[θ[τ]])/Σ (Lz^2/Sin[θ[τ]]^4-a^2 (ε^2+μ));
                                         
DG1={
t'[τ]==ε+(2r[τ](r[τ]^2+a^2)ε-2 a r[τ] Lz)/(Σ Δ),
t[0]==0,
r'[τ]==(pr[τ] Δ)/Σ,
r[0]==r0,
θ'[τ]==pθ[τ]/Σ,
θ[0]==θ0,
φ'[τ]==(2 a r[τ] ε+(Σ-2r[τ])Lz Csc[θ[τ]]^2)/(Σ Δ),
φ[0]==φ0,
pr'[τ]==pr2τ[τ],
pr[0]==pr0,
pθ'[τ]==pθ2τ[τ],
pθ[0]==pθ0,
str'[τ]==vd[τ]/Max[1*^-16, Abs[Sqrt[1-vd[τ]^2]]],
str[0]==0,
vlt'[τ]==vd[τ],
vlt[0]==0
};

(* Kerr-Schild-Koordinaten *)

dr=(pr0 δ)/Ξ; dθ=pθ0/Ξ; dφ=(2a r0 ε+(Ξ-2r0)Lz Csc[θ0]^2)/(Ξ δ)+dr a/δ; dΦ=If[θ0==0, 0, If[θ0==π, 0, dφ]];
φk=φ0+cns/.FindRoot[Sqrt[r0^2+a^2] Cos[φ0+cns]-((r0 Cos[φ0]+a Sin[φ0])),{cns,1}];

DG2={
t''[τ]==(2 ((a^4 Cos[θ[τ]]^4+a^2 Cos[θ[τ]]^2 r[τ]-r[τ]^3-r[τ]^4) r'[τ]^2+r[τ] ((a^2 Cos[θ[τ]]^2-r[τ]^2) t'[τ]^2+r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2-2 a^3 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^3 θ'[τ] φ'[τ]+Sin[θ[τ]]^2 (r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2+a^2 (a^2 Cos[θ[τ]]^2-r[τ]^2) Sin[θ[τ]]^2) φ'[τ]^2+a t'[τ] (a (2 a^2 Cos[θ[τ]]^3 Sin[θ[τ]]+r[τ]^2 Sin[2 θ[τ]]) θ'[τ]+2 (-a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^2 φ'[τ]))+r'[τ] ((a^4 Cos[θ[τ]]^4+2 a^2 Cos[θ[τ]]^2 r[τ]-2 r[τ]^3-r[τ]^4) t'[τ]+a (a r[τ] (2 a^2 Cos[θ[τ]]^3 Sin[θ[τ]]+r[τ]^2 Sin[2 θ[τ]]) θ'[τ]+(-a^4 Cos[θ[τ]]^4-2 a^2 Cos[θ[τ]]^2 r[τ]+2 r[τ]^3+r[τ]^4) Sin[θ[τ]]^2 φ'[τ]))))/(a^2 Cos[θ[τ]]^2+r[τ]^2)^3,
t'[0]==Limit[(ц0 (-dr+a Sin[θ1]^2 dΦ))/(-1+ц0)+\[Sqrt]((1/((-1+ц0)^2 Ξ))(Ξ (dr^2+(-1+ц0) (μ-Ξ dθ^2))+Sin[θ1]^2 dΦ (-2a Ξ dr-(-1+ц0) χ dΦ+ц0^2 a^2 Ξ Sin[θ1]^2 dΦ))), θ1->θ0],
t[0]==0,
r''[τ]==(-8 (a^2 Cos[θ[τ]]^2-r[τ]^2) (a^2 Cos[2 θ[τ]]+r[τ] (2+r[τ])) r'[τ]^2+4 r'[τ] (4 (a^2 Cos[θ[τ]]^2-r[τ]^2) (-2 r[τ]+a^2 Sin[θ[τ]]^2) t'[τ]+2 a^2 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[2 θ[τ]] θ'[τ]-a Sin[θ[τ]]^2 (2 r[τ] (a^2 Cos[θ[τ]]^2 (-4+a^2+a^2 Cos[2 θ[τ]])+2 r[τ] ((2+a^2+a^2 Cos[2 θ[τ]]) r[τ]+r[τ]^3-a^2 Sin[θ[τ]]^2))+a^4 Sin[2 θ[τ]]^2) φ'[τ])+2 (a^2+(-2+r[τ]) r[τ]) (4 (a^2 Cos[θ[τ]]^2-r[τ]^2) t'[τ]^2+4 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2+8 a (-a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^2 t'[τ] φ'[τ]+Sin[θ[τ]]^2 (4 r[τ] ((a^2 Cos[θ[τ]]^2+r[τ]^2)^2-a^2 r[τ] Sin[θ[τ]]^2)+a^4 Sin[2 θ[τ]]^2) φ'[τ]^2))/(8 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3),
r'[0]==dr,
r[0]==r0,
θ''[τ]==(4 a^2 r[τ] Sin[2 θ[τ]] (r'[τ]+t'[τ])^2-8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 r'[τ] θ'[τ]+2 a^2 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[2 θ[τ]] θ'[τ]^2-8 a ((Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[θ[τ]]+r[τ] (a^2+r[τ]^2) Sin[2 θ[τ]]) r'[τ]+r[τ] (a^2+r[τ]^2) Sin[2 θ[τ]] t'[τ]) φ'[τ]+(2 a^6 Cos[θ[τ]]^4+r[τ] (a^4 Cos[θ[τ]]^2 (5+Cos[2 θ[τ]]) r[τ]+2 a^2 (2+Cos[2 θ[τ]]) r[τ]^3+2 r[τ]^5+2 a^2 (a^2 (3+Cos[2 θ[τ]])+4 r[τ]^2) Sin[θ[τ]]^2)) Sin[2 θ[τ]] φ'[τ]^2)/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3),
θ'[0]==dθ,
θ[0]==θ0,
φ''[τ]==If[θ[τ]==0, 0, (4 (a^3 Cos[θ[τ]]^2-a r[τ]^2) r'[τ]^2+4 (a^3 Cos[θ[τ]]^2-a r[τ]^2) t'[τ]^2+4 a r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2-8 (a^2 Cos[θ[τ]]^2+r[τ]^2) (Cot[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2+a^2 r[τ] Sin[2 θ[τ]]) θ'[τ] φ'[τ]+a Sin[θ[τ]]^2 (4 r[τ] ((a^2 Cos[θ[τ]]^2+r[τ]^2)^2-a^2 r[τ] Sin[θ[τ]]^2)+a^4 Sin[2 θ[τ]]^2) φ'[τ]^2+8 a t'[τ] (2 Cot[θ[τ]] r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2) θ'[τ]+a (-a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^2 φ'[τ])+8 r'[τ] ((a^3 Cos[θ[τ]]^2-a r[τ]^2) t'[τ]+a Cot[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2) (a^2 Cos[θ[τ]]^2+r[τ] (2+r[τ])) θ'[τ]-(r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2+a^2 (a^2 Cos[θ[τ]]^2-r[τ]^2) Sin[θ[τ]]^2) φ'[τ]))/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3)],
φ'[0]==dΦ,
φ[0]==φk,
str'[τ]==vd[τ]/Max[1*^-16, Abs[Sqrt[1-vd[τ]^2]]],
str[0]==0,
vlt'[τ]==vd[τ],
vlt[0]==0
};

DGL=If[crd==1, DG1, DG2];

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 8) INTEGRATION |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

sol=

NDSolve[DGL, {t, r, θ, φ, str, vlt, pr, pθ}, {τ, 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 |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

XBL[τ_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]];
YBL[τ_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]];
Z[τ_]:=Evaluate[r[τ] Cos[θ[τ]]/.sol][[1]];

XKS[τ_]:=Evaluate[((r[τ] Cos[φ[τ]]+a Sin[φ[τ]]) Sin[θ[τ]])/.sol][[1]];
YKS[τ_]:=Evaluate[((r[τ] Sin[φ[τ]]-a Cos[φ[τ]]) Sin[θ[τ]])/.sol][[1]];

X[τ_]:=If[crd==1, XBL[τ], XKS[τ]];
Y[τ_]:=If[crd==1, YBL[τ], YKS[τ]];

xBL[τ_]:=Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]];
yBL[τ_]:=Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]];
z[τ_]:=Z[τ];

xKS[τ_]:=Evaluate[((r[τ] Cos[φ[τ]]+A Sin[φ[τ]]) Sin[θ[τ]])/.sol][[1]];
yKS[τ_]:=Evaluate[((r[τ] Sin[φ[τ]]-A Cos[φ[τ]]) Sin[θ[τ]])/.sol][[1]];

x[τ_]:=If[crd==1, xBL[τ], xKS[τ]];
y[τ_]:=If[crd==1, yBL[τ], yKS[τ]];

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[" γ total"], " = ", s[n0[γ[T]]], s["dt/dτ"], s[dp]},
{s[" ς gravt"], " = ", s[n0[ς[T]]], 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[" M irred"], " = ", s[n0[mirr]], s["M"], s[dp]},
{s[" L axial"], " = ", s[n0[Lz]], s["GMm/c"], s[dp]},
{s[" Ř crc.φ"], " = ", s[n0[яi[T]]], s["GM/c²"], s[dp]},
{s[" Σ crc.θ"], " = ", s[n0[Sqrt[Σi[T]]]], s["GM/c²"], 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[N[Q]], s["GMm/c"], s[dp]},
{s[" a SpinP"], " = ", s[n0[a]], s["GM²/c"], s[dp]},
If[dsp==1, {s[" L axial"], " = ", s[n0[Lz]], s["GMm/c"], s[dp]},
{s[" d\.b2r/dτ\.b2"], " = ",  s[n0[Evaluate[r''[T] /. sol][[1]]]], s["c⁴/G/M"], s[dp]}],
If[dsp==1, {s[" L polar"], " = ", s[n0[If[crd==1, pΘ[T], pΘks[T]]]], s["GMm/c"], s[dp]},
{s[" d\.b2φ/dτ\.b2"], " = ", s[n0[Evaluate[φ''[T] /. sol][[1]]]], s["c⁶/G²/M²"], s[dp]}],
If[dsp==1, {s[" p r.mom"], " = ", s[n0[If[crd==1, pR[T], pRks[T]]]], s["mc"], s[dp]},
{s[" d\.b2θ/dτ\.b2"], " = ", s[n0[Evaluate[θ''[T] /. sol][[1]]]], s["c⁶/G²/M²"], 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[ω[T]]], s["c³/G/M"], s[dp]},
{s[" v fdrag"], " = ", s[n0[й[T]]], s["c"], s[dp]},
{s[" Ω fdrag"], " = ", s[n0[Ω[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[crd==1, 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[a==0, {},
Graphics3D[{{PointSize[0.009], Purple, Point[
Xyz[xyZ[{
Sin[-φ0-ц0 a Ξ/χ tk+π/2] Sqrt[x0[A]^2+y0[A]^2],
Cos[-φ0-ц0 a Ξ/χ tk+π/2] Sqrt[x0[A]^2+y0[A]^2],
z0[A]}, w1], w2]]}}]]],
If[crd==1, 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, {}, If[a==0, {},
ParametricPlot3D[
Xyz[xyZ[{
Sin[-φ0-ц0 a Ξ/χ tt+π/2] Sqrt[x0[A]^2+y0[A]^2],
Cos[-φ0-ц0 a Ξ/χ 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, 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, Max[0, T-d1], T},
PlotStyle-> {Thickness[0.004]},
ColorFunction-> Function[{x, y, z, t},
Hue[0, 1, 0.5, 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, 0, 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[" γ total"], " = ", s[n0[γ[tp]]], s["dt/dτ"], s[dp]},
{s[" ς gravt"], " = ", s[n0[ς[tp]]], 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[" M irred"], " = ", s[n0[mirr]], s["M"], s[dp]},
{s[" L axial"], " = ", s[n0[Lz]], s["GMm/c"], s[dp]},
{s[" Ř crc.φ"], " = ", s[n0[яi[tp]]], s["GM/c²"], s[dp]},
{s[" Σ crc.θ"], " = ", s[n0[Sqrt[Σi[tp]]]], s["GM/c²"], 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[N[Q]], s["GMm/c"], s[dp]},
{s[" a SpinP"], " = ", s[n0[a]], s["GM²/c"], s[dp]},
If[dsp==1, {s[" L axial"], " = ", s[n0[Lz]], s["GMm/c"], s[dp]},
{s[" d\.b2r/dτ\.b2"], " = ",  s[n0[Evaluate[r''[tp] /. sol][[1]]]], s["c⁴/G/M"], s[dp]}],
If[dsp==1, {s[" L polar"], " = ", s[n0[If[crd==1, pΘ[tp], pΘks[tp]]]], s["GMm/c"], s[dp]},
{s[" d\.b2φ/dτ\.b2"], " = ", s[n0[Evaluate[φ''[tp] /. sol][[1]]]], s["c⁶/G²/M²"], s[dp]}],
If[dsp==1, {s[" p r.mom"], " = ", s[n0[If[crd==1, pR[tp], pRks[tp]]]], s["mc"], s[dp]},
{s[" d\.b2θ/dτ\.b2"], " = ", s[n0[Evaluate[θ''[tp] /. sol][[1]]]], s["c⁶/G²/M²"], 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[crd==1, 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[a==0, {},
Graphics3D[{{PointSize[0.009], Purple, Point[
Xyz[xyZ[{
Sin[-φ0-ц0 a Ξ/χ т[tp]+π/2] Sqrt[x0[A]^2+y0[A]^2],
Cos[-φ0-ц0 a Ξ/χ т[tp]+π/2] Sqrt[x0[A]^2+y0[A]^2],
z0[A]}, w1], w2]]}}]]],
If[crd==1, 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, {}, If[a==0, {},
ParametricPlot3D[
Xyz[xyZ[{
Sin[-φ0-ц0 a Ξ/χ т[tt]+π/2] Sqrt[x0[A]^2+y0[A]^2],
Cos[-φ0-ц0 a Ξ/χ т[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, 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, Max[0, tp-d1], tp},
PlotStyle-> {Thickness[0.004]},
ColorFunction-> Function[{x, y, z, t},
Hue[0, 1, 0.5, Max[Min[(-tp+(t+d1))/d1, 1], 0]]],
ColorFunctionScaling-> False,
PlotPoints-> Automatic,
MaxRecursion-> mrec]]],
ViewPoint-> {xx, yy, zz}];

Quiet[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}]]

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 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.yukerez.net ||||| Simon Tyran, Vienna ||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

Herleitung aus dem Linienelement, Boyer-Lindquist:

Code: Alles auswählen

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

ClearAll["Global`*"]
℧=0; q=0;
                                               (* 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]]]

Herleitung aus dem Linienelement, Kerr-Schild:

Code: Alles auswählen

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

ClearAll["Global`*"]
                                               (* Input: kovariante metrische Komponenten *)
gtt=-(ц-1);
gtr=-ц;
gtφ=ц a Sin[θ]^2;
gφφ=-χ/Σ Sin[θ]^2;
grφ=a(1+ц)Sin[θ]^2;
grr=-(1+ц);
gθθ=-Σ;
                    (* Koordinaten, Dimensionen, magnetisches Monopol, elektrische Ladung *)
x={t, r, θ, φ}; n=4;
                                                                         "Metrischer Tensor"
metrik={{gtt, gtr, 0, gtφ}, {gtr, grr, 0, grφ}, {0, 0, gθθ, 0}, {gtφ, grφ, 0, gφφ}};
Subscript["g", μσ] -> MatrixForm[metrik]

inversemetrik=Simplify[Inverse[metrik]];
"g"^μσ -> MatrixForm[inversemetrik]
                                                                           (* Abkürzungen *)
Σ=r^2+a^2 Cos[θ]^2;
Δ=r^2-2 r+a^2;
χ=(r^2+a^2)^2-a^2 Sin[θ]^2 Δ;
ц=2r/Σ;
                                                                       "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->φ'[τ];
                                                                      "Bewegungsgleichungen"
geodäsie=Simplify[Table[-Sum[
christoffel[[i, j, k]] x[[j]]' x[[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[FullSimplify[Normal[Solve[
-μ==
gtt Т^2+
grr r\.b4^2+
gθθ θ\.b4^2+
gφφ φ\.b4^2+
2 gtr Т r\.b4+
2 gtφ Т φ\.b4+
2 grφ r\.b4 φ\.b4,
Т,Reals]]][[2]][[1]][[2]]]

Schatten, 2D:

Code: Alles auswählen

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || Mathematica Syntax | BLACK HOLE SHADOWS | kerr.yukterez.net | Version 03.10.2017 || *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

rE=M+Sqrt[M^2-a^2 Cos[θ]^2]; re={Sqrt[rE^2+a^2] Sin[θ],rE Cos[θ]};
rG=M-Sqrt[M^2-a^2 Cos[θ]^2]; rg={Sqrt[rG^2+a^2] Sin[θ],rG Cos[θ]};
rA=M+Sqrt[M^2-a^2]; ra={Sqrt[rA^2+a^2] Sin[θ],rA Cos[θ]};
rI=M-Sqrt[M^2-a^2]; ri={Sqrt[rI^2+a^2] Sin[θ],rI Cos[θ]};
A=Interpolation[{{0,0},{0.1,0.21},{0.2,0.4},{0.3,0.61},{0.4,0.82},{0.5,1.01},{0.6,1.24},{0.7,1.48},
{0.8,1.71},{0.866,1.9},{0.9,2.01},{0.95,2.17},{0.98,2.3},{0.99,2.38},{0.999,2.47},{1,2.5}},
InterpolationOrder -> 1];
B=Interpolation[{{0,Sqrt[27]},{0.1,5.19},{0.2,5.18},{0.3,5.16},{0.4,5.14},{0.5,5.11},{0.6,5.07},
{0.7,5.02},{0.8,4.93},{0.866,4.86},{0.9,4.82},{0.95,4.74},{0.98,4.68},{0.99,4.65},{0.999,4.62},{1,4.6}},
InterpolationOrder -> 1];
pp:=ParametricPlot[{re,rg,ra,ri},{θ,0,2π},
Frame->True,ImageSize->400,PlotRange->{{-8,8},{-6,6}},PlotStyle->{{Black,Thick}}];
cp:=ContourPlot[(x^2+y^2+A[a] x)^2-B[a]^2 (x^2+y^2)==0,{X,-8,8},{y,-6,6},ContourStyle->{Red,Thick}];
a=Sin[w]; M=1; x=-X; Do[Print[Rasterize[Grid[{{Show[pp,cp]},{" a"->N[a]}},
Alignment->Left]]],
{w,0,π/2,π/200}]

Schatten, 3D:

Code: Alles auswählen

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || Mathematica Syntax | BLACK HOLE SHADOWS | kerr.yukterez.net | Version 03.10.2017 || *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

Manipulate[

rE=1+Sqrt[1-a^2 Cos[θ]^2]; re={Sqrt[rE^2+a^2] Sin[θ] Cos[φ], Sqrt[rE^2+a^2] Sin[θ] Sin[φ], rE Cos[θ]};
rG=1-Sqrt[1-a^2 Cos[θ]^2]; rg={Sqrt[rG^2+a^2] Sin[θ] Cos[φ], Sqrt[rG^2+a^2] Sin[θ] Sin[φ], rG Cos[θ]};
rA=1+Sqrt[1-a^2]; ra={Sqrt[rA^2+a^2] Sin[θ] Cos[φ], Sqrt[rA^2+a^2] Sin[θ] Sin[φ], rA Cos[θ]};
rI=1-Sqrt[1-a^2]; ri={Sqrt[rI^2+a^2] Sin[θ] Cos[φ], Sqrt[rI^2+a^2] Sin[θ] Sin[φ], rI Cos[θ]};

rot[{x_, y_, z_}, w_]:={x, y Sin[w]-z Cos[w], y Cos[w]+z Sin[w]};

pp:=Show[
ParametricPlot3D[rot[re, ϑ], {φ, 0, 2 π}, {θ, 0, π},
Mesh -> None, PlotPoints -> 20, PlotStyle -> Directive[Gray, Opacity[0.10]],
ImageSize->400, PlotRange->{{-8, 8}, {-8, 8}, {-6, 6}},
ViewPoint->{0, 20, 0}, ImagePadding->1],
ParametricPlot3D[rot[ra, ϑ], {φ, 0, 2 π}, {θ, 0, π},
Mesh -> None, PlotPoints -> 20, PlotStyle -> Directive[Gray, Opacity[0.15]]],
ParametricPlot3D[rot[ri, ϑ], {φ, 0, 2 π}, {θ, 0, π},
Mesh -> None, PlotPoints -> 20, PlotStyle -> Directive[Gray, Opacity[0.25]]],
ParametricPlot3D[rot[rg, ϑ], {φ, 0, 2 π}, {θ, 0, π},
Mesh -> None, PlotPoints -> 20, PlotStyle -> Directive[Gray, Opacity[0.35]]]];
                       
α=Interpolation[{{0,0},{0.1,0.21},{0.2,0.4},{0.3,0.61},{0.4,0.82},{0.5,1.01},{0.6,1.24},{0.7,1.48},
{0.8,1.71},{0.866,1.9},{0.9,2.01},{0.95,2.17},{0.98,2.3},{0.99,2.38},{0.999,2.47},{1,2.5}},
InterpolationOrder -> 1];

β=Interpolation[{{0,Sqrt[27]},{0.1,5.19},{0.2,5.18},{0.3,5.16},{0.4,5.14},{0.5,5.11},{0.6,5.07},
{0.7,5.02},{0.8,4.93},{0.866,4.86},{0.9,4.82},{0.95,4.74},{0.98,4.68},{0.99,4.65},{0.999,4.62},{1,4.6}},
InterpolationOrder -> 1];

A[a_]:=α[a] Sin[ϑ]+a Sin[ϑ]^3 Cos[ϑ]^2/5;
B[a_]:=β[a]+0.23 Cos[ϑ]^4(1-Sqrt[1-a^4]);

cp:=ContourPlot3D[(x^2+y^2+A[a] x)^2-B[a]^2 (x^2+y^2)==0, {x, -8, 8}, {z, -3, 3}, {y, -6, 6}, ContourStyle->{Black, Thick}, PlotPoints -> 20];

Rasterize[Grid[{{Show[pp, cp]}, {" a"->N[a]}},
Alignment->Left]],
{{ϑ, π/2}, 0, π/2}, {{a, 1}, 0, 1}]

Kerr-Orbits

Verfasst: Sa 25. Jun 2016, 06:27
von Yukterez
a=0.646, x0=7, y0=z0=0, i0=π/2-arctan(5/6), v0=0.4

Bild  ↗

Selbe Bedingungen bei Spinparameter a=0.628

Bild  ↗

Selbe Bedingungen bei Spinparameter a=0.623

Bild  ↗

Prograder gebundener Orbit .
a=0.9, v0=0.4, θ0=π/2, r0=7, i0=arctan(5/6)

Bild  ↗

Mit reflektiertem Inklinationswinkel i0 .
a=0.9, v0=0.4, θ0=π/2, r0=7, i0=π/2-arctan(5/6)

Bild  ↗

schiefer freier Fall

Verfasst: Sa 25. Jun 2016, 11:28
von Yukterez
Freier Fall aus der lokalen Ruhelage .
r0=4GM/c², θ0=π/4, v0=0

Bild  ↗

Plunge Orbit mit neutonischer Orbital- und Fluchtgeschwindigkeit.
v0=vθ0=vz0=√(1GM/r0) → Q=4.5, E=0.749859 (1. Teil)
v0=vθ0=vz0=√(2GM/r0) → Q=18.0, E=1.06046 (2. Teil)
a=0.998, r0=3GM/c², θ0=π/2, Lz=0 (beide Teile)

Bild  ↗

naher Orbit

Verfasst: Sa 25. Jun 2016, 23:19
von Yukterez
Extremer Kerr-Orbit .
r0=4, θ0=π/2, E=0.935711, Lz=1.5, pθ0=2.59808

Bild  ↗

Kerr-Orbits

Verfasst: Sa 2. Jul 2016, 23:20
von Yukterez
Zackiger Orbit (retrograd) .
a=0.95, v0=-0.5, vφ0=-sin(11/50)/2=-0.109115, vθ0=-cos(11/50)/2=-0.487949, E=0.956545, Lz=-0.830327, Q=13.4126, pθ0=-3.66233, pr0=0, r0=6.5, θ0=π/2, i0=π/2-11/50=77.3949°

Bild  ↗

Blumiger Orbit (prograd) .
a=0.9, R0=x0=7, E=0.945711, Lz=0, Q=12.5013

Bild  ↗

Die gestrichelte Linie zeigt die Bahn eines ZAMO auf fixem r0.

a=1.5 (overextremal Kerr)

Verfasst: So 3. Jul 2016, 09:17
von Yukterez
Spinparameter a=1.5 J·c/G/M².
E=0.94104 mc², Lz=2.0127 GMm/c, Q=5.8334 GMm/c

Bild  ↗

Photonenorbits

Verfasst: So 3. Jul 2016, 09:17
von Yukterez
Tatsächlich veschwindender axialer Drehimpuls Lz:
a=1, r0=1+√2, v0=vz0=vθ0=1, vφ0=vr0=0, θ0=π/2:

Bild  ↗

Scheinbar veschwindender, aber negativer axialer Drehimpuls Lz:
a=1, r0=3, v0=1, vφ0=-1/3, vθ0=√8/3, vr0=0, θ0=π/2:

Bild  ↗

Photonenorbits

Verfasst: So 3. Jul 2016, 23:19
von Yukterez
r=3 für alle a (beobachteter Inklinationswinkel: 90°)

Bild  ↗

Zero axial momentum orbits

Verfasst: Mi 14. Jun 2017, 06:00
von Yukterez
Abschuss auf der polaren Achse ohne axialen Drehimpuls (Lz=0) .
Startbedingungen: R0=5, θ0=π/2, v0=vz0=vθ0=51/50·√((1/5)/(1-2/5))

a=0 (Schwarzschild Limit)

Bild  ↗

a=0.1

Bild  ↗

a=0.998 (Thorne Limit)

Bild  ↗

Fluchtgeschwindigkeit

Verfasst: Do 15. Jun 2017, 00:29
von Yukterez
Wurf mit Fluchtgeschwindigkeit aus dem Inneren der Ergosphäre .
a=0.998, r0=1.001 rH, v0=vr0=vesc, vφ0=0, vθ0=0

θ0=π/2:

Bild  ↗

θ0=π/4:

Bild  ↗

θ0=1/1000:

Bild  ↗

weitere Startwinkel zum durchklicken. Fortsetzung auf Seite 2