Kerr Metrik

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

Kerr Metrik

Beitragvon Yukterez » Mi 22. Jun 2016, 05:12

Bild

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

Verwandte Beiträge: Kerr-Newman || Schwarzschild || Geodätengleichung || Gravitationslinsen || Raytracing || Gummituchmodell
Bild

Alle Formeln sind in den 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 Ruhemasse:

Bild
Bild

Die zusammengefassten Terme sind

Bild

Kovariante metrische Koeffizienten:

Bild

Kontravariante Metrik-Komponenten (hochgestellte Buchstaben sind hierbei keine Potenzen sondern Indizes):

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

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

Verzögerte Frame-Dragging Transversalgeschwindigkeit am Äquator des äußeren Horizonts:

Bild

mit dem Horizont- und Ergosphärenradius (Lösung für r bei Δ=0 und gtt=0):

Bild

Von r und θ abhängige verzögerte Frame-Dragging Transversalgeschwindigkeit:

Bild

Auf der äquatorialen Ebene mit θ=π/2:

Bild

Von r und θ abhängige lokale Frame-Dragging Transversalgeschwindigkeit (ab der Ergosphäre >c):

Bild

Auf der äquatorialen Ebene mit θ=π/2:

Bild

Auf die kartesischen Hintergrundkoordinaten projizierte Frame-Dragging Transversalgeschwindigkeit:

Bild

Auf der äquatorialen Ebene mit θ=π/2:

Bild

Gravitative Zeitdilatationskomponente (dt/dτ):

Bild

Axialer und koaxialer Gyrationsradius:

Bild

Axialer und koaxialer Umfang:

Bild
Bild

Beispiel 1:

Bild

Bahn eines Photons durch die Ringsingularität und Austritt im Antiversum; Startbedingungen: r0=7, θ0=π/4, v0=vr0=c (radial einfallendes Photon). Bahnkonstanten: E total=0.8469 hf0, Carter Q=-0.29 GMhf0/c³, Lz=0 GMhf0/c³, Spin: a=0.9 M.

Untere Hälfte: 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).

Obere Hälfte: die selbe Situation in Kerr-Schild-Koordinaten (die Bahn kann damit bis zur Ringsingularität und darüber hinaus fortgesetzt werden, da die infinite Korotation heraustransformiert und die unendlich vielen Umdrehungen in infinitesimaler Eigenzeit am Horizont dadurch vermieden werden). Animationsparameter: affiner Nullgeodätenparameter.
Bild

Beispiel 2:

Bild

Lokal retrograde Kreisbahn eines Testpartikels im inneren eines mit a=0.99 rotierenden schwarzen Lochs (innerhalb des Cauchyhorizonts, aber außerhalb der inneren Ergosphäre sind wieder radial konstante Orbits möglich, für eine nähere Erklärung siehe hier).

Für weitere Beispiele geht es hier entlang.
Bild
Симон Тыран @ facebook || wikipedia || stackexchange || wolfram

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

Kerr Simulator

Beitragvon Yukterez » Sa 25. Jun 2016, 06:27

Bild

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

Code: Alles auswählen

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| Mathematica Syntax | http://kerr.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 = 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 *)
 
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   = Simplify[Limit[pθ0^2+(Lz^2 Csc[ϑ]^2-a^2 (ε^2+μ)) Cos[ϑ]^2, ϑ->θ0]];     (* Carter Q *)
k   = Q+Lz^2+a^2 (ε^2+μ);                                                     (* Carter k *)
μ   = If[v0==1, 0, -1];                                      (* Baryon: μ=-1, Photon: μ=0 *)
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 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 *)
rG = 1-Sqrt[1-a^2 Cos[θ]^2];                                         (* innere Ergosphäre *)
rA = 1+Sqrt[1-a^2];                                                   (* äußerer Horizont *)
rI = 1-Sqrt[1-a^2];                                                   (* innerer Horizont *)
 
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[A_, w1_, w2_] := Xyz[xyZ[
{Sqrt[rG^2+A^2] Sin[θ] Cos[φ], Sqrt[rG^2+A^2] Sin[θ] Sin[φ], rG Cos[θ]}, w1], w2];
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[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];
ц       = 2r[τ]/Σ; ц0=2r0/Ξ;
 
Σ       = 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 δ;
 
т[τ_] := 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 ]/ю[τ]];
 
gs[τ_] := (2 (R[τ]^2-a^2 Cos[Θ[τ]]^2) Sqrt[((a^2+R[τ]^2)^2-a^2 (a^2+(R[τ]-2) R[τ]) Sin[Θ[τ]]^2)/(a^2+2 R[τ]^2+a^2 Cos[2 Θ[τ]])^2])/(R[τ]^2+a^2 Cos[Θ[τ]]^2)^2;
                                                                           (* Gravitation *)
 
ς[τ_] := Sqrt[Χi[τ]/Δi[τ]/Σi[τ]]; ς0 = Sqrt[χ/δ/Ξ];                     (* gravitative ZD *)
ω[τ_] := 2R[τ] a/Χi[τ];           ω0 = 2r0 a/χ; ωd=2r[τ] a/Χ;           (* Frame Dragging *)
Ω[τ_] := ω[τ] 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[Simplify[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];   (* kartesisches R *)
 
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 *)
d1  = 10;                                                                 (* Schweiflänge *)
plp = 50;                                                          (* Flächenplot Details *)
 
 
 
Mrec    = 100;                                           (* Parametric Plot Subdivisionen *)
mrec    = 10;
 
imgsize = 380;                                                               (* Bildgröße *)
w1l     = 0;                                                 (* Startperspektiven, Winkel *)
w2l     = 0;
w1r     = 0;
w2r     = 0;
 
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]},
{s[" т coord"], " = ", s[n0[tk+Sign[1.5-dsp] 2 NIntegrate[R[Ti]/(R[Ti]^2-2 R[Ti]+a^2),{Ti,0,T}]]], s["GM/c³"], s[dp]},
{If[μ==0, s[" affineP"], s[" τ propr"]], " = ", s[n0[T]], s["GM/c³"], s[dp]},
{s[" γ total"], " = ", s[n0[If[dsp==1, γ[T]-2 If[crd==1, 0, R'[T] R[T]/(R[T]^2-2 R[T]+a^2)], γ[T]]]], s[If[dsp==1, "dt/dτ", "dт/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[" L polar"], " = ", s[n0[If[crd==1, pΘ[T], pΘks[T]]]], 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[n0[Q]], s["GMm/c"], s[dp]},
{s[" a SpinP"], " = ", s[n0[a]], s["GM²/c"], s[dp]},
{s[" d\.b2r/dτ\.b2"], " = ",  s[n0[Evaluate[r''[T] /. sol][[1]]]], s["c⁴/G/M"], s[dp]},
{s[" d\.b2φ/dτ\.b2"], " = ", s[n0[Evaluate[φ''[T] /. sol][[1]]]], s["c⁶/G²/M²"], s[dp]},
{s[" d\.b2θ/dτ\.b2"], " = ", s[n0[Evaluate[θ''[T] /. sol][[1]]]], s["c⁶/G²/M²"], s[dp]},
{s[" α dv/dτ"], " = ", s[n0[Evaluate[vlt''[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, 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[" т coord"], " = ", s[n0[т[tp]+Sign[1.5-dsp] 2 NIntegrate[R[Ti]/(R[Ti]^2-2 R[Ti]+a^2),{Ti,0,T}]]], s["GM/c³"], s[dp]},
{s[" γ total"], " = ", s[n0[If[dsp==1, γ[tp]-2 If[crd==1, 0, R'[tp] R[tp]/(R[tp]^2-2 R[tp]+a^2)], γ[T]]]], s[If[dsp==1, "dt/dτ", "dт/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[" L polar"], " = ", s[n0[If[crd==1, pΘ[T], pΘks[T]]]], 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[n0[Q]], s["GMm/c"], s[dp]},
{s[" a SpinP"], " = ", s[n0[a]], s["GM²/c"], s[dp]},
{s[" d\.b2r/dτ\.b2"], " = ",  s[n0[Evaluate[r''[T] /. sol][[1]]]], s["c⁴/G/M"], s[dp]},
{s[" d\.b2φ/dτ\.b2"], " = ", s[n0[Evaluate[φ''[T] /. sol][[1]]]], s["c⁶/G²/M²"], s[dp]},
{s[" d\.b2θ/dτ\.b2"], " = ", s[n0[Evaluate[θ''[T] /. sol][[1]]]], s["c⁶/G²/M²"], s[dp]},
{s[" α dv/dτ"], " = ", s[n0[Evaluate[vlt''[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, 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}];
 
Quiet[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.yukerez.net ||||| Simon Tyran, Vienna ||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)














Solver für die Startbedingungen; Umrechung von lokaler Geschwindigkeit, Eigenzeitableitungen und Erhaltungsgrößen:

Code: Alles auswählen

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE 1: Lokale Geschwindigkeit nach Erhaltungsgrößen ε, Lz, Q ||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 7;
θ0 = π/2;
φ0 = 0;
a  = 9/10;
μ  =-1;
 
(* || 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) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) 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) 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)+(vφ0^2 Cos[θ0]^2 ((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2))/(r0^2+a^2 Cos[θ0]^2)-a^2 Cos[θ0]^2 (-1+v0^2+(Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))^2))/(1+μ v0^2);
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code 1"
Reduce[ε0==ε && L0==Lz && Q0==Q && vr0^2==v0^2-vφ0^2-vθ0^2 && v0>0, {v0,vr0,vφ0,vθ0}]
N[%]

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax |||| http://kerr.yukterez.net |||| Simon Tyran, Vienna  ||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                                                   
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE 2: Lokale Geschwindigkeit nach ersten Eigenzeitableitungen  |||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 7;
θ0 = π/2;
φ0 = 0;
a  = 9/10;
μ  =-1;
 
(* || 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) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) 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) 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)+(vφ0^2 Cos[θ0]^2 ((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2))/(r0^2+a^2 Cos[θ0]^2)-a^2 Cos[θ0]^2 (-1+v0^2+(Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) 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=ε0+(2r0(r0^2+a^2)ε0-2 a r0 L0)/(Ξ δ);
dR=(pr0 δ)/Ξ;
dΘ=pθ0/Ξ;
dΦ=(2 a r0 ε0+(Ξ-2 r0)L0 Csc[θ0]^2)/(Ξ δ);
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code 2"
Reduce[dT==dt && dR==dr && dΘ==dθ && dΦ==dφ && v0>0, {v0,vr0,vφ0,vθ0}]
N[%]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax |||| http://kerr.yukterez.net |||| Simon Tyran, Vienna  ||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                                                   
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE 3: Erste Eigenzeitableitungen nach Erhaltungsgrößen ε, Lz, Q ||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 7;
θ0 = π/2;
φ0 = 0;
a  = 9/10;
μ  =-1;
 
(* || 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) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) 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) 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)+(vφ0^2 Cos[θ0]^2 ((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2))/(r0^2+a^2 Cos[θ0]^2)-a^2 Cos[θ0]^2 (-1+v0^2+(Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) 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=ε0+(2r0(r0^2+a^2)ε0-2 a r0 L0)/(Ξ δ);
dR=(pr0 δ)/Ξ;
dΘ=pθ0/Ξ;
dΦ=(2 a r0 ε0+(Ξ-2 r0)L0 Csc[θ0]^2)/(Ξ δ);
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code 3"
Solve[ε==ε0 && Lz==L0 && Q==Q0 && dt==dT && dr==dR && dθ==dΘ && dφ==dΦ && v0^2==vr0^2+vθ0^2+vφ0^2, {dt,dr,dθ,dφ}]
N[%]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax |||| http://kerr.yukterez.net |||| Simon Tyran, Vienna  ||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                                                   
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE 4: Erste Eigenzeitableitungen nach lokalen Geschwindigkeiten ||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 7;
θ0 = π/2;
φ0 = 0;
a  = 9/10;
μ  =-1;
 
(* || 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) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) 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) 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)+(vφ0^2 Cos[θ0]^2 ((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2))/(r0^2+a^2 Cos[θ0]^2)-a^2 Cos[θ0]^2 (-1+v0^2+(Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) 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=ε0+(2r0(r0^2+a^2)ε0-2 a r0 L0)/(Ξ δ);
dR=(pr0 δ)/Ξ;
dΘ=pθ0/Ξ;
dΦ=(2 a r0 ε0+(Ξ-2 r0)L0 Csc[θ0]^2)/(Ξ δ);
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code 4"
Reduce[dT==dt && dR==dr && dΘ==dθ && dΦ==dφ, {dt,dr,dθ,dφ}]
N[%]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax |||| http://kerr.yukterez.net |||| Simon Tyran, Vienna  ||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                                                   
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE 5: Erhaltungsgrößen ε, Lz, Q nach lokalen Geschwindigkeiten |||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 7;
θ0 = π/2;
φ0 = 0;
a  = 9/10;
μ  =-1;
 
(* || 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) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) 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) 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)+(vφ0^2 Cos[θ0]^2 ((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2))/(r0^2+a^2 Cos[θ0]^2)-a^2 Cos[θ0]^2 (-1+v0^2+(Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) 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=ε0+(2r0(r0^2+a^2)ε0-2 a r0 L0)/(Ξ δ);
dR=(pr0 δ)/Ξ;
dΘ=pθ0/Ξ;
dΦ=(2 a r0 ε0+(Ξ-2 r0)L0 Csc[θ0]^2)/(Ξ δ);
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code 5"
Reduce[ε==ε0 && Lz==L0 && Q==Q0, {ε,Lz,Q}]
N[%]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax |||| http://kerr.yukterez.net |||| Simon Tyran, Vienna  ||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

(* || *)
    (* || *)
        (* || *)
            (* || *)
                 (* || *)
                      (* || *)
                          (* || *)
                              (* || *)
                                  (* ||*)
                                      (* || *)
                                          (* || *)
                                              (* || *)
                                                   (* || *)
                                                        (* || *)
                                                            (* || *)
                                                                (* || *)
                                                                    (* || *)
                                                                        (* ||*)
                                                                            (* || *)
                                                                                (* || *)
                                                                                    (* || *)
                                                                                   
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || CODE 6: Erhaltungsgrößen ε, Lz, Q  nach ersten Eigenzeitableitungen  |||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]
 
(* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 7;
θ0 = π/2;
φ0 = 0;
a  = 9/10;
μ  =-1;
 
(* || 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) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) 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) 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)+(vφ0^2 Cos[θ0]^2 ((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2))/(r0^2+a^2 Cos[θ0]^2)-a^2 Cos[θ0]^2 (-1+v0^2+(Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) 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=ε0+(2r0(r0^2+a^2)ε0-2 a r0 L0)/(Ξ δ);
dR=(pr0 δ)/Ξ;
dΘ=pθ0/Ξ;
dΦ=(2 a r0 ε0+(Ξ-2 r0)L0 Csc[θ0]^2)/(Ξ δ);
 
(* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *)
 
"Code 6"
Solve[ε==ε0 && Lz==L0 && Q==Q0 && dT==dt && dR==dr && dΘ==dθ && dΦ==dφ, {ε,Lz,Q}]
N[%]

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax |||| http://kerr.yukterez.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}]













Bild
Симон Тыран @ facebook || wikipedia || stackexchange || wolfram

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

Photonenorbits

Beitragvon Yukterez » So 3. Jul 2016, 09:17

Bild
Tatsächlich veschwindender axialer Drehimpuls Lz:
Animationsparameter: Koordinatenzeit t, a=1, r0=1+√2, v0=vz0=vθ0=1, vφ0=vr0=0, θ0=π/2:

Bild  ↗

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

Bild  ↗

Animationsparameter: Spin a, r=3, t=300 (beobachteter Inklinationswinkel: 90°)

Bild  ↗
Bild
Симон Тыран @ facebook || wikipedia || stackexchange || wolfram

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

Partikelorbits

Beitragvon Yukterez » Do 15. Jun 2017, 00:29

Bild
Animationsparameter: Koordinatenzeit t, 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  ↗

Freier Fall aus der lokalen Ruhelage .
Animationsparameter: Koordinatenzeit t, 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  ↗

Selbe Situation in horizontüberschreitenden Kerr-Schild Koordinaten: hier entlang  ↗

Extremer Kerr-Orbit .
Animationsparameter: Koordinatenzeit t, r0=4, θ0=π/2, E=0.935711, Lz=1.5, pθ0=2.59808

Bild  ↗

Zackiger Orbit (retrograd) .
Animationsparameter: Koordinatenzeit t, 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.

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))
Animationsparameter: Koordinatenzeit t

a=0 (Schwarzschild Limit)

Bild  ↗

a=0.1

Bild  ↗

a=0.998 (Thorne Limit)

Bild  ↗

Bahndrehimpulsfreier Wurf mit Fluchtgeschwindigkeit aus dem Inneren der Ergosphäre .
Animationsparameter: Koordinatenzeit t, 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.
Bild
Симон Тыран @ facebook || wikipedia || stackexchange || wolfram

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

Kerr-Flächen

Beitragvon Yukterez » Mi 21. Jun 2017, 06:10

Bild
rH: Ereignishorizonte, rE: Ergosphären (in Boyer-Lindquist Koordinaten)

Bild

r,θ,φ (links) vs x,y,z (rechts):

Bild

Farbcode: hellblau = äußere Ergosphäre, cyan = äußerer Horizont, violett = innerer Horizont, Rot = innere Ergosphäre

Bild

Alle relevanten Flächen; die Ringsingularität befindet sich auf der Äquatorebene der inneren Ergosphäre:

Bild

Der Spinparameter in der oberen Illustration beträgt a=0.99Jc/G/M². In den oberen Plots ist die gravitative Masse M gleich 1 gesetzt; für eine Animation mit Mirr=1 klick hier.

Code:

Code: Alles auswählen

(* Perspektive *)
rE = 1 + Sqrt[1 - a^2 Cos[θ]^2]; RE[A_] := {Sqrt[rE^2 + A^2] Sin[θ] Cos[φ], Sqrt[rE^2 + A^2] Sin[θ] Sin[φ], rE Cos[θ]};
rA = 1 + Sqrt[1 - a^2]; RA[A_] := {Sqrt[rA^2 + A^2] Sin[θ] Cos[φ], Sqrt[rA^2 + A^2] Sin[θ] Sin[φ], rA Cos[θ]};
rI = 1 - Sqrt[1 - a^2]; RI[A_] := {Sqrt[rI^2 + A^2] Sin[θ] Cos[φ], Sqrt[rI^2 + A^2] Sin[θ] Sin[φ], rI Cos[θ]};
rG = 1 - Sqrt[1 - a^2 Cos[θ]^2]; RG[A_] := {Sqrt[rG^2 + A^2] Sin[θ] Cos[φ], Sqrt[rG^2 + A^2] Sin[θ] Sin[φ], rG Cos[θ]};
horizons[A_] := Show[
Graphics3D[{}, SphericalRegion -> True, ViewPoint -> {0, -10 Cos[w] + 10 Sin[w], 10 Cos[w] + 10 Sin[w]}, ImageSize -> 340, PlotRange -> 3],
ParametricPlot3D[RE[A], {φ, 0, 2 π}, {θ, 0, π}, Mesh -> None, PlotStyle -> Directive[Blue, Opacity[0.1]]],
ParametricPlot3D[RA[A], {φ, 0, 2 π}, {θ, 0, π}, Mesh -> None, PlotStyle -> Directive[Cyan, Opacity[0.15]]],
ParametricPlot3D[RI[A], {φ, 0, 2 π}, {θ, 0, π}, Mesh -> None, PlotStyle -> Directive[Red, Opacity[0.2]]],
ParametricPlot3D[RG[A], {φ, 0, 2 π}, {θ, 0, π}, Mesh -> None, PlotStyle -> Directive[Red, Opacity[0.5]]]];
a = 0.99;
Do[Print[Rasterize[Grid[{{horizons[a], horizons[0]}}]]], {w, 0, 2 π, π/125}]

Code: Alles auswählen

(* Spinparameter *)
rE = 1 + Sqrt[1 - a^2 Cos[θ]^2]; RE[A_] := {Sqrt[rE^2 + A^2] Sin[θ] Cos[φ], Sqrt[rE^2 + A^2] Sin[θ] Sin[φ], rE Cos[θ]};
rA = 1 + Sqrt[1 - a^2]; RA[A_] := {Sqrt[rA^2 + A^2] Sin[θ] Cos[φ], Sqrt[rA^2 + A^2] Sin[θ] Sin[φ], rA Cos[θ]};
rI = 1 - Sqrt[1 - a^2]; RI[A_] := {Sqrt[rI^2 + A^2] Sin[θ] Cos[φ], Sqrt[rI^2 + A^2] Sin[θ] Sin[φ], rI Cos[θ]};
rG = 1 - Sqrt[1 - a^2 Cos[θ]^2]; RG[A_] := {Sqrt[rG^2 + A^2] Sin[θ] Cos[φ], Sqrt[rG^2 + A^2] Sin[θ] Sin[φ], rG Cos[θ]};
horizons[A_] := Show[
Graphics3D[{}, SphericalRegion -> True, ViewPoint -> {0, 10, 5}, ImageSize -> 340, PlotRange -> 3],
ParametricPlot3D[RE[A], {φ, 0, 2 π}, {θ, 0, π}, Mesh -> None, PlotStyle -> Directive[Blue, Opacity[0.1]]],
ParametricPlot3D[RA[A], {φ, 0, 2 π}, {θ, 0, π}, Mesh -> None, PlotStyle -> Directive[Cyan, Opacity[0.15]]],
ParametricPlot3D[RI[A], {φ, 0, 2 π}, {θ, 0, π}, Mesh -> None, PlotStyle -> Directive[Red, Opacity[0.2]]],
ParametricPlot3D[RG[A], {φ, 0, 2 π}, {θ, 0, π}, Mesh -> None, PlotStyle -> Directive[Red, Opacity[0.5]]]];
Do[Print[Rasterize[Grid[{{horizons[a], horizons[0]}}]]], {a, 0, 1, 1/250}]
Bild
Симон Тыран @ facebook || wikipedia || stackexchange || wolfram

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

Kreisbahnen um rotierende schwarze Löcher

Beitragvon Yukterez » Fr 23. Jun 2017, 09:48

Bild
Die Bedingung für die radiale Beschleunigung und den radialen Impuls auf einer Kreisbahn ist

Bild

Um die lokale Kreisbahngeschwindigkeit (relativ zum ZAMO) für einen Orbit auf gegebenem r mit gegebenem a zu erhalten wird mit v=vφ, vθ=0, θ=π/2 nach vφ aufgelöst. Diese Gleichung hat 2 Lösungen, wovon die kleinere prograd und die größere retrograd ist:

Bild

Der pro- und retrograde Photonenradius liegt damit bei

Bild

und wenn vφ=0, v=vθ gesetzt wird der poloidiale mit Lz=0 (einem lokalen Inklinationswinkel von 90°, siehe Orbit 1) bei

Bild

Wegen der Korotation lokaler Inertialsysteme fahrt das Photon dann trotz dem es keinen axialen Drehimpuls besitzt mit der lokalen Frame-Dragging-Rate die Ф-Achse entlang:

Bild  ↗

Für ein Photon mit genau so viel axialem Drehimpuls um das Frame-Dragging am Äquator zu kompensieren ergibt sich für alle Spinparameter die Boyer-Lindquist-Radialkoordinate

Bild, also der klassische Schwarzschild-Photonenkreisbahnradius (siehe Orbit 2)

Die zu kompensierende Frame-Dragging-Geschwindigkeit beträgt in dem Fall c/3, weshalb der Inklinationswinkel δ0=π/2+ArcSin(1/3) ist, damit die lokale Geschwindigkeitskomponente des Photons gegen die Ф-Achse vФ nach Pythagoras genau 1/3 ist:

Bild  ↗

Auf r=3 ergibt sich für alle a ein beobachteter äquatorialer Inklinationswinkel: Vergleich

Um den von a und r abhängigen Inklinationswinkel für einen Photonenorbit mit Startposition auf der Äquatorebene zu finden wird v0=1, μ=0, vr0=0, pr'0=0 gesetzt und nach δ0 aufgelöst. Dann ergibt sich:

Bild

mit

Bild

Bild

Bild

Bild

Alle geschlossenen Photonenbahnen haben einen konstanten Boyer-Lindquist-Radius, weswegen die Form der abgefahrenen Sphären in kartesischen Koordinaten ein Ellipsoid ergibt.

Animation aller Photonenorbits für den Spinparameter a=1Jc/G/M² (für eine erweiterte Ansicht das Bild anklicken):

Bild  ↗

Herleitung: hier entlang, weiterführender Beitrag: klick

Pro- und retrograde Kreisbahngeschwindigkeit als Funktion von a und r:

Bild

Schwarzschildlimit:

Bild

Astrophysikalisches Limit:

Bild

Extreme Kerr:

Bild

x-Achse: Boyer-Lindquist r, y-Achse: v lokal

Code:

Code: Alles auswählen

ClearAll["Global`*"]

vretro[a_,r_]:=(a^2+2 a Sqrt[r]+r^2)/(Sqrt[a^2+(-2+r) r] (a-r^(3/2)));
vpro[a_,r_]:=(a^2-2 a Sqrt[r]+r^2)/(Sqrt[a^2+(-2+r) r] (a+r^(3/2)));

rh = 1 + Sqrt[1 - a^2];
r1 = 2 (1 + Cos[2/3 ArcCos[-a]]);
r2 = 2 (1 + Cos[2/3 ArcCos[+a]]);

Do[Print[Rasterize[Grid[{{
Plot[{vpro[a, r], vretro[a, r]}, {r, 1, 10},
PlotRange->{{1, 10}, {-1, 1}},
GridLines->{{rh, r1, r2}, {}},
Frame->True, ImageSize->400]},
{"a"->a}}, Alignment->Left]]],
{a, 0, 1, 0.2}]
Bild
Симон Тыран @ facebook || wikipedia || stackexchange || wolfram

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

ZAMO, ZAVO, FREFO & FIDO

Beitragvon Yukterez » Mo 26. Jun 2017, 17:34

Bild
Die mitrotierenden Punkte zeigen lokal ruhende Inertialsysteme LNRFs (locally nonrotating frames) auf denen ZAMOs (zero angular momentum observers) auf konstantem r,θ mit der Frame-Dragging-Rate mitbewegt werden. Darstellung im Inertialsystem eines ZAVOs (zero angular velocity oberservers) at infinity (mit a=1Jc/G/M²):

Bild

ZAMO (magenta), FIDO (cyan) & FREFO (rot) im Bezugssystem des ZAVO:

Bild  ↗

Benötigte Formeln:
Bild

Framedragging Winkelgeschwindigkeit

Bild

Axialer Gyrationsradius

Bild

Kartesischer Umfangsradius

Bild

Gravitative Zeitdilatation

Bild

Die im System des Buchhalters (ZAVO) beobachtete Transversalgeschwindigkeit der lokal stationären Boje (ZAMO) ist

Bild

während die lokal stationäre Boje (ZAMO) misst dass eine zum Buchhalter stationäre und damit relativ zum Strudel bewegte Boje (FIDO) mit

Bild

an ihr vorbeifliegt.
Bild

weiterführender Beitrag: hier entlang
Bild
Симон Тыран @ facebook || wikipedia || stackexchange || wolfram

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

metrische Abstände

Beitragvon Yukterez » Mo 26. Jun 2017, 20:50

Bild
Bestimmung der radialen physikalischen Distanz:
Bild
Das voll ausgeschriebene Linienelement der Kerr-Metrik in BL-Koordinaten ist

Bild                            
Bild
Bild
Bild

Das wird mit a=0 zur Schwarzschild-Metrik in Schwarzschild-Koordinaten

Bild                            
Bild

Für die radiale Strecke wird die Wurzel des roten Terms (g_rr) von r1 bis r2 integriert, für die poloidiale Strecke die des grünen Terms (g_θθ) von θ1 bis θ2, und für die axiale Strecke die des blauen Terms (g_ФФ) von Ф1 bis Ф2. Wenn die Strecke nicht nur auf einer Achse sondern auf mehreren verläuft oder man den Platz auf einer Fläche bzw. den Inhalt eines Volumens betrachtet muss die gekoppelte Differentialgleichung über alle betroffenen Koordinaten integriert werden.

Beispiel: die radiale Strecke vom Schwarzschildradius r1=rs=2GM/c² bis r2=2.1GM/c², also Δr=0.1, ergibt mit a=0 die physikalische Strecke ΔR=0.901826GM/c²:

Bild

Bild

und mit a=Jc/G/M²=1 je nach Polwinkel knapp über oder unter ΔR≈0.2GM/c²:

Bild

Einbettungsdiagramm:

Bild

Wenn man von einer r-Koordinate zur anderen rechnet gilt: je stärker die Rotation, desto geringer die radiale Längenkontraktion, aber desto größer dafür die axiale. Das kommt daher dass der Horizont sich bei steigender Rotation nach hinten verschiebt.

Rechnen wir nicht vom Schwarzschildradius r1=rs=2 weg, sondern vom Ereignishorizont r1=rH=1+√(1-a²), dann werden die Δr=0.1 mit steigendem a zu einem sehr viel höherem ΔR. Bei a=1 würde die radiale Strecke bis zum Ereignishorizont sogar unendlich, siehe Bardeen 1972:

Bild

was wegen dem a≤0.998-Limit und der komischen Zensur in der Natur aber nicht vorkommt. Plot von a=0 bis a=1:

Bild

Einbettungsdiagramm:

Bild

Mit dem praktischen Maximum von a=0.998 beträgt die Strecke bis zum Ereignishorizont in diesem Beispiel

Bild

Selbes Einbettungsdiagramm mit Fokus auf a=0.998:

Bild

ΔR für r1=rH und r2=r1+0.1 (also Δ=0.1) bei a=0.998 und θ=π/2:

Bild

Originalbeitrag: hier entlang
Bild
Симон Тыран @ facebook || wikipedia || stackexchange || wolfram

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

nächstes Kapitel

Beitragvon Yukterez » Di 18. Jul 2017, 09:15

Bild
Bild
Симон Тыран @ facebook || wikipedia || stackexchange || wolfram

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

Referenzen und weiterführende Links

Beitragvon Yukterez » Do 28. Sep 2017, 16:21

verwandte Themen:
Kerr-Newman || Schwarzschild || Geodätengleichung || Gravitationslinsen || Raytracing || Gummituchmodell

Referenzen, weiterführende Literatur und empfohlene Links:
The Kerr Metric || The Kerr-Newman Metric || Kerr Newman Motion || Periodic Table for Black Hole Orbits || Exact Solution for the Kerr Separatrix 1, 2 || The Spinning Black Hole || Madore || Black Hole Penrose Diagrams || Infinite Redshift Surfaces || Universe in Problems || Tapir 26 || Tapir 27 || Odyssey Edu || Horizon Skimming Orbits || Black Holes Introduction || Rigidly rotating ZAMO surfaces || Global Structure of the Kerr Family || Fuerst & Wu || The Angular Momentum of Kerr Black Holes || Roy Kerr 1963 (Original Work) || Raine || James Bardeen 1972 || Shahar Hod 2012 || Hughes 2000 || Claudel et al || Ed Teo || Nülifer et al || Hanauske || Interstellar || Gyoto || Catalogue of Spacetimes || Stellar Charge || Ynogkm || CERN Geodesics Shortcut || Van der Wijk || The interplay between forces in the Kerr–Newman field || Komissarov: Kerr Schild, Electrodynamics || Binary Kerr Merger || De Mees: Rotating Torus || Chrusciel || Tendex & Vortex || Faraoni || Kerr Dual Merger || Multipole Neutron Stars || Multipole Rotating Stars || Shadows || Ring Singularity || BL vs KS Hydrodynamics || Photon motion in Kerr-DeSitter || G-Ray Shadows || Geokerr || Kerr 599 || Kuchelmeister || Null Geodesics || Axialsymmetric Shadows || Nut Lambda Black Holes || GeoVis || Strong Lensing || Hairy BH Notebooks || Interior mass formula || Sebastien Garmier: cuRRay
Bild
Animations by Simon Tyran, Vienna (Yukterez) - reuse permitted under the Creative Commons License CC BY-SA 4.0
Bild
Симон Тыран @ facebook || wikipedia || stackexchange || wolfram


Zurück zu „Yukterez Notizblock“

Wer ist online?

Mitglieder in diesem Forum: Google [Bot] und 2 Gäste