Im Folgenden wird die Geodätengleichung in vorgegebener Schwarzschild Raumzeit betrachtet. Die Geodätengleichung beschreibt wie sich ein Probekörper (Masse Probekörper $m << M$ Masse schwarzes Loch) im Raum bewegt und sagt voraus, dass diese Bewegung sich stets entlang der extremalsten Kurve, in der durch die Metrik beschriebenen gekrümmten Raumzeit, vollzieht.
Die Geodätengleichung lässt sich durch folgendes Variationsprinzip herleiten
\begin{equation} \int_{x^\mu(A)}^{x^\mu(B)} ds= \int_{x^\mu(A)}^{x^\mu(B)} \sqrt{g_{\mu\nu} dx^\mu dx^\nu}=\int_A^B \sqrt{g_{\mu\nu} \frac{dx^\mu}{d\lambda} \frac{dx^\nu}{d\lambda}} d\lambda \rightarrow \hbox{Extremal} \end{equation}
, wobei man, bei festgehaltenen Anfangs- und Endpunkten ($x^\mu(A)$ und $x^\mu(B)$), nach den extremalsten Kurven $x^\mu(\lambda)$ in der Raumzeit sucht. Den Parameter $\lambda$ bezeichnet man als affiner Parameter, welcher z.B. die Eigenzeit $\tau$ des Körpers sein kann. Die Geodätengleichung ergibt sich dann mittels der Euler-Lagrange Gleichungen, wobei man als Lagrangedichte den Ausdruck $L = \sqrt{g_{\mu\nu} \frac{dx^\mu}{d\lambda} \frac{dx^\nu}{d\lambda}}$, bzw. alternativ $L = g_{\mu\nu} \frac{dx^\mu}{d\lambda} \frac{dx^\nu}{d\lambda}$ verwendet.
\begin{equation} \frac{d}{d\lambda} \left( \frac{\partial L}{\partial \frac{d x^\alpha}{d \lambda}} \right) - \frac{\partial L}{\partial x^\alpha} = 0 \quad \xrightarrow{\cdot \frac{1}{2}g^{\mu \alpha }} \quad \frac{d^2 x^\mu}{d\lambda^2} + \Gamma^\mu_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda} ~=~ 0 \end{equation}
$\Gamma^\mu_{\nu \rho}$ sind die Christoffel Symbole zweiter Art und $\lambda$ ein affiner Parameter (z.B. die Eigenzeit $\tau$ des Probekörpers). Die Geodätengleichung ist ein System gekoppelter Differenzialgleichungen, das wir uns im Folgenden näher betrachten werden.
Wir betrachten im Folgenden die Geodätengleichung in vorgegebener Schwarzschild-Raumzeit
$$ g_{\mu\nu}=\left( \begin{array}{ccc} 1-\frac{2\,M}{r} & 0 & 0 & 0\\ 0& -\frac{1}{1-\frac{2\,M}{r}}& 0&0 \\ 0& 0& -r^2& 0\\ 0& 0& 0& -r^2 \hbox{sin}^2(\theta)\\ \end{array} \right) $$
, wobei wir wieder ein sphärisches Koordinatensystem benutzen ($x^\mu=\left(t,r,\theta,\phi \right)$). Die Geodätengleichung stellt ein System gekoppelter nichtlinearer Differenzialgleichungen dar
$$ \begin{eqnarray} && \frac{d^2 t}{d\lambda^2} = - \Gamma^0_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda} \\ && \frac{d^2 r}{d\lambda^2} = - \Gamma^1_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda}\\ && \frac{d^2 \theta}{d\lambda^2} = - \Gamma^2_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda}\\ && \frac{d^2 \phi}{d\lambda^2} = - \Gamma^3_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda} \quad , \end{eqnarray} $$
wobei $\lambda$ ein affiner Parameter (z.B. die Eigenzeit $\tau$), t, r, $\theta$ und $\phi$ die Koordinaten und $\Gamma^\mu_{\nu \rho}$ die Christoffel Symbole zweiter Art darstellen.
Zunächst wird das Python-Modul "EinsteinPy" eingebunden, welches auf dem Modul SymPy basiert und symbolische Berechnungen in der Allgemeinen Relativitätstheorie relativ einfach möglich macht. Dann definieren wir die zugrundeliegenden Koordinaten und die kovariante Form der Raumzeit-Metrik der Schwarzschildmetrik.
from sympy import *
init_printing()
from einsteinpy.symbolic import *
t, r, theta, phi, M = symbols('t, r, theta, phi, M')
Metric = diag((1-2*M/r), -1/(1-2*M/r), -r**2, -r**2*sin(theta)**2).tolist()
g = MetricTensor(Metric, [t, r, theta, phi])
chr = ChristoffelSymbols.from_metric(g)
Geodätische Bewegung eines Probekörpers¶
Um die Geodätengleichung
$$ \frac{d^2 x^\mu}{d\tau^2} + \Gamma^\mu_{\nu \rho} \frac{d x^\nu}{d\tau} \frac{d x^\rho}{d\tau} ~=~ 0 $$
zu berechnen, definieren zunächst den ersten Term der Gleichung $\frac{d^2 x^\mu}{d\tau^2}$:
tau = symbols('tau')
x0 = Function('t')(tau)
x1 = Function('r')(tau)
x2 = Function('theta')(tau)
x3 = Function('phi')(tau)
xmu = GenericVector([x0, x1, x2, x3],[t, r, theta, phi], config='u',parent_metric=Metric)
Geod1 = diff(xmu.tensor(),tau,tau)
und nun den zweiten Term $\Gamma^\mu_{\nu \rho} \frac{d x^\nu}{d\tau} \frac{d x^\rho}{d\tau}$:
Geod2T = tensorproduct(chr.tensor(),diff(xmu.tensor(),tau),diff(xmu.tensor(),tau))
Geod2 = tensorcontraction(Geod2T, (1, 3),(2,4))
Die Geodätengleichung ist ein System von nichtlinearen gekoppelten Differenzialgleichungen zweiter Ordnung. Wir definieren die entsprechenden vier Gleichungen:
Geod = Geod1 + Geod2
GeodEq0 = Eq(Geod[0],0)
GeodEq1 = Eq(Geod[1],0)
GeodEq2 = Eq(Geod[2],0)
GeodEq3 = Eq(Geod[3],0)
GeodEq0
GeodEq1
GeodEq2
GeodEq3
Beispiel einer gebundenen elliptischen Bahn¶
Wir lassen im Folgenden nur ebene Bewegungen zu ($\theta=\frac{\pi}{2}, {\rm cos}(\theta)=0, {\rm sin}(\theta)=1, \frac{d \theta}{d\tau}=0, \frac{d^2 \theta}{d\tau^2}=0$). Unter diesen Annahmen sind nur drei der vier Geodätengleichungen relevant und diese lauten dann:
GeodEq0 = GeodEq0.subs({(theta,pi/2),(diff(x2,tau),0)})
GeodEq1 = GeodEq1.subs({(theta,pi/2),(diff(x2,tau),0)})
GeodEq3 = GeodEq3.subs({(theta,pi/2),(diff(x2,tau),0)})
GeodEq0_a = Eq(diff(x0,tau,tau),solve(GeodEq0,diff(x0,tau,tau))[0])
GeodEq1_a = Eq(diff(x1,tau,tau),solve(GeodEq1,diff(x1,tau,tau))[0])
GeodEq3_a = Eq(diff(x3,tau,tau),solve(GeodEq3,diff(x3,tau,tau))[0])
GeodEq0_a
GeodEq1_a
GeodEq3_a
Um dieses System von drei gekoppelten Differenzialgleichungen zweiter Ordnung mit Python nummerisch lösen zu können, müssen wir es zunächst in ein System von sechs Differenzialgleichungen erster Ordnung umschreiben. Details zur numerischen Lösung von Differenzialgleichungen finden Sie in der Vorlesung Einführung in die Programmierung für Studierende der Physik in den Unterpunkten Differenzialgleichungen: Numerische Lösung von Anfangswertproblemen und Systeme von gekoppelten Differenzialgleichungen und Differenzialgleichungen zweiter Ordnung.
Wir definieren dazu formal
$$ t(\tau) = y_1(\tau) \, , \quad r(\tau) = y_2(\tau), \quad \phi(\tau) = y_3(\tau) \, , \quad \frac{d t(\tau)}{d\tau} = y_4(\tau) \, , \quad \frac{d r(\tau)}{d\tau} = y_5(\tau), \quad \frac{d \phi(\tau)}{d\tau} = y_6 $$
und schreiben unser DGL-System wie folgt um:
y1, y2, y3, y4, y5, y6 = symbols('y_1, y_2, y_3, y_4, y_5, y_6')
SysDGL_0 = Eq(diff(x0,tau,tau),GeodEq0_a.rhs.subs({(r,y2),(diff(x0,tau),y4),(diff(x1,tau),y5),(diff(x3,tau),y6)}))
SysDGL_0
SysDGL_1 = Eq(diff(x1,tau,tau),GeodEq1_a.rhs.subs({(r,y2),(diff(x0,tau),y4),(diff(x1,tau),y5),(diff(x3,tau),y6)}))
SysDGL_1
SysDGL_3 = Eq(diff(x3,tau,tau),GeodEq3_a.rhs.subs({(r,y2),(diff(x0,tau),y4),(diff(x1,tau),y5),(diff(x3,tau),y6)}))
SysDGL_3
Dieses System von drei gekoppelten DGLs 2. Ordnung in der Eigenzeit $\tau$ werden wir später zusätzlich mittels der Programmiersprache C++ lösen. Wir exportieren deshalb schon jetzt das zugrundeliegende System von DGLs als C++ Code:
from sympy.codegen.ast import Assignment
wert = symbols('wert')
print(ccode(Assignment(wert,SysDGL_0.rhs)))
print(ccode(Assignment(wert,SysDGL_1.rhs)))
print(ccode(Assignment(wert,SysDGL_3.rhs)))
wert = 2*M*y_4*y_5/(y_2*(2*M - y_2)); wert = -2*pow(M, 2)*pow(y_5, 2)/(y_2*pow(2*M - y_2, 2)) + 2*pow(M, 2)*pow(y_4, 2)/pow(y_2, 3) + M*pow(y_5, 2)/pow(2*M - y_2, 2) - M*pow(y_4, 2)/pow(y_2, 2) + pow(y_6, 2)*pow(-2*M + y_2, 3)/pow(2*M - y_2, 2); wert = -2*y_5*y_6/y_2;
Nun lösen wir die DGLs mittels Python und setzen im Folgenden $M=1$. Wir benutzen hierfür das Python-Modul SciPy, welches eine breite Kollektion von mathematischen Algorithmen und Funktionen bereitstellt. Im Speziellen werden wir in diesem Semester die Funktion 'solve_ivp(...)' verwenden, die sich im Untermodul 'scipy.integrate' befindet, welches Funktionen zum Lösen von gewöhnlichen Differenzialgleichungen bereitstellt. Zusätzlich berechnen wir die numerische Lösung auch mittels der älteren Methode "integrate.odeint()", die ebenfalls in dem Python-Modul "scipy" definiert ist.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
from scipy import integrate
Wir definieren das System der gekoppelten sechs Differenzialgleichungen:
set_M = 1
f_4 = lambdify((y1,y2,y3,y4,y5,y6), SysDGL_0.rhs.subs(M,set_M))
f_5 = lambdify((y1,y2,y3,y4,y5,y6), SysDGL_1.rhs.subs(M,set_M))
f_6 = lambdify((y1,y2,y3,y4,y5,y6), SysDGL_3.rhs.subs(M,set_M))
def DGLsys(vy, time):
y1,y2,y3,y4,y5,y6 = vy
dy1 = y4
dy2 = y5
dy3 = y6
dy4 = f_4(y1,y2,y3,y4,y5,y6)
dy5 = f_5(y1,y2,y3,y4,y5,y6)
dy6 = f_6(y1,y2,y3,y4,y5,y6)
return np.array([dy1,dy2,dy3,dy4,dy5,dy6])
def DGLsys_ivp(time, vy):
y1,y2,y3,y4,y5,y6 = vy
dy1 = y4
dy2 = y5
dy3 = y6
dy4 = f_4(y1,y2,y3,y4,y5,y6)
dy5 = f_5(y1,y2,y3,y4,y5,y6)
dy6 = f_6(y1,y2,y3,y4,y5,y6)
return np.array([dy1,dy2,dy3,dy4,dy5,dy6])
In Abhängigkeit von den Anfangswerten können unterschiedliche Bahnen der Bewegung entstehen. Wir wählen zunächst als Beispiel einer gebundenen elliptischen Bahn. Zur Zeit $t=\tau=0$ befindet sich der Probekörper bei einem Radius von $r_0=10=5 R_S$ und Winkel $\phi_0=0$. Seine Anfangsgeschwindigkeit in radialer Richtung sei $\left. \frac{dr}{d\tau} \right|_{\tau=0}=0$ und seine Anfangsgeschwindigkeit in $\phi$-Richtung sei $\left. \frac{dr}{d\tau} \right|_{\tau=0}=0.041$. Die weitere nötige Anfangsbedingung ( $\left. \frac{dt}{d\tau} \right|_{\tau=0}$ ) erhalten wir mittels der Bedingung $\left( \frac{ds}{d\tau} \right)^2= g_{\mu\nu} \frac{dx^\mu}{d\tau} \frac{dx^\nu}{d\tau} = g_{\mu\nu} u^\mu u^\nu = 1$, die für massive Probekörper immer erfüllt sein muss ($u^\mu$ ist die Vierergeschwindigkeit des Probekörpers). Die Gleichung für das infinitesimale Weglängenelement $ds^2=1$ in der Ebene lautet:
dt, dr, dtheta, dphi = symbols('dt, dr, d\\theta, d\\phi')
dx = GenericVector([dt, dr, dtheta, dphi],[t, r, theta, phi], config='u',parent_metric=Metric)
dxdx = tensorproduct(dx.tensor(),dx.tensor())
dxdxT = Tensor(dxdx,config="uu")
ds2a = tensorproduct(g.tensor(),dxdxT.tensor())
ds2aT = Tensor(ds2a,config="lluu")
ds2 = tensorcontraction(ds2aT.tensor(), (0, 2),(1, 3))
Eq_ds = Eq(ds2.subs({(theta,pi/2),(dtheta,0)}),1)
Eq_ds
Wir lösen diese Gleichung nach $dt$ auf:
LoesEq_ds = solve(Eq_ds,dt)
LoesEq_ds
Wir benutzen die erste Lösung (positive Lösung) und berechnen für unsere Anfangsbedingung $\left. \frac{dt}{d\tau} \right|_{\tau=0}=1.20835632162041$.
N = 100000
tauval = np.linspace(0, 800, N+1)
r0 = 10.0
t0 = 0.0
phi0 = 0.0
dr0 = 0.0
dphi0 = 0.041
dt0 = float(LoesEq_ds[0].subs({(r,r0),(dr,dr0),(dphi,dphi0),(M,set_M)}))
initialval = np.array([t0,r0,phi0,dt0,dr0,dphi0])
Loes = integrate.odeint(DGLsys, initialval, tauval)
Loes_ivp = integrate.solve_ivp(DGLsys_ivp, [tauval[0], tauval[-1]], initialval, t_eval=tauval, rtol=10**(-10), atol=10**(-10))
print("dt0 =",dt0)
dt0 = 1.208356321620407
rcParams.update({
'figure.figsize' : [6,6],
'axes.titlesize' : 14,
'axes.labelsize' : 16,
'xtick.labelsize' : 14 ,
'ytick.labelsize' : 14
})
plt.cla()
plt.ylabel(r"$\rm Radius \,\,r \,[km]$")
plt.xlabel(r"$\rm Eigenzeit \,\,\tau$")
plt.plot(tauval,Loes[:, 1],c="blue", linewidth=1.5, linestyle='-');
plt.plot(Loes_ivp.t,Loes_ivp.y[1],c="red", linewidth=2, linestyle=':');
Das obere Diagramm zeigt den Abstand des Probekörpers vom schwarzen Loch (Radius r) als Funktion der Eigenzeit $\tau$, wobei die blaue Kurve die Ergebnisse der odeint-Methode und die gepunktete rote Kurve die solve_ivp-Simulation zeigt. Beide Simulationen stimmen gut miteinander überein. Am Ende der Simulation ergeben sich die folgenden Werte für $(t,r,\phi,\frac{d t(\tau)}{d\tau},\frac{d r(\tau)}{d\tau},\frac{d \phi(\tau)}{d\tau})$
solve_ivp-Methode:
print(Loes_ivp.y[0][-1],Loes_ivp.y[1][-1],Loes_ivp.y[2][-1],Loes_ivp.y[3][-1],Loes_ivp.y[4][-1],Loes_ivp.y[5][-1])
908.6596297139737 11.03794693307119 18.670947066104457 1.180602015359073 0.05196781395818743 0.03365171895832977
odeint-Methode:
print(Loes[-1, 0],Loes[-1, 1],Loes[-1, 2],Loes[-1, 3],Loes[-1, 4],Loes[-1, 5])
908.6598050881037 11.037902200863277 18.670989940158336 1.1806030235883063 0.05196848574304601 0.033651962768271884
Die Trajektorie des Körpers wird in dem unteren Bild dargestellt, wobei der rote Punkt den Anfangsort des Probekörpers zeigt.
ax = plt.gca()
ax.cla()
ax.set_ylabel(r"$\rm y \,[km]$")
ax.set_xlabel(r"$\rm x \,[km]$")
ax.plot(Loes[:, 1]*np.cos(Loes[:, 2]),Loes[:, 1]*np.sin(Loes[:, 2]),c="blue", linewidth=1.5, linestyle='-');
ax.scatter(Loes[0, 1]*np.cos(Loes[0, 2]),Loes[0, 1]*np.sin(Loes[0, 2]),c="red", s=60, marker='o');
BH = plt.Circle((0, 0), 2, color='black')
ax.add_patch(BH);
Wir werden in Kürze sehen, dass es Erhaltungsgrößen während der Bewegung des Probekörpers gibt. Während der gesamten Bewegung wird der Drehimpuls $l$ pro Masse $m$, die Teilchenenergie $e$ pro Masse und das Weglängenelement $ds^2$ erhalten. In der numerischen Simulation sind diese Größen ebenfalls in Näherung erhalten, wie man in den folgenden Abbildungen sieht.
import matplotlib.gridspec as gridspec
fig = plt.figure(figsize=(16,5))
gs = gridspec.GridSpec(1, 3, width_ratios=[1,1,1], wspace=0.25)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax3 = plt.subplot(gs[2])
ax1.set_ylabel(r"$\rm Drehimpuls \,\,l$")
ax1.set_xlabel(r"$\rm Eigenzeit \,\,\tau$")
ax2.set_ylabel(r"$\rm Energie \,\,e $")
ax2.set_xlabel(r"$\rm Eigenzeit \,\,\tau$")
ax3.set_ylabel(r"$\rm Weglaengenelement \,\,ds^2$")
ax3.set_xlabel(r"$\rm Eigenzeit \,\,\tau$")
data_l = Loes[:, 1]**2*Loes[:, 5]
data_E = (1-2/Loes[:, 1])*Loes[:, 3]
data_ds2 = (1-2/Loes[:, 1])*Loes[:, 3]**2-Loes[:, 4]**2/(1-2/Loes[:, 1])-Loes[:, 5]**2*Loes[:, 1]**2
ax1.plot(tauval,data_l,c="blue", linewidth=1.5, linestyle='-')
ax2.plot(tauval,data_E,c="blue", linewidth=1.5, linestyle='-')
ax3.plot(tauval,data_ds2,c="blue", linewidth=1.5, linestyle='-');
Klassifizierung unterschiedlicher Bahnbewegungen¶
Literatur:
- General relativity : An introduction for physicists von M. P. Hobson, G. P. Efstathiou und A. N. Lasenby
- Gravity : An introduction to Einstein's general relativity von James B. Hartle
- Allgemeine Relativitätstheorie von Torsten Fließbach
- Relativistic hydrodynamics von Luciano Rezzolla und Olindo Zanotti
Man kann zeigen (siehe z.B. Seite 206 in General relativity : An introduction for physicists by M. P. Hobson, G. P. Efstathiou and A. N. Lasenby), dass sich die erste und vierte Gleichung des Systems der Differenzialgleichungen der Geodätengleichung in die folgenden Gleichungen umschreiben lässt: $$ \begin{eqnarray} \hbox{1. Gleichung:}&& \frac{d}{d\tau} \left[ \left( 1 - \frac{2 M}{r} \right) \frac{dt}{d\tau} \right] = 0 \\ \rightarrow && \left( 1 - \frac{2 M}{r} \right) \frac{dt}{d\tau} = e = \hbox{const} \\ \hbox{4. Gleichung:}&& \frac{d}{d\tau} \left( r^2 \hbox{sin}^2(\theta) \frac{d\phi}{d\tau} \right) = 0 \\ \rightarrow && r^2 \hbox{sin}^2(\theta) \frac{d\phi}{d\tau} = l = \hbox{const} \quad , \end{eqnarray} $$ wobei die während der Bewegungen erhaltenen Größen e (Teilchenenergie pro Masse, $e=\frac{E}{m}$) und l (Drehimpuls pro Masse m, $l=\frac{L}{m}$) sich aus der Definition des Viererimpulses $p_\mu = m \, u_\mu = m \, g_{\mu\nu} \frac{dx^\nu}{d\tau} = m \, \left( g_{tt} \, \frac{dt}{d\tau}, g_{rr} \, \frac{dr}{d\tau}, g_{\theta\theta} \, \frac{d\theta}{d\tau}, g_{\phi\phi} \, \frac{d\phi}{d\tau} \right) $ ergeben (Voraussetzung: diagonale Form der Metrik). Für die Schwarzschildmetrik erhält man:
$$ \begin{eqnarray} && p_t = p_0 = m \frac{dx_0}{d\tau} = m g_{00} \frac{dx^0}{d\tau} = m g_{00} \frac{dt}{d\tau} = m \left( 1 - \frac{2 M}{r} \right) \frac{dt}{d\tau} = m \, e \\ && p_\phi = p_3 = m \frac{dx_3}{d\tau} = m g_{33} \frac{dx^3}{d\tau} = m g_{33} \frac{d\phi}{d\tau} = - m \left( r^2 \hbox{sin}^2(\theta) \right) \frac{d\phi}{d\tau} = - m \, l \end{eqnarray} $$
Die Klassifikation möglicher Bahnen von Probekörpern in vorgegebener Schwarzschild-Raumzeit kann mittels eines definierten effektiven Potenzials illustriert werden. Dieses Potenzial hängt von dem, bei der Bewegung erhaltenen Drehimpuls pro Masse m ab. Die im Zentralfeld möglichen Bewegungen werden mittels zweier erhaltener Größen (l: Drehimpuls pro Masse m und e: Energie pro Masse) charakterisiert. Die Definition des effektiven Potenzial erfolgt mittels der radialen, 2. Geodätengleichung. So definieren z.B. die Literaturangaben 1-3 (siehe Angaben oben) das effektive Potenzial wie folgt: $$ \begin{eqnarray} \hbox{2. Gleichung:}&\rightarrow&\frac{1}{2} \left( \frac{dr}{d\tau} \right)^2 + V(r,M,l) = \frac{1}{2} \left( e^2 -1 \right) \\ \hbox{wobei:} && V(r,M,l) = \frac{l^2}{2 r^2} \left( 1 - \frac{2 M}{r} \right) - \frac{M}{r} \quad . \end{eqnarray} $$ In der Literaturangabe 4 wird das Pozential V(r,M,l) hingegen wie folgt definiert:
$$ \begin{eqnarray} \hbox{2. Gleichung:}&\rightarrow& \left( \frac{dr}{d\tau} \right)^2 + \left( V(r,M,l) \right)^2 = e^2 \\ \hbox{wobei:} && V(r,M,l) = \sqrt{ \left( 1 - \frac{2 M}{r} \right) \left( 1 + \frac{l^2}{r^2} \right) } \quad . \end{eqnarray} $$
l = symbols('l')
VeffFb = -M/r+l**2/(2*r**2)-M*l**2/r**3
VeffRez = sqrt(1-2*M/r)*sqrt(1+l**2/r**2)
setl = r0**2*dphi0
setE = dt0*((1-2*set_M/r0))
print("Drehimpuls l")
print(setl)
print("Energie e")
print(setE)
VeffFbl = lambdify((r),VeffFb.subs({(l,setl),(M,set_M)}))
VeffRezl = lambdify((r),VeffRez.subs({(l,setl),(M,set_M)}))
rval = np.linspace(2.85,30, 1001)
Drehimpuls l 4.1000000000000005 Energie e 0.9666850572963256
fig = plt.figure(figsize=(16,8))
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.21)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_ylabel(r"$\rm V_1(r)$")
ax1.set_xlabel(r"$\rm Radius \, r \,[km]$")
ax2.set_ylabel(r"$\rm V_2(r)$")
ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
ax1.plot(rval,VeffRezl(rval),c="blue", linewidth=1.5, linestyle='-');
ax2.plot(rval,VeffFbl(rval),c="red", linewidth=1.5, linestyle='-');
ax1.scatter(r0,VeffRezl(r0),c="black", s=80, marker='o');
ax2.scatter(r0,VeffFbl(r0),c="black", s=80, marker='o');
Die obere Abbildung zeigt (in der Nomenklatur des 4.Buches (Diagramm links) bzw. des 1. bis 3. Buches (Diagramm rechts) ) das effektive Pozential als Funktion des Radius bei den zuvor gewählten Anfangswerten (schwarze Punkte auf den Kurven).
Wir wollen uns nun die Bewegung des Probekörpers in einer Animation darstellen:
import matplotlib.animation as animation
from IPython.display import HTML
anzframes = 100
step = int(N/anzframes-1)
xylim = 18
fig = plt.figure(figsize=(11, 5))
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.3)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
def init():
ax1.scatter(Loes[0, 1]*np.cos(Loes[0, 2]),Loes[0, 1]*np.sin(Loes[0, 2]),c="red", s=60, marker='o')
BH = plt.Circle((0, 0), 2, color='black')
ax1.add_patch(BH)
ax2.plot(rval,VeffRezl(rval),c="blue", linewidth=1.5, linestyle='-')
ax2.scatter(r0,VeffRezl(r0),c="red", s=80, marker='o')
ax1.set_ylabel(r"$\rm y \,[km]$")
ax1.set_xlabel(r"$\rm x \,[km]$")
ax2.set_ylabel(r"$\rm V(r)$")
ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
ax1.set_xlim(-xylim, xylim)
ax1.set_ylim(-xylim, xylim)
return fig,
def animate(i):
ax1.cla()
ax2.cla()
ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),c="blue", linewidth=0.5, linestyle='-')
ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),c="red", s=60, marker='o')
ax2.plot(rval,VeffRezl(rval),c="blue", linewidth=1.5, linestyle='-')
ax2.scatter(Loes[step*i, 1],VeffRezl(r0),c="red", s=80, marker='o')
BH = plt.Circle((0, 0), 2, color='black')
ax1.add_patch(BH);
ax1.set_ylabel(r"$\rm y \,[km]$")
ax1.set_xlabel(r"$\rm x \,[km]$")
ax2.set_ylabel(r"$\rm V(r)$")
ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
ax1.set_xlim(-xylim, xylim)
ax1.set_ylim(-xylim, xylim)
return fig,
ani = animation.FuncAnimation(fig,animate,init_func=init,frames=anzframes,interval=100)
plt.close(fig)
HTML(ani.to_html5_video())
bzw. ohne Verwendung von ax1.cla() und ax2.cla() in jedem Animationsschritt und mit einer Trennung von Animationsobjekten und statischen Elementen.
# Figure-Setup
fig = plt.figure(figsize=(11, 5))
gs = gridspec.GridSpec(1, 2, width_ratios=[1, 1], wspace=0.3)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
# Statische Elemente der Animation
ax1.set_ylabel(r"$\rm y \,[km]$")
ax1.set_xlabel(r"$\rm x \,[km]$")
ax1.set_xlim(-xylim, xylim)
ax1.set_ylim(-xylim, xylim)
ax2.set_ylabel(r"$\rm V(r)$")
ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
# Statische Potentialkurve
ax2.plot(rval, VeffRezl(rval), c="blue", linewidth=1.5, linestyle='-')
# Statisches Schwarzes Loch
BH = plt.Circle((0, 0), 2, color='black')
ax1.add_patch(BH)
# Animationsobjekte initialisieren
line_orbit, = ax1.plot([], [], c="blue", linewidth=0.5, linestyle='-')
scatter_body = ax1.scatter([], [], c="red", s=60, marker='o')
scatter_pot = ax2.scatter([], [], c="red", s=80, marker='o')
# Simulationswerte der Trajektorie
x_all = Loes[:, 1] * np.cos(Loes[:, 2])
y_all = Loes[:, 1] * np.sin(Loes[:, 2])
def init():
line_orbit.set_data([], [])
scatter_body.set_offsets(np.c_[x_all[0], y_all[0]])
scatter_pot.set_offsets(np.c_[Loes[0, 1], VeffRezl(Loes[0, 1])])
return line_orbit, scatter_body, scatter_pot
def animate(i):
idx = step * i
line_orbit.set_data(x_all[:idx], y_all[:idx])
scatter_body.set_offsets(np.c_[x_all[idx], y_all[idx]])
scatter_pot.set_offsets(np.c_[Loes[idx, 1], VeffRezl(Loes[0, 1])])
return line_orbit, scatter_body, scatter_pot
ani = animation.FuncAnimation(
fig, animate, init_func=init, frames=anzframes,
interval=100, blit=True
)
plt.close(fig)
HTML(ani.to_html5_video())
Wir betrachten uns die gleiche Animation in einem eingebetteten Diagramm der räumlichen Hypersphäre der Mannigfaltigkeit:
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
x,y = symbols('x,y')
zBett = sqrt(8*M*(r-2*M))
R = np.linspace(2.00001, 20, 300)
PHI = np.linspace(0, 2*np.pi, 300)
R, PHI = np.meshgrid(R, PHI)
X = R*np.cos(PHI)
Y = R*np.sin(PHI)
zlam = lambdify((x,y), zBett.subs({(M,set_M),(r,sqrt(x**2+y**2))}))
ZF = np.array(list(np.real(zlam(X,Y))))
anzframes = 50
step = int(N/anzframes-1)
xylim = 18
zlim = 15
fig = plt.figure(figsize=(11,5))
gs = gridspec.GridSpec(1, 2, width_ratios=[2,1], wspace=0.21)
ax1 = plt.subplot(gs[0],projection='3d')
ax2 = plt.subplot(gs[1])
ax1.xaxis.pane.fill = False
ax1.yaxis.pane.fill = False
ax1.zaxis.pane.fill = False
ax1.xaxis.pane.set_edgecolor('white')
ax1.yaxis.pane.set_edgecolor('white')
ax1.zaxis.pane.set_edgecolor('white')
def init():
ax1.plot_surface(X, Y, ZF, cmap=cm.nipy_spectral, linewidth=0, alpha=0.5)
ax1.scatter(Loes[0, 1]*np.cos(Loes[0, 2]),Loes[0, 1]*np.sin(Loes[0, 2]),np.real(zlam(Loes[0, 1]*np.cos(Loes[0, 2]),Loes[0, 1]*np.sin(Loes[0, 2]))),c="black", s=60, marker='o');
ax1.contourf(X, Y, ZF,100, zdir='z', offset=0, cmap=cm.nipy_spectral,alpha=0.5)
ax1.scatter(Loes[0, 1]*np.cos(Loes[0, 2]),Loes[0, 1]*np.sin(Loes[0, 2]),0.1,c="black", s=40, marker='o');
ax2.plot(rval,VeffRezl(rval),c="black", linewidth=1.5, linestyle='-');
ax2.scatter(r0,VeffRezl(r0),c="black", s=80, marker='o');
ax1.set_ylabel(r"$\rm y \,[km]$")
ax1.set_xlabel(r"$\rm x \,[km]$")
ax1.set_zlabel("Visualisierung z");
ax2.set_ylabel(r"$\rm V(r)$")
ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
ax1.view_init(azim=35, elev=20)
ax1.set_xlim(-xylim, xylim)
ax1.set_ylim(-xylim, xylim)
ax1.set_zlim(0, 15)
ax1.grid(False)
return fig,
def animate(i):
ax1.cla()
ax2.cla()
ax1.plot_surface(X, Y, ZF, cmap=cm.nipy_spectral, linewidth=0, alpha=0.5)
ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),np.real(zlam(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]))),c="black", linewidth=1, linestyle='-');
ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),np.real(zlam(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]))),c="black", s=60, marker='o');
ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),0.1,c="black", s=40, marker='o');
ax1.contourf(X, Y, ZF,100, zdir='z', offset=0, cmap=cm.nipy_spectral,alpha=0.5)
ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),0.1,c="black", linewidth=2.5, linestyle='-');
ax2.plot(rval,VeffRezl(rval),c="black", linewidth=1.5, linestyle='-');
ax2.scatter(Loes[step*i, 1],VeffRezl(r0),c="black", s=80, marker='o');
ax1.set_ylabel(r"$\rm y \,[km]$")
ax1.set_xlabel(r"$\rm x \,[km]$")
ax1.set_zlabel("Einbettung z");
ax2.set_ylabel(r"$\rm V(r)$")
ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
ax1.view_init(azim=35, elev=20)
ax1.set_xlim(-xylim, xylim)
ax1.set_ylim(-xylim, xylim)
ax1.set_zlim(0, 15)
ax1.grid(False)
return fig,
ani = animation.FuncAnimation(fig,animate,init_func=init,frames=anzframes,interval=100)
plt.close(fig)
HTML(ani.to_html5_video())
bzw. ohne Verwendung von ax1.cla() und ax2.cla() in jedem Animationsschritt und mit einer Trennung von Animationsobjekten und statischen Elementen.
anzframes = 50
step = int(N/anzframes-1)
xylim = 18
zlim = 15
# Figure-Setup
fig = plt.figure(figsize=(11, 5))
gs = gridspec.GridSpec(1, 2, width_ratios=[2, 1], wspace=0.21)
ax1 = plt.subplot(gs[0], projection='3d')
ax2 = plt.subplot(gs[1])
for axis in [ax1.xaxis, ax1.yaxis, ax1.zaxis]:
axis.pane.fill = False
axis.pane.set_edgecolor('white')
# Achsen-Konfiguration
ax1.set_ylabel(r"$\rm y \,[km]$")
ax1.set_xlabel(r"$\rm x \,[km]$")
ax1.set_zlabel("Einbettung z")
ax1.set_xlim(-xylim, xylim)
ax1.set_ylim(-xylim, xylim)
ax1.set_zlim(0, zlim)
ax1.view_init(azim=35, elev=20)
ax1.grid(False)
ax2.set_ylabel(r"$\rm V(r)$")
ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
ax1.plot_surface(X, Y, ZF, cmap=cm.nipy_spectral, linewidth=0, alpha=0.5)
ax1.contourf(X, Y, ZF, 100, zdir='z', offset=0, cmap=cm.nipy_spectral, alpha=0.5)
ax2.plot(rval, VeffRezl(rval), c="black", linewidth=1.5, linestyle='-')
line_3d, = ax1.plot([], [], [], c="black", linewidth=1, linestyle='-') # Trajektorie auf Fläche
line_shadow, = ax1.plot([], [], [], c="black", linewidth=2.5, linestyle='-') # Schatten auf dem Boden
scatter_3d = ax1.scatter([], [], [], c="black", s=60, marker='o')
scatter_shadow = ax1.scatter([], [], [], c="black", s=40, marker='o')
scatter_pot = ax2.scatter([], [], c="black", s=80, marker='o')
x_all = Loes[:, 1] * np.cos(Loes[:, 2])
y_all = Loes[:, 1] * np.sin(Loes[:, 2])
z_all = np.real(zlam(x_all, y_all))
z_shadow = np.full_like(z_all, 0.1)
def init():
line_3d.set_data([], [])
line_3d.set_3d_properties([])
line_shadow.set_data([], [])
line_shadow.set_3d_properties([])
scatter_3d._offsets3d = ([], [], [])
scatter_shadow._offsets3d = ([], [], [])
scatter_pot.set_offsets(np.c_[[], []])
return fig,
def animate(i):
idx = step * i
line_3d.set_data(x_all[:idx], y_all[:idx])
line_3d.set_3d_properties(z_all[:idx])
line_shadow.set_data(x_all[:idx], y_all[:idx])
line_shadow.set_3d_properties(z_shadow[:idx])
scatter_3d._offsets3d = ([x_all[idx]], [y_all[idx]], [z_all[idx]])
scatter_shadow._offsets3d = ([x_all[idx]], [y_all[idx]], [0.1])
scatter_pot.set_offsets(np.c_[Loes[idx, 1], VeffRezl(Loes[0, 1])])
return fig,
ani = animation.FuncAnimation(
fig, animate, init_func=init, frames=anzframes,
interval=100, blit=True
)
plt.close(fig)
HTML(ani.to_html5_video())
Im Folgenden werden die unterschiedlichen Bahnbewegungen um ein schwarzes Loch klassifiziert. Die einzelnen Bewegungstypen lassen sich physikalisch über die Energie $e$ und den Drehimpuls $l$ des Testteilchens relativ zum effektiven Potenzial $V(r,M,l)$ wie folgt charakterisieren:
A) Kreisförmige Bahnen: Diese Bahnen existieren an den Extrempunkten des effektiven Potenzials. Man unterscheidet hierbei zwischen stabilen kreisförmigen Bahnen (im Minimum des Potenzials, oberhalb des sogenannten ISCO – Innermost Stable Circular Orbit) und instabilen kreisförmigen Bahnen (auf dem Maximum des Potentials).
t0 = 0
phi0 = 0
dr0 = 0
setl = 4.1
setVeffRez = VeffRez.subs({(l,setl),(M,set_M)})
EqCircOrb = Eq(diff(setVeffRez,r),0)
r0A = float(solve(EqCircOrb)[1])
dphi0A = setl/r0A**2
dt0A = float(LoesEq_ds[0].subs({(r,r0A),(dr,dr0),(dphi,dphi0A),(M,set_M)}))
setEA = dt0A*((1-2*set_M/r0A))
initialvalA = np.array([t0,r0A,phi0,dt0A,dr0,dphi0A])
B) Elliptische Bahnen / Rosettenbahnen: Wenn die Energie des Probekörpers größer als das Minimum, aber kleiner als das Maximum des Potenzials ist, resultiert eine gebundene Bewegung. Aufgrund der relativistischen Effekte schließt sich diese Bahn jedoch nicht wie bei Newton (Kepler-Ellipse), sondern führt zu einer markanten Periheldrehung (Einstein-Präzession), wodurch die typische Rosettenform entsteht.
r0B = 8
dphi0B = setl/r0B**2
dt0B = float(LoesEq_ds[0].subs({(r,r0B),(dr,dr0),(dphi,dphi0B),(M,set_M)}))
setEB = dt0B*(1-2*set_M/r0B)
initialvalB = np.array([t0,r0B,phi0,dt0B,dr0,dphi0B])
C) Ungebunden parabolische Bahnen: Diese stellen den Grenzfall dar, bei dem der Probekörper exakt die Fluchtenergie besitzt: $e = \lim\limits_{r \to \infty}V(r,M,l) = 1$ . Das Teilchen kommt theoretisch aus dem Unendlichen mit der Geschwindigkeit Null an, wird vom Schwarzen Loch abgelenkt und entweicht wieder ins Unendliche.
setEC = float(limit_seq(setVeffRez, r))
EqE1 = Eq(setVeffRez,setEC)
r0C = float(solve(EqE1,r)[1])
dphi0C = setl/r0C**2
dt0C = float(LoesEq_ds[0].subs({(r,r0C),(dr,dr0),(dphi,dphi0C),(M,set_M)}))
initialvalC = np.array([t0,r0C,phi0,dt0C,dr0,dphi0C])
D) Ungebunden hyperbolische Bahnen: Hier besitzt das Teilchen eine Überschussenergie: $e > \lim\limits_{r \to \infty}V(r,M,l) = 1$. Es nähert sich mit einer hohen Anfangsgeschwindigkeit aus dem Unendlichen, erfährt in der Nähe des Ereignishorizonts eine starke gravitative Streuung (oft mit extremen relativistischen Ablenkungswinkeln) und verlässt das System wieder.
setED = 1.005
EqE1005 = Eq(setVeffRez,setED)
r0D = solve(EqE1005,r)[2]
r0D = float(re(r0D))
dphi0D = setl/r0D**2
dt0D = float(LoesEq_ds[0].subs({(r,r0D),(dr,dr0),(dphi,dphi0D),(M,set_M)}))
initialvalD = np.array([t0,r0D,phi0,dt0D,dr0,dphi0D])
r0D
E) Einfangbahnen / Capture Orbits: Wenn die Energie oder der Drehimpuls des Teilchens so beschaffen sind, dass es die Potenzialbarriere des Schwarzen Lochs überwindet (Energie $e$ ist größer als das Maximum des Potentials), gibt es keinen klassischen Umkehrpunkt mehr. Das Teilchen spiralisiert unaufhaltsam über den Ereignishorizont hinaus und wird vom Schwarzen Loch dauerhaft verschluckt.
r0E = 225
dphi0E = setl/r0E**2
dr0E=-0.185388257
dt0E = float(LoesEq_ds[0].subs({(r,r0E),(dr,dr0E),(dphi,dphi0E),(M,set_M)}))
setEE = dt0E*((1-2*set_M/r0E))
initialvalE = np.array([t0,r0E,phi0,dt0E,dr0E,dphi0E])
Wir simulieren alle fünf Bahnbewegungen.
N = 50000
tauval = np.linspace(0, 1068, N+1)
LoesA = integrate.odeint(DGLsys, initialvalA, tauval)
LoesB = integrate.odeint(DGLsys, initialvalB, tauval)
LoesC = integrate.odeint(DGLsys, initialvalC, tauval)
LoesD = integrate.odeint(DGLsys, initialvalD, tauval)
LoesE = integrate.odeint(DGLsys, initialvalE, tauval)
Loes_vec = [LoesA,LoesB,LoesC,LoesD,LoesE]
zBett = sqrt(8*M*(r-2*M))
R = np.linspace(2.00001, 40, 1000)
PHI = np.linspace(0, 2*np.pi, 300)
R, PHI = np.meshgrid(R, PHI)
X = R*np.cos(PHI)
Y = R*np.sin(PHI)
zlam = lambdify((x,y), zBett.subs({(M,1),(r,sqrt(x**2+y**2))}))
ZF = np.array(list(np.real(zlam(X,Y))))
rval = np.linspace(2.85,40, 1001)
step = 500
anzframes = 100
xylim = 35
zlim = 20
fig = plt.figure(figsize=(11,5))
gs = gridspec.GridSpec(1, 2, width_ratios=[2,1], wspace=0.21)
ax1 = plt.subplot(gs[0],projection='3d')
ax2 = plt.subplot(gs[1])
ax1.xaxis.pane.fill = False
ax1.yaxis.pane.fill = False
ax1.zaxis.pane.fill = False
ax1.xaxis.pane.set_edgecolor('white')
ax1.yaxis.pane.set_edgecolor('white')
ax1.zaxis.pane.set_edgecolor('white')
def init():
ax1.plot_surface(X, Y, ZF, cmap=cm.nipy_spectral, linewidth=0, alpha=0.5)
ax1.contourf(X, Y, ZF,100, zdir='z', offset=0, cmap=cm.nipy_spectral,alpha=0.5)
ax1.scatter(LoesA[0, 1]*np.cos(LoesA[0, 2]),LoesA[0, 1]*np.sin(LoesA[0, 2]),np.real(zlam(LoesA[0, 1]*np.cos(LoesA[0, 2]),LoesA[0, 1]*np.sin(LoesA[0, 2]))),c="blue", s=60, marker='o');
ax1.scatter(LoesA[0, 1]*np.cos(LoesA[0, 2]),LoesA[0, 1]*np.sin(LoesA[0, 2]),0.1,c="blue", s=40, marker='o');
ax1.scatter(LoesB[0, 1]*np.cos(LoesB[0, 2]),LoesB[0, 1]*np.sin(LoesB[0, 2]),np.real(zlam(LoesB[0, 1]*np.cos(LoesB[0, 2]),LoesB[0, 1]*np.sin(LoesB[0, 2]))),c="red", s=60, marker='o');
ax1.scatter(LoesB[0, 1]*np.cos(LoesB[0, 2]),LoesB[0, 1]*np.sin(LoesB[0, 2]),0.1,c="red", s=40, marker='o');
ax1.scatter(LoesC[0, 1]*np.cos(LoesC[0, 2]),LoesC[0, 1]*np.sin(LoesC[0, 2]),np.real(zlam(LoesC[0, 1]*np.cos(LoesC[0, 2]),LoesC[0, 1]*np.sin(LoesC[0, 2]))),c="green", s=60, marker='o');
ax1.scatter(LoesC[0, 1]*np.cos(LoesC[0, 2]),LoesC[0, 1]*np.sin(LoesC[0, 2]),0.1,c="green", s=40, marker='o');
ax1.scatter(LoesD[0, 1]*np.cos(LoesD[0, 2]),LoesD[0, 1]*np.sin(LoesD[0, 2]),np.real(zlam(LoesD[0, 1]*np.cos(LoesD[0, 2]),LoesD[0, 1]*np.sin(LoesD[0, 2]))),c="grey", s=60, marker='o');
ax1.scatter(LoesD[0, 1]*np.cos(LoesD[0, 2]),LoesD[0, 1]*np.sin(LoesD[0, 2]),0.1,c="grey", s=40, marker='o');
ax1.scatter(LoesE[0, 1]*np.cos(LoesE[0, 2]),LoesE[0, 1]*np.sin(LoesE[0, 2]),np.real(zlam(LoesE[0, 1]*np.cos(LoesE[0, 2]),LoesE[0, 1]*np.sin(LoesE[0, 2]))),c="black", s=60, marker='o');
ax1.scatter(LoesE[0, 1]*np.cos(LoesE[0, 2]),LoesE[0, 1]*np.sin(LoesE[0, 2]),0.1,c="black", s=40, marker='o');
ax2.plot(rval,VeffRezl(rval),c="black", linewidth=1.5, linestyle='-');
ax2.scatter(r0,VeffRezl(r0),c="black", s=80, marker='o');
ax1.set_ylabel(r"$\rm y \,[km]$")
ax1.set_xlabel(r"$\rm x \,[km]$")
ax1.set_zlabel("Visualisierung z");
ax2.set_ylabel(r"$\rm V(r)$")
ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
ax1.view_init(azim=-130, elev=35)
ax1.set_xlim(-xylim, xylim)
ax1.set_ylim(-xylim, xylim)
ax1.set_zlim(0, zlim)
ax2.set_xlim(2, xylim)
ax1.grid(False)
return fig,
def animate(i):
ax1.cla()
ax2.cla()
ax1.plot_surface(X, Y, ZF, cmap=cm.nipy_spectral, linewidth=0, alpha=0.5)
ax1.contourf(X, Y, ZF,100, zdir='z', offset=0, cmap=cm.nipy_spectral,alpha=0.5)
Loes = LoesA
ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),np.real(zlam(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]))),c="blue", linewidth=1, linestyle='-');
ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),0.1,c="blue", linewidth=2.5, linestyle='-');
ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),np.real(zlam(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]))),c="blue", s=60, marker='o');
ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),0.1,c="blue", s=40, marker='o');
Loes = LoesB
ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),np.real(zlam(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]))),c="red", linewidth=1, linestyle='-');
ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),0.1,c="red", linewidth=2.5, linestyle='-');
ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),np.real(zlam(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]))),c="red", s=60, marker='o');
ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),0.1,c="red", s=40, marker='o');
Loes = LoesC
ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),np.real(zlam(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]))),c="green", linewidth=1, linestyle='-');
ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),0.1,c="green", linewidth=2.5, linestyle='-');
ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),np.real(zlam(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]))),c="green", s=60, marker='o');
ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),0.1,c="green", s=40, marker='o');
Loes = LoesD
ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),np.real(zlam(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]))),c="grey", linewidth=1, linestyle='-');
ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),0.1,c="grey", linewidth=2.5, linestyle='-');
ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),np.real(zlam(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]))),c="grey", s=60, marker='o');
ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),0.1,c="grey", s=40, marker='o');
Loes = LoesE
ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),np.real(zlam(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]))),c="black", linewidth=1, linestyle='-');
ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),0.1,c="black", linewidth=2.5, linestyle='-');
ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),np.real(zlam(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]))),c="black", s=60, marker='o');
ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),0.1,c="black", s=40, marker='o');
ax2.plot(rval,VeffRezl(rval),c="black", linewidth=1.5, linestyle='-');
ax2.scatter(LoesA[step*i, 1],VeffRezl(LoesA[0, 1]),c="blue", s=80, marker='o');
ax2.scatter(LoesB[step*i, 1],VeffRezl(LoesB[0, 1]),c="red", s=80, marker='o');
ax2.plot(LoesB[:step*i, 1],VeffRezl(LoesB[0, 1]*LoesB[:step*i, 1]/LoesB[:step*i, 1]),c="red", linewidth=1, linestyle=':');
ax2.scatter(LoesC[step*i, 1],VeffRezl(LoesC[0, 1]),c="green", s=80, marker='o');
ax2.plot(LoesC[:step*i, 1],VeffRezl(LoesC[0, 1]*LoesC[:step*i, 1]/LoesC[:step*i, 1]),c="green", linewidth=1, linestyle=':');
ax2.scatter(LoesD[step*i, 1],VeffRezl(LoesD[0, 1]),c="grey", s=80, marker='o');
ax2.plot(LoesD[:step*i, 1],VeffRezl(LoesD[0, 1]*LoesD[:step*i, 1]/LoesD[:step*i, 1]),c="grey", linewidth=1, linestyle=':');
ax2.scatter(LoesE[step*i, 1],setEE,c="black", s=80, marker='o');
ax2.plot(LoesE[:step*i, 1],setEE*LoesE[:step*i, 1]/LoesE[:step*i, 1],c="black", linewidth=1, linestyle=':');
ax1.set_ylabel(r"$\rm y \,[km]$")
ax1.set_xlabel(r"$\rm x \,[km]$")
ax1.set_zlabel("Einbettung z");
ax2.set_ylabel(r"$\rm V(r)$")
ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
ax1.view_init(azim=-130, elev=35)
ax1.set_xlim(-xylim, xylim)
ax1.set_ylim(-xylim, xylim)
ax1.set_zlim(0, zlim)
ax2.set_xlim(2, xylim)
ax1.grid(False)
return fig,
ani = animation.FuncAnimation(fig,animate,init_func=init,frames=anzframes,interval=100)
plt.close(ani._fig)
HTML(ani.to_html5_video())
bzw. ohne Verwendung von ax1.cla() und ax2.cla() in jedem Animationsschritt und mit einer Trennung von Animationsobjekten und statischen Elementen.
step = 500
anzframes = 100
xylim = 35
zlim = 20
# Figure-Setup
fig = plt.figure(figsize=(11, 5))
gs = gridspec.GridSpec(1, 2, width_ratios=[2, 1], wspace=0.21)
ax1 = plt.subplot(gs[0], projection='3d')
ax2 = plt.subplot(gs[1])
for axis in [ax1.xaxis, ax1.yaxis, ax1.zaxis]:
axis.pane.fill = False
axis.pane.set_edgecolor('white')
ax1.set_ylabel(r"$\rm y \,[km]$")
ax1.set_xlabel(r"$\rm x \,[km]$")
ax1.set_zlabel("Einbettung z")
ax1.set_xlim(-xylim, xylim)
ax1.set_ylim(-xylim, xylim)
ax1.set_zlim(0, zlim)
ax1.view_init(azim=-130, elev=35)
ax1.grid(False)
v_values = VeffRezl(rval)
ax2.set_ylabel(r"$\rm V(r)$")
ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
ax2.set_xlim(np.min(rval), np.max(rval))
# Statische Hintergründe zeichnen
ax1.plot_surface(X, Y, ZF, cmap=cm.nipy_spectral, linewidth=0, alpha=0.5)
ax1.contourf(X, Y, ZF, 100, zdir='z', offset=0, cmap=cm.nipy_spectral, alpha=0.5)
ax2.plot(rval, v_values, c="black", linewidth=1.5, linestyle='-')
# Animationsobjekte
lines_3d = []
lines_shadow = []
scatters_3d = []
scatters_shadow = []
scatters_pot = []
lines_pot = []
colors = ["blue", "red", "green", "gray", "black"]
for c in colors:
l3d, = ax1.plot([], [], [], c=c, linewidth=1.2, linestyle='-')
lsh, = ax1.plot([], [], [], c=c, linewidth=2.0, linestyle='-', alpha=0.6)
lines_3d.append(l3d)
lines_shadow.append(lsh)
scatters_3d.append(ax1.scatter([], [], [], c=c, s=50, marker='o'))
scatters_shadow.append(ax1.scatter([], [], [], c=c, s=30, marker='o'))
scatters_pot.append(ax2.scatter([], [], c=c, s=60, marker='o', zorder=5))
lpot, = ax2.plot([], [], c=c, linewidth=1, linestyle=':', alpha=0.7)
lines_pot.append(lpot)
x_all_vec = []
y_all_vec = []
z_all_vec = []
for Loes in Loes_vec:
x_all = Loes[:, 1] * np.cos(Loes[:, 2])
y_all = Loes[:, 1] * np.sin(Loes[:, 2])
z_all = np.real(zlam(x_all, y_all))
x_all_vec.append(x_all)
y_all_vec.append(y_all)
z_all_vec.append(z_all)
def init():
for j in range(5):
lines_3d[j].set_data([], [])
lines_3d[j].set_3d_properties([])
lines_shadow[j].set_data([], [])
lines_shadow[j].set_3d_properties([])
scatters_3d[j]._offsets3d = ([], [], [])
scatters_shadow[j]._offsets3d = ([], [], [])
scatters_pot[j].set_offsets(np.empty((0, 2)))
lines_pot[j].set_data([], [])
return fig,
def animate(i):
idx = step * i
for j in range(5):
Loes = Loes_vec[j]
x_all = x_all_vec[j]
y_all = y_all_vec[j]
z_all = z_all_vec[j]
lines_3d[j].set_data(x_all[:idx], y_all[:idx])
lines_3d[j].set_3d_properties(z_all[:idx])
lines_shadow[j].set_data(x_all[:idx], y_all[:idx])
lines_shadow[j].set_3d_properties(np.full_like(x_all[:idx], 0.1))
scatters_3d[j]._offsets3d = ([x_all[idx]], [y_all[idx]], [z_all[idx]])
scatters_shadow[j]._offsets3d = ([x_all[idx]], [y_all[idx]], [0.1])
r_hist = Loes[:idx, 1]
if(j < 4):
v_hist = VeffRezl(Loes[0, 1])
else:
v_hist = setEE
if(len(r_hist) > 0):
scatters_pot[j].set_offsets(np.c_[r_hist[-1], v_hist])
lines_pot[j].set_data(r_hist, [v_hist]*len(r_hist))
return fig,
ani = animation.FuncAnimation(
fig, animate, init_func=init, frames=anzframes,
interval=100, blit=True
)
plt.close(fig)
HTML(ani.to_html5_video())
In der oberen Animation erkennt man neben den gebundenen kreisförmigen (A: blau) und elliptischen (B: rot) Bahnen, den parabolischen (C: grün) und hyperbolischen (D: grau) Bahnverläufen auch eine durch das schwarze Loch eingefangene Bahn (E: capture orbit, schwarz). Diese ist in der unteren Abbildung nochmals separat veranschaulicht.
Loes = LoesE
ax = plt.gca()
ax.cla()
ax.set_ylabel(r"$\rm y \,[km]$")
ax.set_xlabel(r"$\rm x \,[km]$")
ax.plot(Loes[:, 1]*np.cos(Loes[:, 2]),Loes[:, 1]*np.sin(Loes[:, 2]),c="black", linewidth=1.5, linestyle='-');
ax.scatter(Loes[0, 1]*np.cos(Loes[0, 2]),Loes[0, 1]*np.sin(Loes[0, 2]),c="black", s=60, marker='o');
ax.axis('equal');
Die Visualisierung kann man auch mittels der "Plotly Python Open Source Graphing Library" erstellen ( siehe https://plotly.com/python/ ). In der folgenden interaktiven Abbildung wird der "Capture orbit" und der "Elliptische Orbit" visualisiert. Zusätzlich wird die Eigen- und Koordinatenzeit ausgegeben, wenn man mit der Maus genau auf eine der beiden Probekörper-Trajektorien fokussiert.
import plotly.graph_objects as go
trichter = go.Surface(x=X, y=Y, z=ZF, colorscale="spectral_r", showscale=True, name="Embedded Surface", hoverinfo='none')
Loes = LoesE
labels = []
for i in range(len(Loes[:, 0])):
labels.append('Eigenzeit '+str(round(tauval[i], 2))+', Koordinatenzeit = '+str(round(Loes[i, 0], 2)))
capture_orbit =go.Scatter3d(
x = Loes[:, 1]*np.cos(Loes[:, 2]),
y = Loes[:, 1]*np.sin(Loes[:, 2]),
z = np.real(zlam(Loes[:, 1]*np.cos(Loes[:, 2]),Loes[:, 1]*np.sin(Loes[:, 2]))),
name = "Capture Orbit",
mode = 'lines',
line = dict(color='black', width=8),
text = labels,
hoverinfo = 'text'
)
Loes = LoesB
labels = []
for i in range(len(Loes[:, 0])):
labels.append('Eigenzeit '+str(round(tauval[i], 2))+', Koordinatenzeit = '+str(round(Loes[i, 0], 2)))
elliptic_orbit =go.Scatter3d(
x = Loes[:, 1]*np.cos(Loes[:, 2]),
y = Loes[:, 1]*np.sin(Loes[:, 2]),
z = np.real(zlam(Loes[:, 1]*np.cos(Loes[:, 2]),Loes[:, 1]*np.sin(Loes[:, 2]))),
name = "Elliptische Bahn",
mode = 'lines',
line = dict(color='red', width=8),
text = labels,
hoverinfo = 'text'
)
layout = go.Layout(autosize=True,
width=700, height=700,
margin=dict(l=15, r=50, b=65, t=90),
scene_aspectmode='cube',
scene = dict(
xaxis_showbackground = False,
yaxis_showbackground = False,
zaxis_showbackground = False,
xaxis_gridcolor = "lightgrey",
yaxis_gridcolor = "lightgrey",
zaxis_gridcolor = "lightgrey",
xaxis_range = [-40,40],
yaxis_range = [-40,40],
zaxis_range = [0,20],
xaxis_title = 'x',
yaxis_title = 'y',
zaxis_title = 'f'),
legend=dict(
yanchor = "top",
y=0.99,
xanchor = "left",
x=0.01)
)
data = [trichter,elliptic_orbit,capture_orbit]
fig = go.Figure(data = data, layout = layout)
fig.show(renderer='iframe')