Nachdem hier immer wieder nach Orbitalen verlangt wird nochmal in aller Deutlichkeit was der Ausdruck "instabiler Orbit" bedeutet. Hier wird ein Partikel mit seiner exakten Orbitalgeschwindigkeit geworfen, allerdings nicht mit exakt 0° sondern mit 1/1000000 (einem Millionstel) Grad ganz leicht nach oben geneigt. Das reicht schon um unseren Wurfgegenstand nach der zweiten Umrundung aus dem System zu katapultieren. Gestartet wird bei r=1.51rs=3.02GM/c²; man achte darauf wie sich die kleinen Änderungen im radialen Abstand immer weiter steigern (Zeile 5):
Demnächst kommen dann die stabilen perihelierenden Orbits bei r>3rs dran. Der upgedatete Code mit ausführlicher Kommentierung und eingebettetem Schwarzschildparaboloid im Hintergrund:
Code:
G = 1; M = 1; c = 1; rs = 2 G M/c^2; (* Einheiten *)
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 *)
T = 2; (* Simulationsdauer *)
vr0 = v0 Sin[φ] κ/J; vθ0 = v0/r0 Cos[φ]/J; (* Längenkontraktion und Tiefenexpansion *)
d1 = T/10; d2 = d1; wp = 30; f = 3; (* Schweifdauer und Frameanzahl *)
para = 20; pstep = 1/2; (* Paraboloid Grid *)
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, T}, WorkingPrecision -> wp,
MaxSteps -> Infinity, Method -> Automatic,
InterpolationOrder -> All];
t[ξ_] := Quiet[ (* Eigenzeit nach Koordinatenzeit *)
χ /.FindRoot[Evaluate[τ[χ] /. sol][[1]] - ξ, {χ, 0}, WorkingPrecision -> wp, Method -> Automatic]];
Τ := Quiet[t[ι]];
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_] = 2 + (n - 2)^2/(w[n, 0] - 2); (* Hintergrundschalen *)
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;
X = 24; Y = 2 X/3; (* Paraboloid Plot *)
Plot[q[x, 0], {x, -X, X},
AspectRatio -> Y/2/X, Frame -> True, PlotRange -> {{-X, X}, {0, Y}}]
Do[Print[ (* Animation nach Eigenzeit *)
Rasterize[Grid[{{Show[Graphics[{
{Black, Circle[{0, 0}, rs]},
{LightGray, Disk[{0, 0}, rs]},
{Lighter[Gray], Dashed, Circle[{0, 0}, r0]}},
Frame -> True, ImageSize -> 400, PlotRange -> 2 r0, ImagePadding -> 1],
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, т},
PlotStyle->{LightGray}],
ParametricPlot[{x[η], y[η]}, {η, 0, т},
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[Evaluate[τ'[т] /. sol][[1]], 8]], s[" dτ/dt"]},
{" ", s["Winkel"], " = ", s[N[Evaluate[(θ[т] /. sol) 180/Pi][[1]], 8]], s[" grad"]},
{" ", s["radialer Abstand"], " = ", s[N[Evaluate[r[т] /. sol][[1]], 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"]}
}, Alignment -> Left]}}, Alignment -> Left]]
], {т, T/f, T, T/f}]
Do[Print[ (* Animation nach Koordinatenzeit *)
Rasterize[Grid[{{Show[Graphics[{
{Black, Circle[{0, 0}, rs]},
{LightGray, Disk[{0, 0}, rs]},
{Lighter[Gray], Dashed, Circle[{0, 0}, r0]}},
Frame -> True, ImageSize -> 400, PlotRange -> 2 r0, ImagePadding -> 1],
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, Τ},
PlotStyle->{LightGray}],
ParametricPlot[{x[η], y[η]}, {η, 0, Τ},
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[Evaluate[τ'[Τ] /. sol][[1]], 8]], s[" dτ/dt"]},
{" ", s["Winkel"], " = ", s[N[Evaluate[(θ[Τ] /. sol) 180/Pi][[1]], 8]], s[" degree"]},
{" ", s["radialer Abstand"], " = ", s[N[Evaluate[r[Τ] /. sol][[1]], 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"]}
}, Alignment -> Left]}}, Alignment -> Left]]
], {ι, и[T]/f, и[T], и[T]/f}]
(* yukterez.net | v4 *)
Invertierend,
Anmelden
Registrieren
FAQ
Suche
, 
hat geschrieben:
, 
, 
, 
, 
