Seite 1 von 1

### Schwarzschild Metric

Verfasst: Do 12. Apr 2018, 22:51

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:

Finkelstein coordinates with the horizon penetrating coordinate time dť→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.