Seite 1 von 1

Schwarzschild Orbits

Verfasst: Sa 21. Mai 2016, 23:08
von Yukterez
Update: ENGLISH VERSION Bild Bild
Bild

Verwandte Beiträge: Kerr-Newman Orbits || Kerr Orbits || Geodätengleichung || Gravitationslinseneffekt
Bild

Die Differentialgleichung lautet





ist die Eigenzeit des Testpartikels, und die Koordinatenzeit eines Beobachters at infinty. Um die Koordinatenzeit eines stationären Beobachters im Abstand vom Schwerpunkt zu erhalten wird einfach durch dividiert, mit dem Schwarzschildradius . Die totale Zeitdilatation ist das multiplikative Produkt aus der gravitativen und der kinematischen Zeitdilatation.

Die transversale und radiale Komponente der lokalen Geschwindigkeit lautet





mit als der initialen transversalen, und als der initialen radialen Geschwindigkeit mit dem lokalen Abschusswinkel (von außen erscheint dieser aufgrund der radialen Längenkontraktion flacher).

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 die einzelnen Komponenten der Energie



sind.

Code: Alles auswählen

(* Schwarzschild Simulator, relativistische Wurfparabel ||| yukterez.net  Version 06, 2017 ||| Syntax: Mathematica *)
(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

ClearAll["Global`*"]                                       (* Variablen frei machen *)
G = 1; M = 1; c = 1; rs = 2 G M/c^2;                       (* Einheiten *)
wp = MachinePrecision;                                     (* numerische Genauigkeit *)
para = 20; pstep = 1/2;                                    (* Paraboloid Grid *)

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

j[v_] := Sqrt[1 - v^2/c^2];    J = j[v0];                  (* Gammafaktor *)
k[r_] := Sqrt[1 - rs/r];       κ = k[r0];                  (* Schwarzschildfaktor *)

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

r0 = 151/100 rs;                                           (* Startradius *)
v0 = Sqrt[G M/r0]/κ;                                       (* lokale Startgeschwindigkeit *)
φ  = Pi/4;                                                 (* Abschusswinkel *)
θ0 = 0;                                                    (* Startwinkel *)
υ  = 2;                                                    (* Simulationsdauer *)

(* г = Sqrt[χ^2 + γ^2]; Θ = ArcSin[γ/г]; *)                (* Kartesisch auf Polar *)

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

vr0 = v0 Sin[φ] κ/J; vθ0 = v0/r0 Cos[φ]/J;                 (* Längenkontraktion und Tiefenexpansion *)
d1 = υ/10; d2 = d1; f = 3;                                 (* Schweifdauer und Frameanzahl *)

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

Needs["DifferentialEquations`NDSolveProblems`"]; Needs["DifferentialEquations`NDSolveUtilities`"];
i[x_] := If[r[t] == 2, 1, If[r[t] < 2, Im[x], x]];

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, υ},
MaxSteps -> Infinity,
Method -> Automatic,
WorkingPrecision -> wp,
InterpolationOrder -> All];

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

t[ξ_] := Quiet[                                            (* Eigenzeit nach Koordinatenzeit *)
χ /.FindRoot[Evaluate[τ[χ] /. sol][[1]] - ξ, {χ, 0}, WorkingPrecision -> wp, Method -> Automatic]];
Τ := Quiet[t[ι]];

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

u[b_] := b - 2;
w[b_] := Quiet[NIntegrate[Sqrt[1/(1 - 2/R)], {R, 2 + 1*^-8, b}]];
q[b_] := Quiet[Sqrt[w[b]^2 - u[b]^2]];
grid[n_] := Quiet[Re[я /. FindRoot[2 + Quiet[NIntegrate[Sqrt[1/(1 - 2/R)], {R, 2, я}]] == n, {я, 3}]]];

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

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]];                       (* radialer Abstand *)
γ[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 *)

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

Xa = 24; Ya = 2 Xa/3;                                      (*Paraboloid Plot*)
Plot[{q[x], q[-x]}, {x, -Xa, Xa}, AspectRatio -> Ya/2/Xa,
Frame -> True, PlotRange -> {{-Xa, Xa}, {0, Ya}}, PlotStyle -> Black]

(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

Do[Print[                                                  (* Animation nach Eigenzeit *)
Rasterize[Grid[{{Show[Graphics[{
{LightGray, Disk[{0, 0}, rs]},
{Lighter[Gray], Dashed, Circle[{0, 0}, r0]}},
Frame -> True, ImageSize -> 360,
PlotRange -> PR, ImagePadding ->  Automatic],
Graphics[Table[{LightGray, Circle[{0, 0}, grid[n]]}, {n, 2 + pstep, para, pstep}]],
Graphics[{PointSize[0.01], Red, Point[{x[т], y[т]}]}],

ParametricPlot[{x[η], y[η]}, {η, 0, Max[1*^-16, т - d1]},
PlotStyle->{Thickness[0.004], Lighter[Lighter[Gray]]}],
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[" dτ/dt"]},
      {s["  Winkel"],              " = ",    s[N[Evaluate[(θ[т] /. sol) 180/Pi][[1]], 8]],    s[" grad"]},
      {s["  radialer Abstand"],    " = ",    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]]
      ], {т, υ/f, υ, υ/f}]
   
(*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)

Do[Print[                                                  (* Animation nach Koordinatenzeit *)
Rasterize[Grid[{{Show[Graphics[{
{LightGray, Disk[{0, 0}, rs]},
{Lighter[Gray], Dashed, Circle[{0, 0}, r0]}},
Frame -> True, ImageSize -> 360,
PlotRange -> PR, ImagePadding -> Automatic],
Graphics[Table[{LightGray, Circle[{0, 0}, grid[n]]}, {n, 2 + pstep, para, pstep}]],
Graphics[{PointSize[0.01], Red, Point[{x[Τ], y[Τ]}]}],

ParametricPlot[{x[η], y[η]}, {η, 0, Max[1*^-16, Τ - d1]},
PlotStyle->{Thickness[0.004], Lighter[Lighter[Gray]]}],
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[" dτ/dt"]},
      {s["  Winkel"],              " = ",    s[N[Evaluate[(θ[Τ] /. sol) 180/Pi][[1]], 8]],    s[" grad"]},
      {s["  radialer Abstand"],    " = ",    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]]
      ], {ι, и[υ]/f, и[υ], и[υ]/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]];                     (* radialer Abstand *)
γ1[t_] := Evaluate[τ'[t] /. sol1][[1]];                    (* Zeitdilatation *)
и1[t_] := Evaluate[τ[t] /. sol1][[1]];                     (* Koordinatenzeit *)

R2[t_] := Evaluate[r[t] /. sol2][[1]];                     (* radialer Abstand *)
γ2[t_] := Evaluate[τ'[t] /. sol2][[1]];                    (* Zeitdilatation *)
и2[t_] := Evaluate[τ[t] /. sol2][[1]];                     (* Koordinatenzeit *)

R3[t_] := Evaluate[r[t] /. sol3][[1]];                     (* radialer Abstand *)
γ3[t_] := Evaluate[τ'[t] /. sol3][[1]];                    (* Zeitdilatation *)
и3[t_] := Evaluate[τ[t] /. sol3][[1]];                     (* Koordinatenzeit *)

R4[t_] := Evaluate[r[t] /. sol4][[1]];                     (* radialer Abstand *)
γ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["    dτ/dt"]},
      {"  ", s["Zeitdilatation 2"], " = ", s[N[Evaluate[τ'[Τ2] /. sol2][[1]], 8]], s["    dτ/dt"]},
      {"  ", s["Zeitdilatation 3"], " = ", s[N[Evaluate[τ'[Τ3] /. sol3][[1]], 8]], s["    dτ/dt"]},
      {"  ", s["Zeitdilatation 4"], " = ", s[N[Evaluate[τ'[Τ4] /. sol4][[1]], 8]], s["    dτ/dt"]},
      {""},
      {"  ", 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["radialer Abstand 1"], " = ", s[N[Evaluate[r[Τ1] /. sol1][[1]], 8]], s["    GM/c²"]},
      {"  ", s["radialer Abstand 2"], " = ", s[N[Evaluate[r[Τ2] /. sol2][[1]], 8]], s["    GM/c²"]},
      {"  ", s["radialer Abstand 3"], " = ", s[N[Evaluate[r[Τ3] /. sol3][[1]], 8]], s["    GM/c²"]},
      {"  ", s["radialer Abstand 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²:

Bild

Startbedingungen: v0 = 1.02-fache Kreisbahngeschwindigkeit im Perihel bei r0 = 5 GM/c² (für die Kerr-Version hier entlang):

Bild

Startbedingungen: v0 = 1.26fache Kreisbahngeschwindigkeit im Perihel bei r0 = 10 rs = 20 GM/c²:

Bild

lokale und dilatierte Geschwindigkeit

Verfasst: Sa 21. Mai 2016, 23:09
von Yukterez
links: lokal, rechts: extern:

Bild

links: lokal (Eigenzeit, Schalengeschwindigkeit), rechts: extern (Koordinatenzeit, verzögerte Geschwindigkeit):

Bild

Newton vs Einstein

Verfasst: Sa 21. Mai 2016, 23:09
von Yukterez
links: Newton, rechts: Einstein (Koordinatenzeit at infinity mit Shapiro-Verzögerung):

Bild

Newton vs Einstein

Verfasst: So 22. Mai 2016, 10:52
von Yukterez
Startbedingungen: v0 = neutonische Kreisbahngeschwindigkeit bei r0 = 10 rs = 20 GM/c² und 0° Abschusswinkel:

Bild

Startbedingungen: v0 = neutonische Kreisbahngeschwindigkeit bei r0 = 10 rs = 20 GM/c² und 45° Abschusswinkel:

Bild

innere Photonensphäre

Verfasst: Do 26. Mai 2016, 09:29
von Yukterez
instabile Orbits um die Photonensphäre:

Bild

instabile Orbits innerhalb der Photonensphäre:

Bild

Schwarzschild Paraboloid

Verfasst: Sa 4. Jun 2016, 08:27
von Yukterez
Die radiale Strecke zwischen den Schwarzschildkoordinaten r1=2 und r2=24 beträgt lokal nicht r2-r1=22, sondern √(528)+ln(23+√(528)) = 26.8 GM/c²:

Bild
x=r, y=√{R²-r²}


Auf der x-Achse sind die Koordinatenabstände (r=U/2/π) in Einheiten von GM/c² angegeben, während in der Tat so viele stationäre Lineale dazwischen passen wie sie auf dem Paraboloid Platz haben. Die Koordinaten bis 3GM/c² in der Draufsicht:

Bild
x=x, y=y.


unter Newton und im euklidschen Raum sähe das entsprechende Raster so aus (U=2πr):

Bild

Die x- und y-Achsen sind in beiden Bildern von (-3..+3)GM/c² abgebildet, mit einem lokalen Schalenabstand von GM/c²/5..

Anhang

Verfasst: Mo 6. Jun 2016, 12:32
von Yukterez