Das ist die deutschsprachige Version. 
English versions are available on en.yukterez.net and yukipedia.Verwandte Beiträge: Kerr-Newman Orbits || Kerr Orbits || Geodäsie || Gravitationslinsen || Raytracing || Flamms Paraboloid
Metrischer Tensor in Schwarzschild Koordinaten {t,r,θ,Ф}:
Alternative Kerr-Schild Koordinaten {u,r,θ,Ф} mit der horizontüberschreitenden Koordinatenzeit u wobei du=dt+2dr/(r-rs):
Bewegungsgleichung in Schwarzschildkoordinaten; aufgrund der Kugelsymmetrie kann auf eine Winkelkoordinate reduziert werden so dass√(dθ²+sin²θ dФ²)→dθ:
mit v⊥=v cos φ als der initialen transversalen, und v∥=v sin φ als der initialen radialen Geschwindigkeit 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 dilatierte 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
sind konstant, wobei
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²:
Alternative Kerr-Schild Koordinaten {u,r,θ,Ф} mit der horizontüberschreitenden Koordinatenzeit u wobei du=dt+2dr/(r-rs):
Bewegungsgleichung in Schwarzschildkoordinaten; aufgrund der Kugelsymmetrie kann auf eine Winkelkoordinate reduziert werden so dass
mit v⊥=v cos φ als der initialen transversalen, und v∥=v sin φ als der initialen radialen Geschwindigkeit 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 dilatierte 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
sind konstant, wobei
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²:
