Kerr Orbits

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

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

Alle Formeln sind in natürlichen Einheiten:

$\95dpi \rm{G=M=c=1}$ , d.h. alle Längen haben die Einheit $\95dpi \rm \ {G M/c^2}$ und Zeiten $\95dpi \rm \ {G M/c^3}$ .

Metrik und Koordinaten

Die der Kürze wegen zusammengefassten Terme sind:

$\95dpi \rm \Sigma =a^2 \cos ^2 \theta +r^2 \ , \ \ \Delta =a^2+r^2-2 r \ , \ \ \chi =\left(a ^2+r^2\right)^2-a ^2 \ \sin ^2 \theta \ \Delta$

Die kovarianten metrischen Koeffizienten sind:

$\95dpi {g_{\rm t t} = {\rm \frac{\Delta-a^2\ \sin^2 \theta}{\Sigma}} \ , \ \ g_{\rm t \phi} = {\rm \frac{2 \ a \ r}{\Sigma} \ \sin^2\theta} \ , \ \ g_{\rm \phi \phi} = {\rm -\frac{\chi \ \sin^2 \theta}{\Sigma}} \ , \ \ g_{\rm r r} = {\rm -\frac{\Sigma}{\Delta}} \ , \ \ g_{\rm \theta \theta} = -\Sigma}$

Kontravariante Metrik-Komponenten:

$\95dpi {g^{\rm t t} = {\rm \frac{\chi}{\Delta \ \Sigma}} \ , \ \ g^{\rm t \phi} = {\rm \frac{2 \ a \ r}{\Delta \ \Sigma}} \ , \ \ g^{\rm \phi \phi} = {\rm -\frac{\Delta-a^2 \ \sin^2 \theta}{\Delta \ \Sigma \ \sin^2 \theta}} \ , \ \ g^{\rm r r} = {\rm -\frac{\Delta}{\Sigma}} \ , \ \ g^{\rm \theta \theta} = -\frac{1}{\Sigma}}$

Dabei steht $\95dpi \rm{a}$ für den Spinparameter. Mit der Transformationsregel in kartesische Koordinaten:

$\95dpi \rm x = \sqrt {r^2 + a^2} \sin\theta \ \cos\phi \ , \ y = \sqrt {r^2 + a^2} \sin\theta \ \sin\phi \ , \ z = r \cos\theta \quad$

lautet das Linienelement in Boyer-Lindquist-Koordinaten:

$\95dpi {\rm d \tau^2 = \left( 1 - \frac{2 \ r}{\Sigma } \right) d t^2 - \frac{\Sigma}{\Delta} \ d r^2 - \Sigma \ d \theta^2 - \frac{\chi}{\Sigma} \ \sin^2 \theta \ d\phi^2 + \frac{4 \ r \ a \sin^2 \theta }{\Sigma } \ dt \ d \phi}$

und mit der Transformation:

$\95dpi {\rm x = \sqrt {r^2 + a^2} \sin\theta \ \cos(\psi-arctan(a/r)) \ , \ y = \sqrt {r^2 + a^2} \sin\theta \ \sin(\psi-arctan(a/r)) \ , \ z = r \cos\theta }$

mit der Koordinatenzeit T und dem Azimuthalwinkel ψ:

$\95dpi {\rm d T=d t+2 \d r \ r/\Delta \ , \ \ T=t+2 \int \frac{r \ d r}{r^2+a^2-2 \ r} \ , \ \ d \psi=d \phi+ d \phi \ a/ \Delta \ , \ \ \psi=\phi+a \ \int \frac{d r}{a^2+r^2-2 r}}$

lautet das Linienelement in Kerr-Schild-Koordinaten:

$\95dpi {\rm {\rm d}\tau^2={\rm d}x^2+{\rm d}y^2+{\rm d}z^2-{\rm d}T^2+\frac{r_{\rm s} \ r^3}{r^4+a^2 \ z^2} \ \left({\rm d}T+\frac{r \ (x \ {\rm d}x+y \ {\rm d}y)}{r^2+a^2}+\frac{a \ (y \ {\rm d}x-x \ {\rm d}y)}{r^2+a^2}+\frac{z \ {\rm d}z}{r}\right)^2}$

Bewegungsgleichungen

Koordinatenzeitableitung nach der Eigenzeit (dt/dτ):

$\95dpi \rm \dot t = \frac{2 \ E \ r \ \left(a^2+r^2\right)-2 \ a \ L_z \ r}{\Delta \ \Sigma }+E = \frac{\varsigma}{\sqrt{1-v^2}}$

Radialkoordinatenableitung (dr/dτ):

$\95dpi \rm \dot r = \frac{\Delta \ p_r }{\Sigma }$

Radiale Impulskomponentenableitung:

$\95dpi \rm \dot p_r = \frac{(r-1) \left(\mu \ \left(a^2+r^2\right)-k\right)+2 \ E^2 \ r \left(a^2+r^2\right)-2 \ a \ E \ L_z +\Delta \ \mu \ r}{\Delta \ \Sigma }-\frac{2 \ p_r ^2 \ (r-1)}{\Sigma }$

Radialimpulskomponente:

$\95dpi \rm p_r = \frac{v_{r}}{\sqrt{1+\mu \ v^2}} \sqrt{\frac{\Sigma}{\Delta}}$

Breitengradableitung (dθ/dτ):

$\95dpi \rm \dot \theta = \frac{ p_{\theta}}{\Sigma }$

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

$\95dpi \rm \dot p_{\theta} = \sin \theta \ \cos \theta \left(\frac{L_z^2}{\sin ^4 \theta }-a^2 \left(E ^2+\mu \right)\right)/ \ \Sigma$

Latitudinaldrehimpulskomponente:

$\95dpi \rm p_{\theta} = \frac{v_{\theta} \ \sqrt{\Sigma}}{\sqrt{1+\mu \ v^2}}$

Längengradableitung (dФ/dτ):

$\95dpi \rm \dot \phi = \frac{2 \ a \ E \ r+ L_z \ \csc ^2 \theta \ (\Sigma -2 r)}{\Delta \ \Sigma }$

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

$\95dpi \rm \dot p_{\phi} = 0$

Longitudinaldrehimpulskomponente:

$\95dpi \rm p_{\phi} = L_z$

Erhaltungsgröße Carter-Konstante:

$\95dpi \rm Q = p_{\theta}^{2} + \cos^{2}\theta \left( a^{2}(\mu^{2} - E^{2}) + \frac{L_z^2}{ \sin^2\theta}\right) = a^2 \ (\mu^{2}-E^2) \ \sin^2 \delta + L_z^2 \ \tan^2 \delta$

Erhaltungsgröße Carter k:

$\95dpi \rm k = a^2 \left(E ^2+\mu \right)+ L_z ^2+Q$

Erhaltungsgröße Gesamtenergie:

$\95dpi {\rm E = \sqrt{\frac{(\Sigma - 2 \ r) \left(\dot{\theta}^2 \ \Delta \ \Sigma +\dot{r}^2 \ \Sigma -\Delta \ \mu \right)}{\Delta \ \Sigma }+\dot{\phi}^2 \ \Delta \ \sin ^2 \theta } = \sqrt{\frac{\Delta \ \Sigma}{(1+\mu \ v^2) \ \chi}} + \Omega \ L_z}$

Erhaltungsgröße Drehimpuls entlang Ф:

$\95dpi \rm L_z = \frac{\sin ^2 \theta \ (\dot{\phi} \ \Delta \ \Sigma - 2 \ a \ E \ r)}{\Sigma -2 \ r} = \frac{v_{\phi} \ \bar R_{\phi} }{\sqrt{1+\mu \ v^2}}$

Frame-Dragging Winkelableitung (dФ/dτ):

$\95dpi \omega = -\frac{g_{t\phi}}{g_{\phi\phi}} = \frac{\rm 2 \ a \ r}{\chi}$

Gravitative Zeitdilatationskomponente (dt/dτ):

$\95dpi \varsigma = \sqrt{g^{\rm t t}} = \sqrt{\frac{\chi}{\Delta \ \Sigma}}$

Axialer und koaxialer Gyrationsradius:

$\95dpi {\rm \bar R_{\phi}} = \sqrt{-g_{\phi \phi}} = \sqrt{\frac{\chi}{\Sigma}} \ \sin \theta \ , \ \ \rm \bar{R}_{\theta} = \sqrt{-g_{\theta \theta}} = \sqrt{\Sigma}$

Axialer und koaxialer Umfang:

$\95dpi \rm U_{\phi} = 2 \ \pi \ \bar R_{\phi} \ , \ \ \rm U_{\theta} = \int_0^{2 \pi } \sqrt{\Sigma } \ {\rm d}\theta$

Lokale 3er-Geschwindigkeit auf der r-Achse:

$\95dpi \rm \frac{v_{r}}{\sqrt{1+\mu \ v^2}} = \dot r \ \sqrt{\frac{\Sigma}{\Delta}}$

Lokale 3er-Geschwindigkeit auf der θ-Achse:

$\95dpi \rm \frac{v_{\theta} \ \sqrt{\Sigma}}{\sqrt{1+\mu \ v^2}} = \dot \theta \ \Sigma$

Lokale 3er-Geschwindigkeit auf der Ф-Achse:

$\95dpi \rm \frac{v_{\phi} \ \bar R_{\phi} }{\sqrt{1+\mu \ v^2}} = \dot \phi \ \Sigma - \omega \ \bar R_{\phi} \ \varsigma$

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)

$\95dpi \rm v_r=v \ \sin \alpha \ , \ \ v_{\phi}=v \ \cos \alpha \ \cos \delta \ , \ \ v_{\theta}=v \ \cos \alpha \ \sin \delta$

Aus der Unendlichkeit beobachtete Geschwindigkeit:

$\95dpi \rm \beta = \sqrt{(dx/dt)^2+(dy/dt)^2+(dz/dt)^2}$

Radiale Fluchtgeschwindigkeit:

$\95dpi \rm \varsigma = \frac{1}{\sqrt{1-v_{\rm esc}^2 }} \ \to \ v_{\rm esc}=\frac{\sqrt{\varsigma ^2-1}}{ \varsigma}$

Code

Code: Alles auswählen

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| Mathematica Syntax | http://kerr.yukterez.net | Version: 21.08.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 [BL]: A=0, kartesisch [KS]: A=a *)

tmax=300;                                                                    (* Eigenzeit *)
Tmax=300;                                                              (* Koordinatenzeit *)
tMax=Min[Tmax, т[plunge-1*^-6]];                                      (* Integrationsende *)

r0=7;                                                                      (* Startradius *)
θ0=π/2;                                                                    (* Breitengrad *)
φ0=0;                                                                       (* Längengrad *)
a=9/10;                                                                  (* Spinparameter *)
μ=-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 *)

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

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

т[τ_]:=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/χ;                 (* 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 *)
v[τ_]:=If[μ==0, 1, Abs[Re[-((\[Sqrt](-a^4(ε-Lz ω[τ])^2-2 a^2R[τ]^2 (ε-Lz ω[τ])^2-
       R[τ]^4(ε-Lz ω[τ])^2+Δi[τ](Σi[τ]+a^2 Sin[Θ[τ]]^2 (ε-
       Lz ω[τ])^2)))/(Sqrt[-(a^2+R[τ]^2)^2+
       a^2 Sin[Θ[τ]]^2 Δi[τ]](ε - Lz ω[τ])))]]];          (* lokale Dreiergeschwindigkeit *)
pΘ[τ_]:=Evaluate[pθ[τ] /. sol][[1]];
pR[τ_]:=Evaluate[pr[τ] /. sol][[1]];
sh[τ_]:=Re[Sqrt[ß[τ]^2-Ω[τ]^2]];
epot[τ_]:=ε+μ-ekin[τ];                                             (* potentielle Energie *)
ekin[τ_]:=If[μ==0, ς[τ], 1/Sqrt[1-v[τ]^2]-1];                       (* kinetische Energie *)

                                                               (* beobachtete Inklination *)
ink0:=б/. Solve[Z'[0]/ю[0] Cos[б]==-Y'[0]/ю[0] Sin[б]&&б>0&&б<2π&&б<δp[r0, a], б][[1]];

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 7) DIFFERENTIALGLEICHUNG |||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

dp= \!\(\*SuperscriptBox[\(Y\),\(Y\)]\); n0[z_]:=Chop[N[z]];
                                         
DGL={
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'[τ]==1/(Σ Δ) (((r[τ]^2+a^2)μ-k)(r[τ]-1)+r[τ] Δ μ+
2r[τ](r[τ]^2+a^2) ε^2-2 a ε Lz)-(2pr[τ]^2 (r[τ]-1))/Σ,
pr[0]==pr0,
pθ'[τ]==(Sin[θ[τ]]Cos[θ[τ]])/Σ (Lz^2/Sin[θ[τ]]^4-a^2 (ε^2+μ)),
pθ[0]==pθ0
};

DG2={
t'[τ]==ε+(2r[τ](r[τ]^2+a^2)ε-2 a r[τ] Lz)/(Σ Δ)+2r'[τ] r[τ]/Δ,
t[0]==0,
r'[τ]==(pr[τ] Δ)/Σ,
r[0]==r0,
θ'[τ]==pθ[τ]/Σ,
θ[0]==θ0,
φ'[τ]==(2 a r[τ] ε+(Σ-2r[τ])Lz Csc[θ[τ]]^2)/(Σ Δ)+r'[τ] a/Δ,
φ[0]==φ0,
pr'[τ]==1/(Σ Δ) (((r[τ]^2+a^2)μ-k)(r[τ]-1)+r[τ] Δ μ+
2r[τ](r[τ]^2+a^2) ε^2-2 a ε Lz)-(2pr[τ]^2 (r[τ]-1))/Σ,
pr[0]==pr0,
pθ'[τ]==(Sin[θ[τ]]Cos[θ[τ]])/Σ (Lz^2/Sin[θ[τ]]^4-a^2 (ε^2+μ)),
pθ[0]==pθ0
};

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

sol=NDSolve[DGL, {t, r, θ, φ, 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 |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

X[τ_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]];            (* kartesisch *)
Y[τ_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]];
Z[τ_]:=Evaluate[r[τ] Cos[θ[τ]]/.sol][[1]];

x[τ_]:=Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]];       (* Plotkoordinaten *)
y[τ_]:=Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]];
z[τ_]:=Z[τ];

XYZ[τ_]:=Sqrt[X[τ]^2+Y[τ]^2+Z[τ]^2]; XY[τ_]:=Sqrt[X[τ]^2+Y[τ]^2];  (* kartesischer Radius *)

Xyz[{x_, y_, z_}, α_]:={x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z};      (* Rotationsmatrix *)
xYz[{x_, y_, z_}, β_]:={x Cos[β]+z Sin[β], y, z Cos[β]-x Sin[β]};
xyZ[{x_, y_, z_}, ψ_]:={x, y Cos[ψ]-z Sin[ψ], y Sin[ψ]+z Cos[ψ]};

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 10) PLOT EINSTELLUNGEN |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

PR=1.2r0;                                                                   (* Plot Range *)
VP={r0, r0, r0};                                                      (* Perspektive x,y,z*)
d1=10;                                                                    (* Schweiflänge *)
plp=50;                                                            (* Flächenplot Details *)
w1l=0; w2l=0; w1r=0; w2r=0;                                          (* Startperspektiven *)
Mrec=100; mrec=10;                                       (* Parametric Plot Subdivisionen *)
imgsize=380;                                                                 (* Bildgröße *)

s[text_]:=Style[text, FontSize->font]; font=11;                            (* Anzeigestil *)

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

display[T_]:=Grid[{
{s[" t coord"], " = ", s[n0[tk]], s["GM/c³"], s[dp]},
{If[μ==0, s[" affineP"], s[" τ propr"]], " = ", s[n0[T]], s["GM/c³"], s[dp]},
{s[" γ 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[" Ř crc.φ"], " = ", s[n0[яi[T]]], s["GM/c²"], s[dp]},
{s[" Σ crc.θ"], " = ", s[n0[Sqrt[Σi[T]]]], s["GM/c²"], s[dp]},
{s[" Δ rad.f"], " = ", 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[" L axial"], " = ", s[n0[Lz]], s["GMm/c"], s[dp]},
{s[" L polar"], " = ", s[n0[pΘ[T]]], s["GMm/c"], s[dp]},
{s[" p r.mom"], " = ", s[n0[pR[T]]], s["mc"], s[dp]},
{s[" R carts"], " = ", s[n0[XYZ[T]]], s["GM/c²"], s[dp]},
{s[" x carts"], " = ", s[n0[X[T]]], s["GM/c²"], s[dp]},
{s[" y carts"], " = ", s[n0[Y[T]]], s["GM/c²"], s[dp]},
{s[" z carts"], " = ", s[n0[Z[T]]], s["GM/c²"], s[dp]},
{s[" ω 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 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 *)
Rasterize[
Show[Graphics3D[{
{PointSize[0.009], Red, Point[
Xyz[xyZ[{x[T], y[T], z[T]}, w1], w2]]}},
ImageSize-> imgsize,
PlotRange-> PR,
SphericalRegion->False,
ImagePadding-> 1],
horizons[A, None, w1, w2],
If[a==0, {},
Graphics3D[{{PointSize[0.009], Purple, Point[
Xyz[xyZ[{
Sin[φ0-ω0 tk+π/2] Sqrt[x0[A]^2+y0[A]^2],
Cos[φ0-ω0 tk+π/2] Sqrt[x0[A]^2+y0[A]^2],
z0[A]}, w1], w2]]}}]],
If[tk==0, {}, If[a==0, {},
ParametricPlot3D[
Xyz[xyZ[{
Sin[φ0-ω0 tt+π/2] Sqrt[x0[A]^2+y0[A]^2],
Cos[φ0-ω0 tt+π/2] Sqrt[x0[A]^2+y0[A]^2],
z0[A]}, w1], w2],
{tt, Max[0, tk-199/100 π/ω0], tk},
PlotStyle -> {Thickness[0.001], Dashed, Purple},
PlotPoints-> Automatic,
MaxRecursion-> mrec]]],
If[tk==0, {},
Block[{$RecursionLimit = Mrec},
ParametricPlot3D[
Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2], {tt, 0, 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}]];

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/3}]

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 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[" Ř crc.φ"], " = ", s[n0[яi[tp]]], s["GM/c²"], s[dp]},
{s[" Σ crc.θ"], " = ", s[n0[Sqrt[Σi[tp]]]], s["GM/c²"], s[dp]},
{s[" Δ rad.f"], " = ", 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[" L axial"], " = ", s[n0[Lz]], s["GMm/c"], s[dp]},
{s[" L polar"], " = ", s[n0[pΘ[tp]]], s["GMm/c"], s[dp]},
{s[" p r.mom"], " = ", s[n0[pR[tp]]], s["mc"], s[dp]},
{s[" R carts"], " = ", s[n0[XYZ[tp]]], s["GM/c²"], s[dp]},
{s[" x carts"], " = ", s[n0[X[tp]]], s["GM/c²"], s[dp]},
{s[" y carts"], " = ", s[n0[Y[tp]]], s["GM/c²"], s[dp]},
{s[" z carts"], " = ", s[n0[Z[tp]]], s["GM/c²"], s[dp]},
{s[" ω 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 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 *)
Rasterize[
Show[Graphics3D[{
{PointSize[0.009], Red, Point[
Xyz[xyZ[{x[tp], y[tp], z[tp]}, w1], w2]]}},
ImageSize-> imgsize,
PlotRange-> PR,
SphericalRegion->False,
ImagePadding-> 1],
horizons[A, None, w1, w2],
If[a==0, {},
Graphics3D[{{PointSize[0.009], Purple, Point[
Xyz[xyZ[{
Sin[φ0-ω0 т[tp]+π/2] Sqrt[x0[A]^2+y0[A]^2],
Cos[φ0-ω0 т[tp]+π/2] Sqrt[x0[A]^2+y0[A]^2],
z0[A]}, w1], w2]]}}]],
If[tk==0, {}, If[a==0, {},
ParametricPlot3D[
Xyz[xyZ[{
Sin[φ0-ω0 т[tt]+π/2] Sqrt[x0[A]^2+y0[A]^2],
Cos[φ0-ω0 т[tt]+π/2] Sqrt[x0[A]^2+y0[A]^2],
z0[A]}, w1], w2],
{tt, Max[0, д[т[tp]-199/100 π/ω0]], tp},
PlotStyle -> {Thickness[0.001], Dashed, Purple},
PlotPoints-> Automatic,
MaxRecursion-> 12]]],
If[tk==0, {},
Block[{$RecursionLimit = Mrec},
ParametricPlot3D[
Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2], {tt, 0, 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}]];

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) ABSPANN MIT PERSPEKTIVISCHER ROTATION ||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

tMAX=д[Tmax];

display[T_]:=Grid[{
{s[" t coord"], " = ", s[n0[т[T]]], 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[" Ř crc.φ"], " = ", s[n0[яi[T]]], s["GM/c²"], s[dp]},
{s[" Σ crc.θ"], " = ", s[n0[Sqrt[Σi[T]]]], s["GM/c²"], s[dp]},
{s[" Δ rad.f"], " = ", 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[" L axial"], " = ", s[n0[Lz]], s["GMm/c"], s[dp]},
{s[" L polar"], " = ", s[n0[pΘ[T]]], s["GMm/c"], s[dp]},
{s[" p r.mom"], " = ", s[n0[pR[T]]], s["mc"], s[dp]},
{s[" R carts"], " = ", s[n0[XYZ[T]]], s["GM/c²"], s[dp]},
{s[" x carts"], " = ", s[n0[X[T]]], s["GM/c²"], s[dp]},
{s[" y carts"], " = ", s[n0[Y[T]]], s["GM/c²"], s[dp]},
{s[" z carts"], " = ", s[n0[Z[T]]], s["GM/c²"], s[dp]},
{s[" ω 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 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}];

plot1c[{xx_, yy_, zz_, tk_, w1_, w2_}]:=
Rasterize[Show[Graphics3D[{
{PointSize[0.009], Red, Point[Xyz[xyZ[{x[tk], y[tk], z[tk]}, w1], w2]]}},
ImageSize -> imgsize,
PlotRange -> {{-PR, PR}, {-PR, PR}, {-PR, PR}},
SphericalRegion -> False,
ImagePadding -> 1],
horizons[A, None, w1, w2],
If[a==0, {},
Graphics3D[{{PointSize[0.009], Purple, Point[Xyz[xyZ[{
Sin[φ0-ω0 т[tk]+π/2] Sqrt[x0[A]^2+y0[A]^2],
Cos[φ0-ω0 т[tk]+π/2] Sqrt[x0[A]^2+y0[A]^2],
z0[A]}, w1], w2]]}}]],
If[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-> 12]],
Block[{$RecursionLimit = Mrec},
ParametricPlot3D[
Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2],
{tt, 0, tk},
PlotStyle -> {Thickness[0.003], Gray},
PlotPoints -> Automatic,
MaxRecursion -> mrec]],
ViewPoint -> {xx, yy, zz}]];

Do[Print[Rasterize[Grid[{{
plot1c[{0, -Infinity, 0, tMAX, w1l+0, w2l+Ц}],
plot1c[{0, 0, Infinity, tMAX, w1r+Ц, w2r+0}],
display[tMAX]}},
Alignment -> Left]]],
{Ц, 0, π, π/5}]

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 14) EXPORTOPTIONEN |||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

(* Export als HTML Dokument *)
(* Export["dateiname.html", EvaluationNotebook[], "GraphicsOutput" -> "PNG"] *)
(* Export direkt als Bildsequenz *)
(* Do[Export["dateiname" <> ToString[Ц*5/π] <> ".png", Rasterize[...]   ], {Ц, 0, 0, Pi}] *)

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||| http://kerr.yukerez.net ||||| Simon Tyran, Vienna ||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

Herleitung aus dem Linienelement:

Code: Alles auswählen

koord={t,r,θ,φ};
Σ=r^2+a^2 Cos[θ]^2;
Δ=r^2-2 r+a^2;
χ=(r^2+a^2)^2-a^2 Sin[θ]^2 Δ;
tt=1-2r/Σ;
rr=-Σ/Δ;
θθ=-Σ;
φφ=-χ/Σ  Sin[θ]^2;
tφ=4a r Sin[θ]^2/Σ;
n=4;
metrik={{tt,0,0,tφ/2},{0,rr,0,0},{0,0,θθ,0},{tφ/2,0,0,φφ}};
metrik//MatrixForm
inversemetrik=Simplify[Inverse[metrik]];
inversemetrik//MatrixForm
christoffel:=Simplify[Table[(1/2)Sum[(inversemetrik[[i,s]])
(D[metrik[[s,j]],koord[[k]]]+D[metrik[[s,k]],koord[[j]]] -D[metrik[[j,k]],koord[[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[x_]:=(((((((x/.t->t[τ])/.r->r[τ])/.θ->θ[τ])/.φ->φ[τ])/.Derivative[1][t[τ]]-> t'[τ])/.Derivative[1][r[τ]]->r'[τ])/.Derivative[1][θ[τ]]-> θ'[τ])/.Derivative[1][φ[τ]]->φ'[τ]
list[x_]:=x[[1]]==x[[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]]]
geodäsie=Simplify[Table[-Sum[christoffel[[i,j,k]] koord[[j]]' koord[[k]]',{j,1,n},{k,1,n}],{i,1,n}]];
bewegungsgleichung=Table[{koord[[i]]''[τ]==rplc[geodäsie[[i]]]},{i,1,n}];
bewegungsgleichung[[1]][[1]]
bewegungsgleichung[[2]][[1]]
bewegungsgleichung[[3]][[1]]
bewegungsgleichung[[4]][[1]]

Anwendung

a=0.646, x0=7, y0=z0=0, i0=π/2-arctan(5/6), v0=0.4

↗

Selbe Bedingungen bei Spinparameter a=0.628

↗

Selbe Bedingungen bei Spinparameter a=0.623

↗

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

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

↗

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

↗

wikipedia · Beitragvon **Yukterez** » Sa 25. Jun 2016, 11:28

Freier Fall aus der lokalen Ruhelage .
r0=4GM/c², θ0=π/4, v0=0

↗

Plunge Orbit mit neutonischer Orbital- und Fluchtgeschwindigkeit.
v0=vθ0=vz0=√(1GM/r0) (1. Teil)
v0=vθ0=vz0=√(2GM/r0) (2. Teil)
r0=3GM/c², θ0=π/4 (beide Teile)

↗

wikipedia · Beitragvon **Yukterez** » Sa 25. Jun 2016, 23:19

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

↗

wikipedia · Beitragvon **Yukterez** » Sa 2. Jul 2016, 23:20

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°

↗

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

↗

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

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

Spinparameter a=1.5 J·c/G/M².
E=0.94104 mc², Lz=2.0127 GMm/c, Q=5.8334 GMm/c

↗

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

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

↗

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:

↗

wikipedia · Beitragvon **Yukterez** » So 3. Jul 2016, 23:19

r=3 für alle a (beobachteter Inklinationswinkel: 90°)

↗

wikipedia · Beitragvon **Yukterez** » Mi 14. Jun 2017, 06:00

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)