Kerr Metric

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

Kerr Metric

Beitragvon Yukterez » Do 12. Apr 2018, 22:24

Bild

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

Bild
Bild

Here we use natural units of G=M=c=1, so lengths are in GM/c² and times in GM/c³. The metric signature is time-positive (+,-,-,-). a is the spin parameter (for black holes 0≤a≤M), M the mass equivalent of the total energy of the black hole, and Mirr its irreducible mass:

Bild
Bild

Shorthand terms:

Bild

Covariant metric coeffizients:

Bild

Contravariant components (superscripted letters are not powers, but indices):

Bild

The dimensionless spin parameter is a=Jc/G/M². Transformation into cartesian coordinates:

Bild

Line element in Boyer Lindquist coordinates:

Bild

Metric tensor (t,r,θ,Ф):

Bild

Bild

With the transformation:

Bild

where T is a finkelsteinlike time coordinate (radially infalling photons move with dr/dt=1) and ψ the flattened azimuthal angle:

Bild

the metric in Kerr Schild coordinates (T,r,θ,ψ) is:

Bild
Bild

Equations of motion in Boyer Lindquist coordinates:
Bild

Coordinate time by proper time (dt/dτ):

Bild

First proper time derivative of the radial coordinate (dr/dτ):

Bild

Radial momentum derivative:

Bild

Radial momentum:

Bild

Derivative of the poloidial component of motion (dθ/dτ):

Bild

Derivative of the poloidial angular momentum (dpθ/dτ):

Bild

Axial angular momentum:

Bild

Derivative of the axial component of motion (dФ/dτ):

Bild

Axial angular momentum derivative (pФ/dτ):

Bild

Longitudinal component of the angular momentum:

Bild

Constant of motion, Carter's constant:

Bild

Constant of motion, Carter k:

Bild

Constant of motion, total energy:

Bild

Constant of motion, axial angular momentum:

Bild

Local 3-velocity component along the r-axis:

Bild

Local 3-velocity component along the θ-axis:

Bild

Local 3-velocity component along the Ф-axis:

Bild

Local 3-velocity, total:

Bild

For massive testparticles μ=-1 and for photons μ=-0. δ is the inclination angle. With α as the vertical launch anglel the components of the local velocity (relative to a ZAMO) are

Bild

Shapirodelayed and frame dragged velocity as observed at infinity:

Bild

Radial escape velocity:

Bild
Bild

Frame-Dragging angular velocity oberserved at infinity (dФ/dt):

Bild

Delayed Frame-Dragging transverse velocity at the equator of the outer horizon:

Bild

with the horizons and ergospheres (solution for r at Δ=0 and gtt=0):

Bild

r and θ dependend delayed Frame-Dragging transverse velocities:

Bild

at the equatorialen plane at θ=π/2:

Bild

r und θ dependend local Frame-Dragging transverse velocities (greater than c inside of the ergosphere):

Bild

at the equatorialen plane at θ=π/2:

Bild

Cartesian projection of the Frame-Dragging transverse velocity:

Bild

at the equatorialen plane at θ=π/2:

Bild

Gravitational time dilation component relative to a ZAMO (dt/dτ):

Bild

Axial and coaxial radius of gyration:

Bild

Axial and coaxial circumference:

Bild
Bild

For images and animations see the german version of this site.
Bild
Симон Тыран @ vk || wikipedia || stackexchange || wolfram

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

Kerr Metric

Beitragvon Yukterez » Mi 18. Apr 2018, 01:02

Bild

Simulator code, Boyer-Lindquist & Kerr-Schildcoordinates, particles & photons:
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 for the initial conditions;converseion between local velocity, proper time derivatives and constants of motion:

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














Derivation from the line element, Boyer-Lindquist coordinates:

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














Derivation from the line element, Kerr-Schild coordinates:

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














Shadow, 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}]














Shadow, 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
Симон Тыран @ vk || wikipedia || stackexchange || wolfram


Zurück zu „Yukterez Notepad“

Wer ist online?

Mitglieder in diesem Forum: 0 Mitglieder und 1 Gast