Schon wieder eine ganze Seite vollgespammt; eine gute Gelegenheit das nächste Kapitel zu beginnen. Nachdem wir die Wurfparabel nach Newton und Schwarzschild gesehen haben ist es jetzt an der Zeit noch einen Schritt weiter zu gehen und stark rotierende schwarze Löcher im extremen Feld zu betrachten. Zuerst der mathematische Teil:
- Code: Alles auswählen
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;
tmax=300; (* Eigenzeit *)
Tmax=300; (* Koordinatenzeit *)
r0=7; (* Startradius *)
θ0=π/2; (* Breitengrad *)
φ0=0; (* Längengrad *)
a=0.998; (* Spinparameter *)
μ=-1; (* Baryon: μ=-1, Photon: μ=0 *)
v0=1/4; (* Anfangsgeschwindigkeit *)
α0=Pi/4; (* vertikaler Abschusswinkel *)
ψ0=Pi/4; (* Bahninklinationswinkel *)
vr0=v0 Sin[α0]; (* radiale Geschwindigkeitskomponente *)
vφ0=v0 Cos[α0] Sin[ψ0]; (* longitudinale Geschwindigkeitskomponente *)
vθ0=v0 Cos[α0] Cos[ψ0]; (* latitudinale Geschwindigkeitskomponente *)
ε=Sqrt[δ Ξ/((a^2+r0^2)^2-a^2 δ Sin[θ0]^2)]/Sqrt[1-v0^2]+Lz щ;
Lz=vφ0 Sqrt[Ы^2/(1-v0^2)];
pθ0=vθ0 Sqrt[(Ы^2+z0^2)/(1-v0^2)];
pr0=vr0 Sqrt[(Ξ/δ)/(1-v0^2)]; (* Energie und Drehimpulskomponenten *)
j[v_]:=Sqrt[1-v^2]; (* Lorentzfaktor *)
щ=2r0 a/((r0^2+a^2)^2-a^2 (r0^2+a^2-2 r0)Sin[θ0]^2); (* Frame Drag *)
я=Sqrt[((r[τ]^2+a^2)^2-a^2 Δ Sin[θ[τ]]^2)/(r[τ]^2 +a^2 Cos[θ[τ]]^2)]Sin[θ[τ]];
яi[τ_]:=Sqrt[((R[τ]^2+a^2)^2-a^2 Δi[τ] Sin[Θ[τ]]^2)/(R[τ]^2 +a^2 Cos[Θ[τ]]^2)]Sin[Θ[τ]];
Ы=Sqrt[((r0^2+a^2)^2-a^2 δ Sin[θ0]^2)/(r0^2 +a^2 Cos[θ0]^2)]Sin[θ0];
Σ=r[τ]^2+a^2 Cos[θ[τ]]^2; (* zusammengefasste Terme *)
Σ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;
k=pθ[τ]^2+Lz^2 Csc[θ[τ]]^2+a^2 (ε^2 Sin[θ[τ]]^2+μ); (* Carter k *)
x0=Sqrt[r0^2+a^2] Sin[θ0] Cos[φ0];
y0=Sqrt[r0^2+a^2] Sin[θ0] Sin[φ0];
z0=r0 Cos[θ0]; (* kartesische Koordinaten *)
dp= \!\(\*SuperscriptBox[\(Y\),\(Y\)]\); n0[z_]:=If[z<1*^-10, 0, N[z]];
n0[z_]:=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
}; (* Differentialgleichung *)
sol=NDSolve[DGL, {t, r, θ, φ, pr, pθ, ν}, {τ, 0, tmax},
WorkingPrecision-> wp,
MaxSteps-> Infinity,
Method-> mta,
InterpolationOrder-> All]; (* Integrator *)
X[τ_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]];
Y[τ_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]];
Z[τ_]:=Evaluate[r[τ] Cos[θ[τ]]/.sol][[1]];
XYZ[τ_]:=Sqrt[X[τ]^2+Y[τ]^2+Z[τ]^2]; XY[τ_]:=Sqrt[X[τ]^2+Y[τ]^2];
Xyz[{x_, y_, z_}, α_]:={x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z};
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[ψ]};
rE=1+Sqrt[1-a^2 Cos[θ]^2]; (* äußere Ergosphäre *)
RE[A_]:=
{Sqrt[rE^2+A^2] Sin[θ]Cos[φ], Sqrt[rE^2+A^2] Sin[θ]Sin[φ], rE Cos[θ]};
rG=1-Sqrt[1-a^2 Cos[θ]^2]; (* innere Ergosphäre *)
RG[A_]:=
{Sqrt[rG^2+A^2] Sin[θ]Cos[φ], Sqrt[rG^2+A^2] Sin[θ]Sin[φ], rG Cos[θ]};
rA=1+Sqrt[1-a^2]; (* äußerer Horizont *)
RA[A_]:=
{Sqrt[rA^2+A^2] Sin[θ]Cos[φ], Sqrt[rA^2+A^2] Sin[θ]Sin[φ], rA Cos[θ]};
rI=1-Sqrt[1-a^2]; (* innerer Horizont *)
RI[A_]:=
{Sqrt[rI+A^2] Sin[θ]Cos[φ], Sqrt[rI+A^2] Sin[θ]Sin[φ], rI Cos[θ]};
horizons[A_, mesh_]:=Show[
ParametricPlot3D[RE[A], {φ, 0, 2π}, {θ, 0, π},
Mesh->mesh, PlotStyle->Directive[Blue, Opacity[0.15]]],
ParametricPlot3D[RA[A], {φ, 0, 2π}, {θ, 0, π},
Mesh->None, PlotStyle->Directive[Cyan, Opacity[0.15]]],
ParametricPlot3D[RI[A], {φ, 0, 2π}, {θ, 0, π},
Mesh->None, PlotStyle->Directive[Red, Opacity[0.25]]],
ParametricPlot3D[RG[A], {φ, 0, 2π}, {θ, 0, π},
Mesh->None, PlotStyle->Directive[Red, Opacity[0.35]]]];
BLKS:=Grid[{{horizons[a, 35], horizons[0, 35]}}];
т[τ_]:=Evaluate[t[τ]/.sol][[1]]; (* Koordinatenzeit nach Eigenzeit *)
д[ξ_] :=Quiet[Ξ /.FindRoot[т[Ξ]-ξ, {Ξ, 0}]];
T :=Quiet[д[tk]]; (* Eigenzeit nach Koordinatenzeit *)
γ[τ_]:=Evaluate[t'[τ]/.sol][[1]]; (* Anzeige im Display *)
R[τ_]:=Evaluate[r[τ]/.sol][[1]];
Φ[τ_]:=Evaluate[φ[τ]/.sol][[1]];
Θ[τ_]:=Evaluate[θ[τ]/.sol][[1]];
ß[τ_]:=Sqrt[X'[τ]^2+Y'[τ]^2+Z'[τ]^2 ]/γ[τ];
ς[τ_]:=Sqrt[((a^2+R[τ]^2)^2-a^2 (a^2+(R[τ]-2)R[τ])Sin[Θ[τ]]^2)/((a^2+
(R[τ]-2)R[τ])(a^2 Cos[Θ[τ]]^2+R[τ]^2))];
Λ[τ_]:=R[τ]^2+a^2-2 R[τ];
Υ[τ_]:=(R[τ]^2+a^2)^2-a^2 Λ[τ]Sin[Θ[τ]]^2;
ρ[τ_]:=R[τ]^2+a^2 Cos[Θ[τ]]^2;
ω[τ_]:=2R[τ] a/Υ[τ];
Ω[τ_]:=ω[τ] Sqrt[X[τ]^2+Y[τ]^2];
ж[τ_]:=Sqrt[ς[τ]^2-1]/ς[τ];
v[τ_]:=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 ω[τ])))]];
pΘ[τ_]:=Evaluate[pθ[τ] /. sol][[1]];
pR[τ_]:=Evaluate[pr[τ] /. sol][[1]];
sh[τ_]:=Sqrt[ß[τ]^2-Ω[τ]^2];
epot[τ_]:=ε-1-ekin[τ];
ekin[τ_]:=1/Sqrt[1-v[τ]^2];
(* Plot nach Koordinatenzeit *)
display[T_]:=Grid[{
{s[" t coord"], " = ", s[n0[tk]], s["GM/c³"], s[dp]},
{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]]], s["rad"], s[dp]},
{s[" θ lattd"], " = ", s[n0[Θ[T]]], s["rad"], s[dp]},
{s[" r'τ.Σ/Δ"], " = ", s[N[R'[T] Σi[T]/Δi[T]]], s["c"], s[dp]},
{s[" φ'τ*rgy"], " = ", s[n0[Φ'[T] яi[T]]], s["c"], s[dp]},
{s[" θ'τ*rgy"], " = ", s[n0[Θ'[T] Sqrt[яi[T]^2+Z[T]^2]]], s["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[" 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[" v escpe"], " = ", s[n0[ж[T]]], s["c"], s[dp]},
{s[" v delay"], " = ", s[n0[sh[T]]], s["c"], s[dp]},
{s[" ω fdrag"], " = ", s[n0[ω[T]]], s["c³/G/M"], s[dp]},
{s[" v fdrag"], " = ", s[n0[ω[T] яi[T] ς[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 local"], " = ", s[n0[v[T]]], s["c"], s[dp]},
{s[" "], s[" "], s[" "], s[" "]}},
Alignment-> Left, Spacings-> {0, 0}];
PR=1.2r0; (* Plot Range *)
VP={r0, r0, r0}; (* Perspektive x,y,z*)
d1=50; (* Schweiflänge *)
mrec=10; (* Parametric Plot Subdivisionen *)
imgsize=380; (* Bildgröße *)
s[text_]:=Style[text, FontSize->font]; font=11; (* Anzeigestil *)
A=a; (* pseudosphärisch [BL]: A=0, kartesisch [KS]: A=a *)
plot0[{xx_, yy_, zz_, tk_}]:=
Rasterize[
Show[Graphics3D[{
{PointSize[0.007], Red, Point[{X[T], Y[T], Z[T]}]}},
ImageSize-> imgsize,
PlotRange-> PR,
SphericalRegion->False,
ImagePadding-> 1],
horizons[A, None],
ViewPoint-> {xx, yy, zz}]];
plot1[{xx_, yy_, zz_, tk_}]:=
Rasterize[
Show[Graphics3D[{
{PointSize[0.007], Red, Point[{X[T], Y[T], Z[T]}]}},
ImageSize-> imgsize,
PlotRange-> PR,
SphericalRegion->False,
ImagePadding-> 1],
horizons[A, None],
ParametricPlot3D[{X[tt], Y[tt], Z[tt]}, {tt, 0, Max[1*^-16, T-d1]},
PlotStyle-> {Thickness[0.003], Gray},
PlotPoints-> Automatic,
MaxRecursion-> mrec],
ParametricPlot3D[{X[tt], Y[tt], Z[tt]}, {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[{{
plot0[{0, -Infinity, 0, tk}], plot0[{0, 0, Infinity, tk}], display[Quiet[д[tk]]]
}}, Alignment->Left]]],
{tk, 0, 0, 1}]
Do[
Print[Rasterize[Grid[{{
plot1[{0, -Infinity, 0, tk}], plot1[{0, 0, Infinity, tk}], display[Quiet[д[tk]]]
}}, Alignment->Left]]],
{tk, Tmax, Tmax, 10}]
(* Plot nach Eigenzeit *)
display[T_]:=Grid[{
{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]]], s["rad"], s[dp]},
{s[" θ lattd"], " = ", s[n0[Θ[tp]]], s["rad"], s[dp]},
{s[" r'τ.Σ/Δ"], " = ", s[N[R'[tp] Σi[tp]/Δi[tp]]], s["c"], s[dp]},
{s[" φ'τ*rgy"], " = ", s[n0[Φ'[tp] яi[tp]]], s["c"], s[dp]},
{s[" θ'τ*rgy"], " = ", s[n0[Θ'[tp] Sqrt[яi[tp]^2+Z[tp]^2]]], s["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[" 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[" v escpe"], " = ", s[n0[ж[tp]]], s["c"], s[dp]},
{s[" v delay"], " = ", s[n0[sh[tp]]], s["c"], s[dp]},
{s[" ω fdrag"], " = ", s[n0[ω[tp]]], s["c³/G/M"], s[dp]},
{s[" v fdrag"], " = ", s[n0[ω[tp] яi[tp] ς[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 local"], " = ", s[n0[v[tp]]], s["c"], s[dp]},
{s[" "], s[" "], s[" "], s[" "]}},
Alignment-> Left, Spacings-> {0, 0}];
PR=1.2r0; (* Plot Range *)
VP={r0, r0, r0}; (* Perspektive x,y,z *)
d1=50; (* Schweiflänge *)
mrec=10; (* Parametric Plot Subdivisionen *)
imgsize=380; (* Bildgröße *)
s[text_]:=Style[text, FontSize->font]; font=11; (* Anzeigestil *)
A=a; (* pseudosphärisch [BL]: A=0, kartesisch [KS]: A=a *)
plot0[{xx_, yy_, zz_, tk_}]:=
Rasterize[
Show[Graphics3D[{
{PointSize[0.007], Red, Point[{X[tp], Y[tp], Z[tp]}]}},
ImageSize-> imgsize,
PlotRange-> PR,
SphericalRegion->False,
ImagePadding-> 1],
horizons[A, None],
ViewPoint-> {xx, yy, zz}]];
plot1[{xx_, yy_, zz_, tk_}]:=
Rasterize[
Show[Graphics3D[{
{PointSize[0.007], Red, Point[{X[tp], Y[tp], Z[tp]}]}},
ImageSize-> imgsize,
PlotRange-> PR,
SphericalRegion->False,
ImagePadding-> 1],
horizons[A, None],
ParametricPlot3D[{X[tt], Y[tt], Z[tt]}, {tt, 0, Max[1*^-16, tp-d1]},
PlotStyle-> {Thickness[0.003], Gray},
PlotPoints-> Automatic,
MaxRecursion-> mrec],
ParametricPlot3D[{X[tt], Y[tt], Z[tt]}, {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[{{
plot0[{0, -Infinity, 0, tp}], plot0[{0, 0, Infinity, tp}], display[tp]
}}, Alignment->Left]]],
{tp, 0, 0, 1}]
Do[
Print[Rasterize[Grid[{{
plot1[{0, -Infinity, 0, tp}], plot1[{0, 0, Infinity, tp}], display[tp]
}}, Alignment->Left]]],
{tp, tmax, tmax, 10}]
(* http://kerr.yukerez.net *) (* Code by Simon Tyran, Vienna *)
Animationen folgen später. Die Bedienung erfolgt auf die selbe Weise wie auch beim Schwarzschildsimulator; die Anfangsbedingungen werden in der Form {r, θ, ψ} für die Position und {vr, vθ, vψ} für die einzelnen Geschwindigkeitskomponenten bzw. {v0, φ0, δ0} für die automatische Zerlegung in ebendiese nach vertikalem und horizontalem Abschusswinkel eingegeben.
Weiterführende Lektüre zum tieferen Verständnis und Inspirationsquelle für zukünftige Updates:
https://arxiv.org/pdf/0811.3815v1.pdfFür den Fall dass JuRo & Chief die Seite wieder vollgespammt haben bevor ich dazu gekommen sein werde die Plots zu erstellen binde ich hier unten schon mal zwei vorläufig unsichtbare Platzhalter-GIFs ein die dann bei Gelegenheit aktualisiert werden.
Fleißig,
Edit: Code upgedated (Display erweitert, numerische Genauigkeit erhöht, Umschaltung von pseudosphärischer auf kartesische Darstellung ermöglicht, etc)Симон Тыран ↯ Veni, vidi, didici ✲ ♥ yukterez.net