This is the english version. 
Deutschsprachige Version auf schwarzschild.yukterez.net.
Shadow and horizon of a black hole, the observer is at a distance of r=50GM/c². Zoom out: [-], ray animation: play, raytracing code: click
left: lightrays in flat Minkoswki space, right: Schwarzschild. Distance of the light source to the black hole: 20GM/c²
Orbit with perihelion shift; initial conditions: r0=5, θ0=π/2, v0=vz0=vθ0=51/50·√((1/5)/(1-2/5)). Geodesic solver: geodesics.txt
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):
Gullstrand Painlevé (Raindrop) coordinates, horizon penetrating coordinate time defined by freefallers from infinity:
Equations of motion in Schwarzschild coordinates:
For a purely radial motion the equation of motion simplifies to
τ is 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. v⊥=v cos ζ (the transverse), and v∥=v sin ζ (the radial component of the local velocity). ζ 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 in the frame of an observer at infinity
are conserved. The rest, kinetic and potential energy (defined as the difference between local and total energy) are
The required radial velocity to get from r0 to r1 is
and the escape velocity from r0 to infinity
The physical distance between r1 and r2 is
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):
Gullstrand Painlevé (Raindrop) coordinates, horizon penetrating coordinate time defined by freefallers from infinity:
Equations of motion in Schwarzschild coordinates:
For a purely radial motion the equation of motion simplifies to
τ is 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. v⊥=v cos ζ (the transverse), and v∥=v sin ζ (the radial component of the local velocity). ζ 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 in the frame of an observer at infinity
are conserved. The rest, kinetic and potential energy (defined as the difference between local and total energy) are
The required radial velocity to get from r0 to r1 is
and the escape velocity from r0 to infinity
The physical distance between r1 and r2 is
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.
