This is the english version. 
Deutschsprachige Version auf schwarzschild.yukterez.net.Metric tensor in Schwarzschild coordinates {t,r,θ,φ}; superscripted letters are not powers but indices:
Alternative Kerr-Schild form {u,r,θ,φ} with the horizon penetrating finkelsteinlike null coordinate time u, where du=dt+2dr/(r-rs):
Equations of motion in Schwarzschild coordinates:
τ ist the proper time of the test particle, and t the coordinate time of an observer at infinty. To get the shell time of a stationary fiducional observer at r=R, τ gets divided by √(1-rs/r), where rs=2GM/c² is the Schwarzschildradius. The total time dilation is the product of the gravitational and the kinetic component.
The transverse and radial components of the local velocity are
with v⊥=v cos φ as the transverse, and v∥=v sin φ as the radial component. φ is the vertical launch angle (because of the radial length contraction φ at small r looks flatter when viewed from infinity).
Transformation between local and observed (shapirodelayed) velocities:
With Pythagoras we get the total velocity:
The orbital angular momentum
and the total energy of the test particle
are conserved, where
Codes in Mathematica Syntax:
Different integration methods (for orbits which don't scratch the horizon Automatic is just fine):
Multiple test particles:
Free faller versus photon:
The comments in the code are in german, but you can translate them with Google Translate or any other translator if anything is unclear. For images and animations see the german version of this site.
Alternative Kerr-Schild form {u,r,θ,φ} with the horizon penetrating finkelsteinlike null coordinate time u, where du=dt+2dr/(r-rs):
Equations of motion in Schwarzschild coordinates:
τ ist the proper time of the test particle, and t the coordinate time of an observer at infinty. To get the shell time of a stationary fiducional observer at r=R, τ gets divided by √(1-rs/r), where rs=2GM/c² is the Schwarzschildradius. The total time dilation is the product of the gravitational and the kinetic component.
The transverse and radial components of the local velocity are
with v⊥=v cos φ as the transverse, and v∥=v sin φ as the radial component. φ is the vertical launch angle (because of the radial length contraction φ at small r looks flatter when viewed from infinity).
Transformation between local and observed (shapirodelayed) velocities:
With Pythagoras we get the total velocity:
The orbital angular momentum
and the total energy of the test particle
are conserved, where
Codes in Mathematica Syntax:
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 *)
Different integration methods (for orbits which don't scratch the horizon Automatic is just fine):
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}Multiple test particles:
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 *)Free faller versus 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}]]The comments in the code are in german, but you can translate them with Google Translate or any other translator if anything is unclear. For images and animations see the german version of this site.