Das ist die deutschsprachige Version. 
English versions are available on en.yukterez.net and yukipedia.Verwandte Beiträge: Kerr Newman || Kerr || Reissner Nordström || Geodäsie || Gravitationslinsen || Raytracing || Flamms Paraboloid
Metrischer Tensor in Schwarzschild Koordinaten {t,r,θ,Ф}:
Erweiterte Finkelsteinkoordinaten {ť,r,θ,Ф} mit der horizontüberschreitenden Koordinatenzeit ť wobei dť=dt+2dr/(r-rs):
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).
τ ist die Eigenzeit des Testpartikels, und t die Koordinatenzeit eines Beobachters at infinty. Um die Schalenzeit eines stationären Beobachters auf r=R zu erhalten wird τ einfach durch √(1-rs/r) dividiert, 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 kinetische Energie im System eines lokalen Schalenbeobachters ist
Die radiale Fluchtgeschwindigkeit um von r0 nach r1 zu gelangen ist
womit die radiale Fluchtgeschwindigkeit in die Unendlichkeit (r1→∞)
ist. Physikalischer Abstand zwischen zwei verschiedenen r:
Bei einer rein radialen Bewegung vereinfacht sich die Bewegungsgleichung zu
Simulator (alte Version in Schwarzschildkoordinaten; für horizontüberschreitende Finkelsteinkoordinaten kann im aktuelleren Kerr-Simulator crd=2 und a=0 gesetzt werden):
Verschiedene Integrationsmethoden (bei Orbits in einem Mindestabstand vom zweifachen EH reicht die Methode Automatic):
Mehrere Testpartikel:
Freifall vs Photon:
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²:
Erweiterte Finkelsteinkoordinaten {ť,r,θ,Ф} mit der horizontüberschreitenden Koordinatenzeit ť wobei dť=dt+2dr/(r-rs):
Bewegungsgleichung in Schwarzschildkoordinaten; aufgrund der Kugelsymmetrie kann auf eine Winkelkoordinate reduziert werden so dass
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).
τ ist die Eigenzeit des Testpartikels, und t die Koordinatenzeit eines Beobachters at infinty. Um die Schalenzeit eines stationären Beobachters auf r=R zu erhalten wird τ einfach durch √(1-rs/r) dividiert, 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 kinetische Energie im System eines lokalen Schalenbeobachters ist
Die radiale Fluchtgeschwindigkeit um von r0 nach r1 zu gelangen ist
womit die radiale Fluchtgeschwindigkeit in die Unendlichkeit (r1→∞)
ist. Physikalischer Abstand zwischen zwei verschiedenen r:
Bei einer rein radialen Bewegung vereinfacht sich die Bewegungsgleichung zu
Simulator (alte Version in Schwarzschildkoordinaten; für horizontüberschreitende Finkelsteinkoordinaten kann im aktuelleren Kerr-Simulator crd=2 und a=0 gesetzt werden):
Code: Alles auswählen
(*||| Schwarzschild Simulator |||||||| yukterez.net ||||||||||| Version 13.02.2018 |||||||| Syntax: Mathematica |||*)
(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
ClearAll["Global`*"] (* Variablen frei machen *)
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: Alles auswählen
Method -> {"EventLocator", "Event" -> (r[t] - 2.000001)}Code: Alles auswählen
Method -> {"StiffnessSwitching", Method -> {"ExplicitRungeKutta", Automatic}}Code: Alles auswählen
Method -> {"ImplicitRungeKutta", DifferenceOrder" -> 20}Mehrere Testpartikel:
Code: Alles auswählen
(* 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: Alles auswählen
ClearAll["Global`*"]; \
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]; \[Kappa] = k[r0];
r0 = 10;
v0 = 0;
\[CurlyPhi] = 0;
\[Theta]0 = Pi/2;
\[Upsilon] = 35;
vr0 = v0 Sin[\[CurlyPhi]] \[Kappa]/J; v\[Theta]0 = v0/r0 Cos[\[CurlyPhi]]/J;
Needs["DifferentialEquations`NDSolveProblems`"]; \
Needs["DifferentialEquations`NDSolveUtilities`"];
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] \[Theta]'[t]^2 - (3 G M)/c^2 \[Theta]'[t]^2, r'[0] == vr0,
r[0] == r0, \[Theta]''[t] == -((2 r'[t] \[Theta]'[t])/r[t]),
\[Theta]'[0] == v\[Theta]0, \[Theta][0] == \[Theta]0,
\[Tau]'[t] == Sqrt[c^2 r[t] + r[t] r'[t]^2 - c^2 rs + r[t]^3 \[Theta]'[t]^2 - r[t]^2 rs \[Theta]'[t]^2]/(c Sqrt[r[t] - rs] Sqrt[1 - rs/r[t]]), \[Tau][0] == 0,
cl'[t] == ((r'[t]/k[r[t]])^2 + (\[Theta]'[t] r[t])^2)/c^2,
cl[0] == 0}, {r, \[Theta], \[Tau], cl}, {t, 0, \[Upsilon]},
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[\[Tau]_] := Quiet[Evaluate[r[\[Tau]] /. 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[\[Tau][T] /. solf][[1]]];
t[\[Xi]_] := Quiet[
\[Chi] /. FindRoot[
Evaluate[\[Tau][\[Chi]] /. sol][[1]] - \[Xi], {\[Chi], 0},
WorkingPrecision -> wp, Method -> Automatic]];
\[CapitalTau] := Quiet[t[\[Iota]]];
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[\[Theta][t] /. sol]] Evaluate[r[t] /. sol])[[1]]
y[t_] := (Cos[Evaluate[\[Theta][t] /. sol]] Evaluate[r[t] /. sol])[[1]]
R[t_] := Evaluate[r[t] /. sol][[1]];
\[Gamma][t_] := Evaluate[\[Tau]'[t] /. sol][[1]];
и[t_] := Evaluate[\[Tau][t] /. sol][[1]];
cr\[Theta][t_] := Evaluate[cl'[t] /. sol][[1]]; (*Celerität*)
vr\[Theta][t_] := cr\[Theta][t]/Sqrt[1 + cr\[Theta][t]^2];
clr[t_] := Evaluate[r'[t] /. sol][[1]];
cl\[Theta][t_] := R[t] Evaluate[\[Theta]'[t] /. sol][[1]];
vr[t_] := clr[t]/\[Gamma][t]/k[R[t]]^2;
vt[t_] := cl\[Theta][t]/\[Gamma][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[\[Iota] + 30.553877]]], 0}]}],
Graphics[{PointSize[0.01], Green, Point[{rp[\[Iota]], 0}]}]]},
{Grid[{
{s[" Zeit A"], " = ", s[N[30.553877 + \[Iota], 8]]},
{s[" Zeit B"], " = ", s[N[(30.553877 + \[Iota]) Sqrt[1 - 2/10], 8]]},
{s[" Zeit C"], " = ", s[N[Quiet[t[\[Iota] + 30.553877]], 8]]}},
Alignment -> Left, Spacings -> {0, 1/2}]}}, Alignment -> Left]]],
{\[Iota], 0, 40,0.1}]]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²:
