Schwarzschild Metrik (yukterez.net Backup)

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

Verwandte Beiträge: Kerr Newman || Kerr || De Sitter || Geodäten || Gravitationslinsen || Raytracing || Flamms Paraboloid



Freier Fall in ein schwarzes Loch mit v=-c√(rₛ/r) aus der Perspektive der Freifallers. Winkeldurchmesser des Schattens: klick, Diskussion: 1 und 2


Vollpanorama des mit v=-c√(rₛ/r) hineinfallenden Beobachters im Moment in dem der Horizont überschritten wird


Rot/Blauverschiebungsprofil für den mit der negativen Fluchtgeschwindigkeit fallenden Beobachter im oberen Bild


Aberration in den Augen von drei Beobachtern auf der selben Position r=6GM/c² mit verschiedenen Geschwindigkeitsvektoren


Optische Verzerrung einer Erdkugel mit r=1.0001rs, Beobachter auf R=17.5rs: durch den Gravitationslinseneffekt ist auch die Rückseite sichtbar


links: Strahlenbündel im flachen Minkowskiraum, rechts: Schwarzschild. Entfernung der Lichtquelle von schwarzen Loch: 20GM/c²


Strahlenbündel mit Stoßparameter b=√27≈5.2GM/c² (der scheinbare Radius des schwarzen Lochs im System eines Beobachters at infinity)


Orbit mit Periheldrehung; Startbedingungen: r₀=5, θ0=π/2, v₀=vz₀=vθ₀=51/50·√((1/5)/(1-2/5)).

▥ Metrischer Tensor in Droste Koordinaten {t,r,θ,Ф}:

  

▦ Eingehende Kerr-Schild/Finkelsteinkoordinaten mit ť wobei dt→dť+dr(rs/r)/(1-rs/r)/c:

  

▧ Raindrop (Doran/Gullstrand-Painlevé) Koordinaten, Koordinatenzeit defininert durch die Eigenzeit von Freifallern, dt→dτ+dr√(rs/r)/(1-rs/r)/c:

  

Bewegungsgleichung in Schwarzschildkoordinaten; aufgrund der Kugelsymmetrie kann auf eine Winkelkoordinate reduziert werden so dass √(dθ²+sin²θ dФ²)→dθ. Damit lautet die Lösung der geodätischen Gleichung:







mit v⊥=v cos ζ als die transversale, und v∥=v sin ζ als die radiale lokale Dreiergeschwindigkeit mit dem lokalen Abschusswinkel ζ (von außen erscheint dieser aufgrund der radialen Längenkontraktion flacher). Bei einer rein radialen Bewegung vereinfacht sich die Bewegungsgleichung zu



τ ist die Eigenzeit des Testpartikels, und t die Koordinatenzeit eines Beobachters at infinty. Um die Schalenzeit eines stationären Beobachters auf r zu erhalten wird t einfach mit √(1-rs/r) multipliziert, mit dem Schwarzschildradius rs=2GM/c². Die totale Zeitdilatation ist das multiplikative Produkt aus der gravitativen und der kinematischen Zeitdilatation.

Um von der lokalen Geschwindigkeit auf die beobachtete zu transformieren gilt für die radiale und transversale Komponente:



Der Betrag der lokalen bzw. verzögerten Geschwindigkeit ist dann nach Pythagoras



Der Drehimpuls



und die Gesamtenergie at infinity (die relativistische Masse mal der gravitativen Zeitdilatation)



sind konstant. Die lokale kinetische Energie und die potentielle Energie als Differenz der lokalen Energie zu jener at infinity ist



Effektives Potential:



Die radiale Fluchtgeschwindigkeit um von r₀ nach r₁ zu gelangen ist



womit die radiale Fluchtgeschwindigkeit in die Unendlichkeit (r₁→∞)



ist. Lokale Kreisbahngeschwindigkeit, bei der Photonensphäre auf r=3rs/2 die Lichtgeschwindigkeit:



Fallzeit aus der Ruhelage von r₀ bis r (Eigenzeit):



Koordinatenzeit für den freien Fall von r₀ bis r:



Physikalischer Abstand zwischen zwei verschiedenen r im System des externen stationären Koordinatenbuchhalters:



Abstand vom Horizont zur Singularität im System eines mit der negativen Fluchtgeschwindigkeit einfallenden Freifallers:



in Droste Koordinaten mit grr=-1/(1-rs/r) und v=-√(rs/r) wobei γ=1/√(1-v²/c²) bzw. in Raindrop Koordinaten mit grr=-1 und v=0 wobei γ=1. Bei einem Freifall aus der Ruhelage knapp über dem Horizont ist die aufintegrierte Strecke d=πGM/c². Für die Gleichungen in anderen Koordinaten siehe hier. Simulator (mit a=℧=0 reduziert sich der Kerr Newman Simulator auf Schwarzschild): klick, für die veraltete Version in 2D das untere Feld aufklappen.

Simulator (alte Version; im aktuelleren Kerr-Simulator mit erweitertem Display kann mit a=0, crd=1 in Schwarzschild- und mit a=0, crd=2 in horizontüberschreitenden Finkelsteinkoordinaten simuliert werden):

Code:

(*||| Schwarzschild Simulator |||||||| yukterez.net ||||||||||| Version 13.02.2018 |||||||| Syntax: Mathematica |||*)
(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

G = 1; M = 1; c = 1; rs = 2 G M/c^2;                       (* Einheiten *)
wp = MachinePrecision;                                     (* numerische Genauigkeit *)
set= {"GlobalAdaptive", "MaxErrorIncreases" -> 100, Method -> "GaussKronrodRule"};
para = 20; pstep = 1/3;                                    (* Paraboloid Grid *)

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

j[v_] := Sqrt[1 - v^2/c^2];    J = j[v0];                  (* Gammafaktor *)
k[r_] := Sqrt[1 - rs/r];       κ = k[r0];                  (* Schwarzschildfaktor *)

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

r0 = 3/2 rs;                                               (* Startradius *)
v0 = 7/10 c;                                               (* lokale Startgeschwindigkeit *)
ζ  = Pi/4;                                                 (* Abschusswinkel *)
θ0 = 0;                                                    (* Startwinkel *)
τM = 20 G M/c^3;                                           (* Simulationsdauer *)

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

τMax = Min[τM, plunge-0.0001]; TMax=Min[τM, и[plunge-0.0001]];
(* г = Sqrt[χ^2 + γ^2]; Θ = ArcSin[γ/г]; *)                (* Kartesisch auf Polar *)
vr0 = v0 Sin[ζ] κ/J; vθ0 = v0/r0 Cos[ζ]/J;                 (* Längenkontraktion und Tiefenexpansion *)
d1 = τM/10; d2 = d1; f = 3;                                (* Schweifdauer und Frameanzahl *)

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

sol = NDSolve[{                                            (* Differentialgleichung *)
r''[t] == -((G M)/r[t]^2) + r[t] θ'[t]^2 - (3 G M)/c^2 θ'[t]^2,
r'[0]  == vr0,
r[0]   == r0,
θ''[t] == -((2 r'[t] θ'[t])/r[t]),
θ'[0]  == vθ0,
θ[0]   == θ0,
τ'[t]  == Sqrt[c^2 r[t] + r[t] r'[t]^2 - c^2 rs + r[t]^3 θ'[t]^2 - r[t]^2 rs θ'[t]^2]/(c Sqrt[r[t] - rs] Sqrt[1 - rs/r[t]]),
τ[0]   == 0,
cl'[t] == ((r'[t] / k[r[t]])^2 + (θ'[t] r[t])^2)/c^2,
cl[0]  == 0
}, {r, θ, τ, cl}, {t, 0, τM},
MaxSteps -> Infinity,
Method -> Automatic,
WorkingPrecision -> wp,
InterpolationOrder -> All,
StepMonitor :> (laststep=plunge; plunge=t;
stepsize=plunge-laststep;), Method->{"EventLocator",
"Event" :> (If[stepsize<1*^-5, 0, 1])}];

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

t[ξ_] := Quiet[                                            (* Eigenzeit nach Koordinatenzeit *)
χ /.FindRoot[Evaluate[τ[χ] /. sol][[1]] - ξ, {χ, 0}, WorkingPrecision -> wp, Method -> Automatic]];
Τ := Quiet[t[ι]];

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

grr[r_]:=If[Abs[r]>2, (1-2/Abs[r])^-1, (1-(2Abs[r]^2)/2^3)^-1];
gtt[r_]:=If[Abs[r]>2, (1-2/Abs[r])^-1, 4(3Sqrt[1-2/2]-Sqrt[1-(2Abs[r]^2)/2^3])^-2];

Ȓ[я_]:=Quiet[NIntegrate[Sqrt[Abs[grr[r]]], {r, 0, Abs[я]}, Method -> set, MaxRecursion -> 100]];
Ř[я_]:=Quiet[NIntegrate[Sqrt[Abs[grr[r]-1]], {r, Abs[я], para}, Method -> set, MaxRecursion -> 100]];
ř[я_]:=Quiet[Ř[0]-Ř[я]];
grid[я_]:=Quiet[x/.FindRoot[Ȓ[x]-я, {x, 1}]];

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

x[t_] := (Sin[Evaluate[θ[t] /. sol]] Evaluate[r[t] /. sol])[[1]]
y[t_] := (Cos[Evaluate[θ[t] /. sol]] Evaluate[r[t] /. sol])[[1]]

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

R[t_] := Evaluate[r[t] /. sol][[1]];                       (* Radialkoordinate *)
γ[t_] := Evaluate[τ'[t] /. sol][[1]];                      (* Zeitdilatation *)
и[t_] := Evaluate[τ[t] /. sol][[1]];                       (* Koordinatenzeit *)

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

crθ[t_] := Evaluate[cl'[t] /. sol][[1]];                   (* Celerität *)
vrθ[t_] := crθ[t]/Sqrt[1 + crθ[t]^2];
clr[t_] := Evaluate[r'[t] /. sol][[1]];
clθ[t_] := R[t] Evaluate[θ'[t] /. sol][[1]];

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

vr[t_] := clr[t]/γ[t]/k[R[t]]^2;                           (* lokale Geschwindigkeit, radial *)
vt[t_] := clθ[t]/γ[t]/k[R[t]];                             (* lokale Geschwindigkeit, transversal *)
vp[t_] := Sqrt[vr[t]^2 + vt[t]^2];                         (* lokale Geschwindigkeit, total *)
vc[t_] := Sqrt[vr[t]^2 k[R[t]]^2 + vt[t]^2] k[R[t]];       (* Koordinatengeschwindigkeit, total *)

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

s[text_] := Style[text, FontSize -> font];  font = 11;     (* Stil der Anzeigetafel *)
PR = 2 r0;                                                 (* Plot Range *)

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

Do[Print[                                                  (* Animation nach Eigenzeit *)
Rasterize[Grid[{{Show[

Graphics[{Table[{LightGray, Circle[{0, 0}, grid[n]]}, {n, pstep, para, pstep}]},
Frame -> True, ImageSize -> 360,
PlotRange -> PR, ImagePadding ->  Automatic],
Graphics[{Magenta, Opacity[0.1], Annulus[{0, 0}, {2, para}]}],
Graphics[{White, Opacity[0.7], Disk[{0, 0}, 2]}],
Graphics[{Cyan, Opacity[0.1], Disk[{0, 0}, 2]}],
Graphics[{Black, Dashed, Circle[{0, 0}, r0]}],

Graphics[{PointSize[0.01], Red, Point[{x[т], y[т]}]}],

ParametricPlot[{x[η], y[η]}, {η, 0, Max[1*^-16, т - d1]},
PlotStyle->{Thickness[0.004], LightGray}],
ParametricPlot[{x[η], y[η]}, {η, Max[1*^-16, т - d2], т},
ColorFunction -> Function[{x, y, η},
Hue[0, 1, 0.5, Max[Min[(-т + (η + d2))/d2, 1], 0]]],
ColorFunctionScaling -> False]]},

      {Grid[{
      {s["  Eigenzeit"],           " = ",    s[N[т, 8]],                                      s[" GM/c³"]},
      {s["  Koordinatenzeit"],     " = ",    s[N[Evaluate[τ[т] /. sol][[1]], 8]],             s[" GM/c³"]},
      {s["  Zeitdilatation"],      " = ",    s[N[γ[т], 8]],                                   s[" dt/dτ"]},
      {s["  Winkel"],              " = ",    s[N[Evaluate[(θ[т] /. sol) 180/Pi][[1]], 8]],    s[" Grad"]},
      {s["  Radialkoordinate"],    " = ",    s[N[R[т] , 8]],                                  s[" GM/c²"]},
      {s["  x-Achse"],             " = ",    s[N[x[т], 8]],                                   s[" GM/c²"]},
      {s["  y-Achse"],             " = ",    s[N[y[т], 8]],                                   s[" GM/c²"]},
      {s["  v lokal"],             " = ",    s[N[vp[т], 8]],                                  s[" c"]},
      {s["  v extern"],            " = ",    s[N[vc[т], 8]],                                  s[" c"]},
      {s["  kinetische Energie"],  " = ",    s[N[1/Sqrt[1 - vp[т]^2] - 1, 8]],                s[" mc²"]},
      {s[" "],                     "   ",    s["                 "],                          s[" "]}
      }, Alignment -> Left, Spacings -> {0, 1/2}]}}, Alignment -> Left]]
      ], {т, τMax/f, τMax, τMax/f}]
   
(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

Do[Print[                                                  (* Animation nach Koordinatenzeit *)
Rasterize[Grid[{{Show[

Graphics[{Table[{LightGray, Circle[{0, 0}, grid[n]]}, {n, pstep, para, pstep}]},
Frame -> True, ImageSize -> 360,
PlotRange -> PR, ImagePadding ->  Automatic],
Graphics[{Magenta, Opacity[0.1], Annulus[{0, 0}, {2, para}]}],
Graphics[{White, Opacity[0.7], Disk[{0, 0}, 2]}],
Graphics[{Cyan, Opacity[0.1], Disk[{0, 0}, 2]}],
Graphics[{Black, Dashed, Circle[{0, 0}, r0]}],

Graphics[{PointSize[0.01], Red, Point[{x[Τ], y[Τ]}]}],

ParametricPlot[{x[η], y[η]}, {η, 0, Max[1*^-16, Τ - d1]},
PlotStyle->{Thickness[0.004], LightGray}],
ParametricPlot[{x[η], y[η]}, {η, Max[1*^-16, Τ - d2], Τ},
ColorFunction -> Function[{x, y, η},
Hue[0, 1, 0.5, Max[Min[(-Τ + (η + d2))/d2, 1], 0]]],
ColorFunctionScaling -> False]]},

      {Grid[{
      {s["  Eigenzeit"],           " = ",    s[N[Τ, 8]],                                      s[" GM/c³"]},
      {s["  Koordinatenzeit"],     " = ",    s[N[ι, 8]],                                      s[" GM/c³"]},
      {s["  Zeitdilatation"],      " = ",    s[N[γ[Τ], 8]],                                   s[" dt/dτ"]},
      {s["  Winkel"],              " = ",    s[N[Evaluate[(θ[Τ] /. sol) 180/Pi][[1]], 8]],    s[" Grad"]},
      {s["  Radialkoordinate"],    " = ",    s[N[R[Τ], 8]],                                   s[" GM/c²"]},
      {s["  x-Achse"],             " = ",    s[N[x[Τ], 8]],                                   s[" GM/c²"]},
      {s["  y-Achse"],             " = ",    s[N[y[Τ], 8]],                                   s[" GM/c²"]},
      {s["  v lokal"],             " = ",    s[N[vp[Τ], 8]],                                  s[" c"]},
      {s["  v extern"],            " = ",    s[N[vc[Τ], 8]],                                  s[" c"]},
      {s["  kinetische Energie"],  " = ",    s[N[1/Sqrt[1 - vp[Τ]^2] - 1, 8]],                s[" mc²"]},
      {s[" "],                     "   ",    s["                 "],                          s[" "]}
      }, Alignment -> Left, Spacings -> {0, 1/2}]}}, Alignment -> Left]]
      ], {ι, TMax/f, TMax, TMax/f}]
   
(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
(* Syntax: Mathematica || schwarzschild.yukterez.net || kerr.yukterez.net ||| Simon Tyran - [Симон Тыран] - Vienna *)


Verschiedene Integrationsmethoden (bei Orbits in einem Mindestabstand vom zweifachen EH reicht die Methode Automatic):

Code:

Method -> {"EventLocator", "Event" -> (r[t] - 2.000001)}

Code:

Method -> {"StiffnessSwitching", Method -> {"ExplicitRungeKutta", Automatic}}

Code:

Method -> {"ImplicitRungeKutta", DifferenceOrder" -> 20}

Mehrere Testpartikel:

Code:

(* Schwarzschild Simulator für mehrere Testpartikel ||| schwarzschild.yukterez.net  5 2015 ||| Syntax: Mathematica *)
(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////*)

G = 1; M = 1; c = 1; rs = 2 G M/c^2;                       (* Einheiten *)
wp = MachinePrecision;                                     (* Genauigkeit *)
para = 60; pstep = 1/2;                                    (* Paraboloid Grid *)

(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////*)

j[v_] := Sqrt[1 - v^2/c^2];                                (* Gammafaktor *)
k[r_] := Sqrt[1 - rs/r];                                   (* Schwarzschildfaktor *)

(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////*)

r01 = Sqrt[x01^2 + y01^2];                                 (* Startradius *)
y01 = +1 rs; x01 = +8 rs;                                  (* Position *)
v01 = 0;                                                   (* lokale Startgeschwindigkeit *)
ζ1  = 0;                                                   (* Abschusswinkel *)
θ01 = ArcSin[y01/r01];                                     (* Startwinkel *)

r02 = Sqrt[x02^2 + y02^2];                                 (* Startradius *)
y02 = -1 rs; x02 = +8 rs;                                  (* Position *)
v02 = 0;                                                   (* lokale Startgeschwindigkeit *)
ζ2  = 0;                                                   (* Abschusswinkel *)
θ02 = ArcSin[y02/r02];                                     (* Startwinkel *)

r03 = Sqrt[x03^2 + y03^2];                                 (* Startradius *)
y03 = +1 rs; x03 = +6 rs;                                  (* Position *)
v03 = 0;                                                   (* lokale Startgeschwindigkeit *)
ζ3  = 0;                                                   (* Abschusswinkel *)
θ03 = ArcSin[y03/r03];                                     (* Startwinkel *)

r04 = Sqrt[x04^2 + y04^2];                                 (* Startradius *)
y04 = -1 rs; x04 = +6 rs;                                  (* Position *)
v04 = 0;                                                   (* lokale Startgeschwindigkeit *)
ζ4  = 0;                                                   (* Abschusswinkel *)
θ04 = ArcSin[y04/r04];                                     (* Startwinkel *)

(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////*)

T  = 100;                                                  (* Simulationsdauer *)

(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////*)

sol1 = NDSolve[{                                           (* Differentialgleichung *)
r''[t] == -((G M)/r[t]^2) + r[t] θ'[t]^2 - (3 G M)/c^2 θ'[t]^2,
r'[0]  == 0,
r[0]   == r01,
θ''[t] == -((2 r'[t] θ'[t])/r[t]),
θ'[0]  == 0,
θ[0]   == θ01,
τ'[t]  == Sqrt[c^2 r[t] + r[t] r'[t]^2 - c^2 rs + r[t]^3 θ'[t]^2 - r[t]^2 rs θ'[t]^2]/(c Sqrt[r[t] - rs] Sqrt[1 - rs/r[t]]),
τ[0]   == 0,
cl'[t] == ((r'[t] / k[r[t]])^2 + (θ'[t] r[t])^2)/c^2,
cl[0]  == 0
}, {r, θ, τ, cl}, {t, 0, T}, WorkingPrecision -> wp,
MaxSteps -> Infinity, Method -> Automatic,
InterpolationOrder -> All];

sol2 = NDSolve[{                                            (* Differentialgleichung *)
r''[t] == -((G M)/r[t]^2) + r[t] θ'[t]^2 - (3 G M)/c^2 θ'[t]^2,
r'[0]  == 0,
r[0]   == r02,
θ''[t] == -((2 r'[t] θ'[t])/r[t]),
θ'[0]  == 0,
θ[0]   == θ02,
τ'[t]  == Sqrt[c^2 r[t] + r[t] r'[t]^2 - c^2 rs + r[t]^3 θ'[t]^2 - r[t]^2 rs θ'[t]^2]/(c Sqrt[r[t] - rs] Sqrt[1 - rs/r[t]]),
τ[0]   == 0,
cl'[t] == ((r'[t] / k[r[t]])^2 + (θ'[t] r[t])^2)/c^2,
cl[0]  == 0
}, {r, θ, τ, cl}, {t, 0, T}, WorkingPrecision -> wp,
MaxSteps -> Infinity, Method -> Automatic,
InterpolationOrder -> All];

sol3 = NDSolve[{                                            (* Differentialgleichung *)
r''[t] == -((G M)/r[t]^2) + r[t] θ'[t]^2 - (3 G M)/c^2 θ'[t]^2,
r'[0]  == 0,
r[0]   == r03,
θ''[t] == -((2 r'[t] θ'[t])/r[t]),
θ'[0]  == 0,
θ[0]   == θ03,
τ'[t]  == Sqrt[c^2 r[t] + r[t] r'[t]^2 - c^2 rs + r[t]^3 θ'[t]^2 - r[t]^2 rs θ'[t]^2]/(c Sqrt[r[t] - rs] Sqrt[1 - rs/r[t]]),
τ[0]   == 0,
cl'[t] == ((r'[t] / k[r[t]])^2 + (θ'[t] r[t])^2)/c^2,
cl[0]  == 0
}, {r, θ, τ, cl}, {t, 0, T}, WorkingPrecision -> wp,
MaxSteps -> Infinity, Method -> Automatic,
InterpolationOrder -> All];

sol4 = NDSolve[{                                            (* Differentialgleichung *)
r''[t] == -((G M)/r[t]^2) + r[t] θ'[t]^2 - (3 G M)/c^2 θ'[t]^2,
r'[0]  == 0,
r[0]   == r04,
θ''[t] == -((2 r'[t] θ'[t])/r[t]),
θ'[0]  == 0,
θ[0]   == θ04,
τ'[t]  == Sqrt[c^2 r[t] + r[t] r'[t]^2 - c^2 rs + r[t]^3 θ'[t]^2 - r[t]^2 rs θ'[t]^2]/(c Sqrt[r[t] - rs] Sqrt[1 - rs/r[t]]),
τ[0]   == 0,
cl'[t] == ((r'[t] / k[r[t]])^2 + (θ'[t] r[t])^2)/c^2,
cl[0]  == 0
}, {r, θ, τ, cl}, {t, 0, T}, WorkingPrecision -> wp,
MaxSteps -> Infinity, Method -> Automatic,
InterpolationOrder -> All];

(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////*)

t1[ξ_] := Quiet[                                           (* Eigenzeit nach Koordinatenzeit *)
χ /.FindRoot[Evaluate[τ[χ] /. sol1][[1]] - ξ, {χ, 0},
WorkingPrecision -> wp, Method -> Automatic]];
Τ1 := Quiet[t1[ι]];
t2[ξ_] := Quiet[                                           (* Eigenzeit nach Koordinatenzeit *)
χ /.FindRoot[Evaluate[τ[χ] /. sol2][[1]] - ξ, {χ, 0},
WorkingPrecision -> wp, Method -> Automatic]];
Τ2 := Quiet[t2[ι]];
t3[ξ_] := Quiet[                                           (* Eigenzeit nach Koordinatenzeit *)
χ /.FindRoot[Evaluate[τ[χ] /. sol3][[1]] - ξ, {χ, 0},
WorkingPrecision -> wp, Method -> Automatic]];
Τ3 := Quiet[t3[ι]];
t4[ξ_] := Quiet[                                           (* Eigenzeit nach Koordinatenzeit *)
χ /.FindRoot[Evaluate[τ[χ] /. sol4][[1]] - ξ, {χ, 0},
WorkingPrecision -> wp, Method -> Automatic]];
Τ4 := Quiet[t4[ι]];

(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////*)

u[x_, y_] = Max[2, Sqrt[x^2 + y^2]];                       (* flamm'sches Paraboloid *)
w[x_, y_] = 2 + Integrate[Sqrt[1/(1 - 2/R)], {R, 2, u[x, y]}];
q[x_, y_] = Sqrt[w[x, y]^2 - u[x, y]^2];
grid[n_]  = я /. FindRoot[2 + Sqrt[(-2 + я) я] + Log[-1 + я + Sqrt[(-2 + я) я]] == R, {я, 0}];

(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////*)

x1[t_] := (Sin[Evaluate[θ[t] /. sol1]] Evaluate[r[t] /. sol1])[[1]]
y1[t_] := (Cos[Evaluate[θ[t] /. sol1]] Evaluate[r[t] /. sol1])[[1]]

x2[t_] := (Sin[Evaluate[θ[t] /. sol2]] Evaluate[r[t] /. sol2])[[1]]
y2[t_] := (Cos[Evaluate[θ[t] /. sol2]] Evaluate[r[t] /. sol2])[[1]]

x3[t_] := (Sin[Evaluate[θ[t] /. sol3]] Evaluate[r[t] /. sol3])[[1]]
y3[t_] := (Cos[Evaluate[θ[t] /. sol3]] Evaluate[r[t] /. sol3])[[1]]

x4[t_] := (Sin[Evaluate[θ[t] /. sol4]] Evaluate[r[t] /. sol4])[[1]]
y4[t_] := (Cos[Evaluate[θ[t] /. sol4]] Evaluate[r[t] /. sol4])[[1]]

(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////*)

R1[t_] := Evaluate[r[t] /. sol1][[1]];                     (* Radialkoordinate *)
γ1[t_] := Evaluate[τ'[t] /. sol1][[1]];                    (* Zeitdilatation *)
и1[t_] := Evaluate[τ[t] /. sol1][[1]];                     (* Koordinatenzeit *)

R2[t_] := Evaluate[r[t] /. sol2][[1]];                     (* Radialkoordinate *)
γ2[t_] := Evaluate[τ'[t] /. sol2][[1]];                    (* Zeitdilatation *)
и2[t_] := Evaluate[τ[t] /. sol2][[1]];                     (* Koordinatenzeit *)

R3[t_] := Evaluate[r[t] /. sol3][[1]];                     (* Radialkoordinate *)
γ3[t_] := Evaluate[τ'[t] /. sol3][[1]];                    (* Zeitdilatation *)
и3[t_] := Evaluate[τ[t] /. sol3][[1]];                     (* Koordinatenzeit *)

R4[t_] := Evaluate[r[t] /. sol4][[1]];                     (* Radialkoordinate *)
γ4[t_] := Evaluate[τ'[t] /. sol4][[1]];                    (* Zeitdilatation *)
и4[t_] := Evaluate[τ[t] /. sol4][[1]];                     (* Koordinatenzeit *)

(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////*)

crθ1[t_] := Evaluate[cl'[t] /. sol1][[1]];                 
vrθ1[t_] := crθ1[t]/Sqrt[1 + crθ1[t]^2];
clr1[t_] := Evaluate[r'[t] /. sol1][[1]];
clθ1[t_] := R1[t] Evaluate[θ'[t] /. sol1][[1]];

crθ2[t_] := Evaluate[cl'[t] /. sol2][[1]];                 
vrθ2[t_] := crθ2[t]/Sqrt[1 + crθ2[t]^2];
clr2[t_] := Evaluate[r'[t] /. sol2][[1]];
clθ2[t_] := R2[t] Evaluate[θ'[t] /. sol2][[1]];

crθ3[t_] := Evaluate[cl'[t] /. sol3][[1]];                 
vrθ3[t_] := crθ3[t]/Sqrt[1 + crθ3[t]^2];
clr3[t_] := Evaluate[r'[t] /. sol3][[1]];
clθ3[t_] := R3[t] Evaluate[θ'[t] /. sol3][[1]];

crθ4[t_] := Evaluate[cl'[t] /. sol4][[1]];                 
vrθ4[t_] := crθ4[t]/Sqrt[1 + crθ4[t]^2];
clr4[t_] := Evaluate[r'[t] /. sol4][[1]];
clθ4[t_] := R4[t] Evaluate[θ'[t] /. sol4][[1]];

(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////*)

vr1[t_] := clr1[t]/γ1[t]/k[R1[t]]^2;                        (* lokale Geschwindigkeit, radial *)
vt1[t_] := clθ1[t]/γ1[t]/k[R1[t]];                          (* lokale Geschwindigkeit, transversal *)
vp1[t_] := Sqrt[vr1[t]^2 + vt1[t]^2];                       (* lokale Geschwindigkeit, total *)
vc1[t_] := Sqrt[vr1[t]^2 k[R1[t]]^2 + vt1[t]^2] k[R1[t]];   (* Koordinatengeschwindigkeit, total *)

vr2[t_] := clr2[t]/γ2[t]/k[R2[t]]^2;                        (* lokale Geschwindigkeit, radial *)
vt2[t_] := clθ2[t]/γ2[t]/k[R2[t]];                          (* lokale Geschwindigkeit, transversal *)
vp2[t_] := Sqrt[vr2[t]^2 + vt2[t]^2];                       (* lokale Geschwindigkeit, total *)
vc2[t_] := Sqrt[vr2[t]^2 k[R2[t]]^2 + vt2[t]^2] k[R2[t]];   (* Koordinatengeschwindigkeit, total *)

vr3[t_] := clr3[t]/γ3[t]/k[R3[t]]^2;                        (* lokale Geschwindigkeit, radial *)
vt3[t_] := clθ3[t]/γ3[t]/k[R3[t]];                          (* lokale Geschwindigkeit, transversal *)
vp3[t_] := Sqrt[vr3[t]^2 + vt3[t]^2];                       (* lokale Geschwindigkeit, total *)
vc3[t_] := Sqrt[vr3[t]^2 k[R3[t]]^2 + vt3[t]^2] k[R3[t]];   (* Koordinatengeschwindigkeit, total *)

vr4[t_] := clr4[t]/γ4[t]/k[R4[t]]^2;                        (* lokale Geschwindigkeit, radial *)
vt4[t_] := clθ4[t]/γ4[t]/k[R4[t]];                          (* lokale Geschwindigkeit, transversal *)
vp4[t_] := Sqrt[vr4[t]^2 + vt4[t]^2];                       (* lokale Geschwindigkeit, total *)
vc4[t_] := Sqrt[vr4[t]^2 k[R4[t]]^2 + vt4[t]^2] k[R4[t]];   (* Koordinatengeschwindigkeit, total *)

(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////*)

s[text_] := Style[text, FontSize -> font];  font = 11;

Do[Print[                                            (* Animation nach Koordinatenzeit *)
  Rasterize[Grid[{{Show[Graphics[{
      {LightGray, Disk[{0, 0}, rs]}},
       Frame -> True, ImageSize -> 400, PlotRange -> 1.2 Max[r01, r02, r03, r04], ImagePadding -> 1],
       Graphics[Table[{Lighter[Gray], Circle[{0, 0}, grid[n]]}, {n, 2 + pstep, para, pstep}]],
      
       Graphics[{PointSize[0.01], Red, Point[{x1[Τ1], y1[Τ1]}]}],
       Graphics[{PointSize[0.01], Green, Point[{x2[Τ2], y2[Τ2]}]}],
       Graphics[{PointSize[0.01], Blue, Point[{x3[Τ3], y3[Τ3]}]}],
       Graphics[{PointSize[0.01], Pink, Point[{x4[Τ4], y4[Τ4]}]}]
      
      ]},
        {Grid[{
      {"  ", s["Koordinatenzeit"], " = ", s[N[ι, 8]], s["    GM/c³"]},
      {""},
      {"  ", s["Eigenzeit 1"], " = ", s[N[Τ1, 8]], s["    GM/c³"]},
      {"  ", s["Eigenzeit 2"], " = ", s[N[Τ2, 8]], s["    GM/c³"]},
      {"  ", s["Eigenzeit 3"], " = ", s[N[Τ3, 8]], s["    GM/c³"]},
      {"  ", s["Eigenzeit 4"], " = ", s[N[Τ4, 8]], s["    GM/c³"]},
      {""},
      {"  ", s["Zeitdilatation 1"], " = ", s[N[Evaluate[τ'[Τ1] /. sol1][[1]], 8]], s["    dt/dτ"]},
      {"  ", s["Zeitdilatation 2"], " = ", s[N[Evaluate[τ'[Τ2] /. sol2][[1]], 8]], s["    dt/dτ"]},
      {"  ", s["Zeitdilatation 3"], " = ", s[N[Evaluate[τ'[Τ3] /. sol3][[1]], 8]], s["    dt/dτ"]},
      {"  ", s["Zeitdilatation 4"], " = ", s[N[Evaluate[τ'[Τ4] /. sol4][[1]], 8]], s["    dt/dτ"]},
      {""},
      {"  ", s["Winkel 1"], " = ", s[N[Evaluate[(θ[Τ1] /. sol1) 180/Pi][[1]], 8]], s["    grad"]},
      {"  ", s["Winkel 2"], " = ", s[N[Evaluate[(θ[Τ2] /. sol2) 180/Pi][[1]], 8]], s["    grad"]},
      {"  ", s["Winkel 3"], " = ", s[N[Evaluate[(θ[Τ3] /. sol3) 180/Pi][[1]], 8]], s["    grad"]},
      {"  ", s["Winkel 4"], " = ", s[N[Evaluate[(θ[Τ4] /. sol4) 180/Pi][[1]], 8]], s["    grad"]},
      {""},
      {"  ", s["Radialkoordinate"], " = ", s[N[Evaluate[r[Τ1] /. sol1][[1]], 8]], s["    GM/c²"]},
      {"  ", s["Radialkoordinate 2"], " = ", s[N[Evaluate[r[Τ2] /. sol2][[1]], 8]], s["    GM/c²"]},
      {"  ", s["Radialkoordinate 3"], " = ", s[N[Evaluate[r[Τ3] /. sol3][[1]], 8]], s["    GM/c²"]},
      {"  ", s["Radialkoordinate 4"], " = ", s[N[Evaluate[r[Τ4] /. sol4][[1]], 8]], s["    GM/c²"]},
      {""},
      {"  ", s["x-Achse 1"], " = ", s[N[x1[Τ1], 8]], s["    GM/c²"]},
      {"  ", s["x-Achse 2"], " = ", s[N[x2[Τ2], 8]], s["    GM/c²"]},
      {"  ", s["x-Achse 3"], " = ", s[N[x3[Τ3], 8]], s["    GM/c²"]},
      {"  ", s["x-Achse 4"], " = ", s[N[x4[Τ4], 8]], s["    GM/c²"]},
      {""},
      {"  ", s["y-Achse 1"], " = ", s[N[y1[Τ1], 8]], s["    GM/c²"]},
      {"  ", s["y-Achse 2"], " = ", s[N[y2[Τ2], 8]], s["    GM/c²"]},
      {"  ", s["y-Achse 3"], " = ", s[N[y3[Τ3], 8]], s["    GM/c²"]},
      {"  ", s["y-Achse 4"], " = ", s[N[y4[Τ4], 8]], s["    GM/c²"]},
      {""},
      {"  ", s["v lokal 1"], " = ", s[N[vp1[Τ1], 8]], s["    c"]},
      {"  ", s["v lokal 2"], " = ", s[N[vp2[Τ2], 8]], s["    c"]},
      {"  ", s["v lokal 3"], " = ", s[N[vp3[Τ3], 8]], s["    c"]},
      {"  ", s["v lokal 4"], " = ", s[N[vp4[Τ4], 8]], s["    c"]},
      {""},
      {"  ", s["v extern 1"], " = ", s[N[vc1[Τ1], 8]], s["    c"]},
      {"  ", s["v extern 2"], " = ", s[N[vc2[Τ2], 8]], s["    c"]},
      {"  ", s["v extern 3"], " = ", s[N[vc3[Τ3], 8]], s["    c"]},
      {"  ", s["v extern 4"], " = ", s[N[vc4[Τ4], 8]], s["    c"]}
      
        }, Alignment -> Left, Spacings -> {0, 1/2}]}}, Alignment -> Left]]
        ], {ι, 1, 142, 1}]
      
(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////*)
                                                                                                   (* yukterez.net *)

Freifall vs Photon:

Code:

G = 1; M = 1; c = 1; rs =  2 G M/c^2;
wp = MachinePrecision;
para = 24; pstep =  1/2;

j[v_] := Sqrt[1 - v^2/c^2]; J = j[v0];
k[r_] := Sqrt[1 - rs/r]; κ = k[r0];

r0 = 10;
v0 = 0;
φ = 0;
θ0 = Pi/2;
υ = 35;

vr0 = v0 Sin[φ] κ/J; vθ0 = v0/r0 Cos[φ]/J;

i[x_] := If[r[t] == 2, 1, If[r[t] < 2, Im[x], x]];

sol = NDSolve[
{r''[t] == -((G M)/r[t]^2) + r[t] θ'[t]^2 - (3 G M)/c^2 θ'[t]^2, r'[0] == vr0,
r[0] == r0, θ''[t] == -((2 r'[t] θ'[t])/r[t]),
θ'[0] == vθ0, θ[0] == θ0,
τ'[t] == Sqrt[c^2 r[t] + r[t] r'[t]^2 - c^2 rs + r[t]^3 θ'[t]^2 - r[t]^2 rs θ'[t]^2]/(c Sqrt[r[t] - rs] Sqrt[1 - rs/r[t]]), τ[0] == 0,
cl'[t] == ((r'[t]/k[r[t]])^2 + (θ'[t] r[t])^2)/c^2,
cl[0] == 0}, {r, θ, τ, cl}, {t, 0, υ},
MaxSteps -> Infinity, Method -> Automatic, WorkingPrecision -> wp, InterpolationOrder -> All];

solp = NDSolve[
{r'[t] == -(1 - 2/r[t]), r[0] == 10}, r, {t, 0, 1000},
MaxSteps -> Infinity, Method -> Automatic, WorkingPrecision -> 32, InterpolationOrder -> All];

rp[τ_] := Quiet[Evaluate[r[τ] /. solp][[1]]];
тp = Quiet[tt /. FindRoot[rp[tt] - r1, {tt, 1}, WorkingPrecision -> 32]];

rf[t_] := Quiet[Evaluate[r[t] /. solf][[1]]];
T = Quiet[t /. FindRoot[rf[t] - r1, {t, 1}]];
тf = Quiet[Evaluate[τ[T] /. solf][[1]]];

t[ξ_] := Quiet[
χ /. FindRoot[
Evaluate[τ[χ] /. sol][[1]] - ξ, {χ, 0},
WorkingPrecision -> wp, Method -> Automatic]];
Τ:= Quiet[t[ι]];

u[b_] = b - 2;
Quiet[w[b_] = NIntegrate[Sqrt[1/(1 - 2/R)], {R, 2, b}];
q[b_] = Sqrt[w[b]^2 - u[b]^2];]
grid[n_] = я /. FindRoot[2 + Sqrt[(-2 + я) я] + Log[-1 + я + Sqrt[(-2 + я) я]] == n, {я, 0}];

x[t_] := (Sin[Evaluate[θ[t] /. sol]] Evaluate[r[t] /. sol])[[1]]
y[t_] := (Cos[Evaluate[θ[t] /. sol]] Evaluate[r[t] /. sol])[[1]]

R[t_] := Evaluate[r[t] /. sol][[1]];
γ[t_] := Evaluate[τ'[t] /. sol][[1]];
и[t_] := Evaluate[τ[t] /. sol][[1]];

crθ[t_] := Evaluate[cl'[t] /. sol][[1]]; (*Celerität*)
vrθ[t_] := crθ[t]/Sqrt[1 + crθ[t]^2];
clr[t_] := Evaluate[r'[t] /. sol][[1]];
clθ[t_] := R[t] Evaluate[θ'[t] /. sol][[1]];

vr[t_] := clr[t]/γ[t]/k[R[t]]^2;
vt[t_] := clθ[t]/γ[t]/k[R[t]];
vp[t_] := Sqrt[vr[t]^2 + vt[t]^2];
vc[t_] := Sqrt[vr[t]^2 k[R[t]]^2 + vt[t]^2] k[R[t]];

s[text_] := Style[text, FontSize -> 11]; PR = 12;

Quiet[Do[Print[
Rasterize[Grid[{{Show[
Graphics[{{Gray, Disk[{0, 0}, rs]}, {Black, Dashed,
Circle[{0, 0}, r0]}}, Frame -> True, ImageSize -> 500, PlotRange -> PR, ImagePadding -> Automatic],
Graphics[Table[{Gray, Circle[{0, 0}, grid[n]]}, {n, 2 + pstep, para, pstep}]],
Graphics[{PointSize[0.01], Green, Point[{2, 0}]}],
Graphics[{PointSize[0.01], Blue, Point[{10, 0}]}],
Graphics[{PointSize[0.01], Red, Point[{x[Quiet[t[ι + 30.553877]]], 0}]}],
Graphics[{PointSize[0.01], Green, Point[{rp[ι], 0}]}]]},
{Grid[{
{s["  Zeit A"], " = ", s[N[30.553877 + ι, 8]]},
{s["  Zeit B"], " = ", s[N[(30.553877 + ι) Sqrt[1 - 2/10], 8]]},
{s["  Zeit C"], " = ", s[N[Quiet[t[ι + 30.553877]], 8]]}},
Alignment -> Left, Spacings -> {0, 1/2}]}}, Alignment -> Left]]],
{ι, 0, 40,0.1}]]

Winkeldurchmesser des Schattens:

Code:

(* Angular diameter of a Schwarzschild black hole shadow, schwarzschild.yukterez.net *)

r0 = Infinity;
rs = 2; χ0 = rs/r0; χ = rs/r; p = 4/27;

α[n_] := ArcCos[(χ^2 Sqrt[1-χ0] Sqrt[χ-χ0]+n Sqrt[p(χ^3-χ^2+p)])/(χ^3-χ0 χ^2+p)] 360/π;
δ = If[D[α[-1], r]>0, Evaluate[α[+1]], Evaluate[α[-1]]];

Plot[δ, {r, 0, 10}, Frame->True,
AspectRatio->1/3, ImageSize->500,
GridLines->{Table[n, {n, 1, 10}], {45, 90, 135}},
PlotRange->{{0, 10}, {0, 180}}]

N@Block[{r=10}, δ]
N@Block[{r=0.0000001}, δ]

links: lokal, rechts: extern:



links: lokal (Eigenzeit, Schalengeschwindigkeit), rechts: extern (Koordinatenzeit, verzögerte Geschwindigkeit):

links: Newton, rechts: Einstein (Koordinatenzeit at infinity mit Shapiro-Verzögerung):

Startbedingungen: v0 = neutonische Kreisbahngeschwindigkeit bei r0 = 10 rs = 20 GM/c² und 0° Abschusswinkel:



Startbedingungen: v0 = neutonische Kreisbahngeschwindigkeit bei r0 = 10 rs = 20 GM/c² und 45° Abschusswinkel:

instabile Orbits um die Photonensphäre:



instabile Orbits innerhalb der Photonensphäre:

Probe für MTW 32.1, Animation:



Standbilder: Startbedingung: r0=10, v0=0:



Bei τ=33.7, t=∞ ist der kollabierende Stern auf seinen Schwarzschildradius geschrumpft:



Bei τ=35.1 und t=42.4 ist die Singularität erreicht:



τ,t und τ,r Diagramme:



Benötigte Formeln (G=M=c=1):



Diskussion: Uwudl

Kollaps von 5 Kugelblitzschalen mit jeweils M=1 (Gesamtmasse=5M):



Bewegungsgleichungen:







usw. (siehe Code).

Einfallender Partikelstream, Finkelsteinkoordinaten vs Schwarzschildkoordinaten:



Diskussion: Uwudl, mehr zum Thema: hier entlang.

Zwei mit v=c√(rs/r) ins SL fallende Freifaller (rot) senden bei r=rs/2 ein Signal (grün):



weiterführende Beiträge: SE (Archiv) und Mahag

Startbedingungen: v0 = 1.111fache Kreisbahngeschwindigkeit im Perihel bei r0 = 2.8 rs = 5.6 GM/c²:



Startbedingungen: v0 = 1.02-fache Kreisbahngeschwindigkeit im Perihel bei r0 = 5 GM/c² (für die Kerr-Version hier entlang):



Startbedingungen: v0 = 1.26fache Kreisbahngeschwindigkeit im Perihel bei r0 = 10 rs = 20 GM/c²: