Meervoudige massa’s in de ORT#
De voorgaande hoofdstukken beschrijven de zwaartekracht van één bolvormige massa met de formule:
In dit hoofdstuk breiden we het model uit naar meerdere massa’s. De sleutel is het onderscheid tussen twee grootheden:
\(v_{grav}\) — een scalair: de totale gravitatie-energieparameter
\(\vec{v}_{netto}\) — een vector: de netto gerichte stroming
De algemene formule:
import sys, pathlib
sys.path.insert(0, str(pathlib.Path().resolve().parent / 'shared'))
from ort_core import C, G, MultiGravity, SphericalMassDistribution, M_SUN
import math
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
%matplotlib inline
1 Wat is \(v_{grav}\)?#
1.1 Energieparameter, geen snelheid#
De grootheid \(v_{grav} = \sqrt{2GM/r}\) wordt vaak “ontsnappingssnelheid” genoemd, maar dat is het niet. Het is:
De snelheid die Newtons kinetische energie (\(\frac{1}{2}m_0 v^2\)) toekent aan de potentiële energie. In de ORT hoort bij dezelfde energie een lagere snelheid (\(v_{esc} = c\sin\theta\) met \(E_k = m_0c^2(1/\cos\theta - 1)\)).
Toch werkt \(v_{grav} = \sqrt{2GM/r}\) exact:
Dit matcht de Schwarzschild-metriek. De reden: de Gullstrand-Painlevé coördinaten gebruiken precies dezelfde grootheid als “riviersnelheid” en reproduceren de Schwarzschild-geometrie exact.
De Schwarzschild-straal zelf is gedefinieerd door dezelfde relatie: \(\frac{1}{2}m_0 c^2 = GM m_0/r_s\), dus \(r_s = 2GM/c^2\).
1.2 Energieën tellen op#
Omdat \(v_{grav}^2 \propto E_{pot}\) en potentiële energieën optellen:
Dit is een scalaire som. De \(v_{grav}^2\)-waarden van individuele massa’s tellen op als energieën, niet als vectoren.
De gravitatievelden tellen op als vectoren:
Dit onderscheid — scalaire potentiaal vs. vectorieel veld — is de sleutel.
2 De formule voor meerdere massa’s#
2.1 \(v_{grav}\) en \(\vec{v}_{netto}\)#
Voor meerdere massa’s zijn er twee grootheden:
\(v_{grav}\) |
\(\vec{v}_{netto}\) |
|
|---|---|---|
Type |
scalair |
vector |
Definitie |
\(\sqrt{2\Phi}\) |
richting en grootte uit \(\vec{g}\) |
Bron |
potentiaal (alle massa) |
veld (netto kracht) |
Bolschil binnen |
\(> 0\) |
\(= 0\) |
Bolschil buiten |
$= |
\vec{v}_{netto} |
De formule:
Voor één massa geldt \(|\vec{v}_{netto}| = v_{grav}\), en de formule wordt:
De bekende éénmassa-formule.
2.2 \(\vec{v}_{netto}\) uit het gravitatieveld#
De richting van \(\vec{v}_{netto}\) is de richting van \(\vec{g}\).
De grootte volgt uit \(\vec{g}\) en de gradiënt van \(|\vec{g}|\). Voor één massa geldt:
Algemeen: bereken \(\vec{g}\) en de richtingsafgeleide \(d|\vec{g}|/ds\) langs \(\vec{g}\) (de getijdekracht). Dan:
Alles is berekenbaar uit de massaposities:
\(\vec{g} = \sum G M_i / r_i^2 \cdot \hat{r}_i\)
\(d|\vec{g}|/ds\) numeriek of via de getijdentensor
2.3 De bolschilstelling#
Binnen een uniforme bolschil:
\(\vec{g} = 0\) (Newton: geen nettokracht)
Dus \(\vec{v}_{netto} = 0\) en de kruisterm verdwijnt
Maar \(\Phi \neq 0\), dus \(v_{grav} > 0\)
De formule reduceert tot:
Voor een stilstaand object: \(c_{local}^2 = c^2 - v_{grav}^2\). Isotrope verlaging van \(c_{local}\), geen richtingseffect.
Buiten de bolschil werkt die als een puntmassa: \(|\vec{v}_{netto}| = v_{grav}\), en de oorspronkelijke formule geldt.
| | \(v_{grav}\) | \(|\vec{v}_{netto}|\) | Effect | |—|—|—|—| | Buiten schil | \(\sqrt{2GM/r}\) | \(= v_{grav}\) | Anisotroop (= puntmassa) | | Binnen schil | \(\sqrt{2GM/R}\) | \(0\) | Isotroop (geen kracht) |
2.4 Fysische interpretatie: de rivier#
De analogie met een rivier:
Eén massa: een rivier met een duidelijke stroomrichting. \(v_{grav} = |\vec{v}_{netto}|\): alle energie zit in de coherente stroom.
Binnen een bolschil: water onder druk vanuit alle richtingen. Dezelfde energie (\(v_{grav} > 0\)), maar geen netto stroming (\(\vec{v}_{netto} = 0\)).
Twee gelijke massa’s, middenpunt: twee rivieren die frontaal botsen. De stromen cancelen (\(\vec{v}_{netto} = 0\)), maar de energie blijft (\(v_{grav} > 0\)).
De kruisterm \(2\,(\vec{v}_{ruimte} \cdot \vec{v}_{netto})\) is het effect van de netto stroming op een bewegend object. In stilstaand water (isotroop) is er geen richtingseffect — alleen de druk (tijddilatatie en isotrope ruimte-uitrekking).
3 Visualisaties#
def plot_decomposition(mg, title, x_range, y_range, N=200):
"""Plot v_grav, |v_netto|, en het verschil (isotroop deel) voor een MultiGravity configuratie."""
x = np.linspace(x_range[0], x_range[1], N)
y = np.linspace(y_range[0], y_range[1], N)
data = mg.decompose_grid(x, y)
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
fig.suptitle(title, fontsize=16, fontweight='bold')
v_grav_c = data['v_grav'] / C
v_netto_c = data['v_stroom'] / C # v_stroom in code = v_netto
v_iso_c = data['v_iso'] / C
c_local = np.sqrt(np.maximum(0, 1 - (data['v_grav']/C)**2))
extent = [x_range[0], x_range[1], y_range[0], y_range[1]]
# Panel 1: v_grav/c (scalair)
ax = axes[0, 0]
im = ax.imshow(v_grav_c, origin='lower', extent=extent, cmap='hot', aspect='equal')
plt.colorbar(im, ax=ax, label='$v_{grav}/c$')
ax.set_title('$v_{grav}/c$ (scalair, totale energie)')
# Panel 2: |v_netto|/c met richtingspijlen (lengte = grootte)
ax = axes[0, 1]
im = ax.imshow(v_netto_c, origin='lower', extent=extent, cmap='Blues', aspect='equal')
plt.colorbar(im, ax=ax, label='$|v_{netto}|/c$')
skip = max(1, N // 20)
X, Y = np.meshgrid(x, y)
sx = data['stroom_x']
sy = data['stroom_y']
ax.quiver(X[::skip, ::skip], Y[::skip, ::skip],
sx[::skip, ::skip], sy[::skip, ::skip],
color='navy', alpha=0.6)
ax.set_title('$|\\vec{v}_{netto}|/c$ (vector, netto stroming)')
# Panel 3: isotroop deel = v_grav² - |v_netto|²
ax = axes[1, 0]
im = ax.imshow(v_iso_c, origin='lower', extent=extent, cmap='Greens', aspect='equal')
plt.colorbar(im, ax=ax, label='isotroop deel / $c$')
ax.set_title('$\\sqrt{v_{grav}^2 - |v_{netto}|^2}\\;/\\;c$ (isotroop deel)')
# Panel 4: c_local/c
ax = axes[1, 1]
im = ax.imshow(c_local, origin='lower', extent=extent, cmap='coolwarm', aspect='equal')
plt.colorbar(im, ax=ax, label='$c_{local}/c$')
ax.set_title('$c_{local}/c$')
for ax in axes.flat:
for pos in mg.positions:
ax.plot(pos[0], pos[1], 'ko', markersize=8)
ax.set_xlabel('x [m]')
ax.set_ylabel('y [m]')
plt.tight_layout()
return fig
def plot_profile(mg, title, x_range, N=2000, exclude_radius=0):
"""1D-profiel langs de x-as."""
x = np.linspace(x_range[0], x_range[1], N)
v_grav = np.full(N, np.nan)
v_netto = np.full(N, np.nan)
v_iso = np.full(N, np.nan)
g_x = np.full(N, np.nan)
for i, xi in enumerate(x):
too_close = any(abs(xi - pos[0]) < exclude_radius for pos in mg.positions)
if too_close:
continue
d = mg.decompose(xi, 0)
v_grav[i] = math.sqrt(d['v_grav2']) / C
v_netto[i] = math.sqrt(d['v_stroom2']) / C
v_iso[i] = math.sqrt(d['v_iso2']) / C
g_x[i] = d['g_field'][0]
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8), sharex=True)
fig.suptitle(title, fontsize=14, fontweight='bold')
ax1.plot(x, v_grav, 'r-', lw=2, label='$v_{grav}/c$ (scalair)')
ax1.plot(x, v_netto, 'b-', lw=2, label='$|v_{netto}|/c$ (gericht)')
ax1.plot(x, v_iso, 'g--', lw=2, label='isotroop deel / $c$')
ax1.set_ylabel('v / c')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)
for pos in mg.positions:
ax1.axvline(pos[0], color='gray', ls='--', alpha=0.5)
ax2.axvline(pos[0], color='gray', ls='--', alpha=0.5)
ax2.plot(x, g_x, 'k-', lw=2)
ax2.set_ylabel('$g_x$ [m/s²]')
ax2.set_xlabel('x [m]')
ax2.axhline(0, color='gray', lw=0.5)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
return fig
3.1 Eén massa (referentie)#
Voor één massa geldt \(|\vec{v}_{netto}| = v_{grav}\) overal. Het isotrope deel is nul: alle gravitatie-energie zit in de coherente stroom.
M = 10 * M_SUN
rs = 2 * G * M / C**2
mg1 = MultiGravity([M], [(0, 0)])
print(f'Eén massa: 10 M_zon, rs = {rs:.0f} m')
# Check op r = 10·rs
r_test = 10 * rs
d = mg1.decompose(r_test, 0)
print(f'\nBij r = 10·rs:')
print(f' v_grav/c = {math.sqrt(d["v_grav2"])/C:.6f}')
print(f' |v_netto|/c = {math.sqrt(d["v_stroom2"])/C:.6f}')
print(f' isotroop/c = {math.sqrt(d["v_iso2"])/C:.6f}')
print(f' c_local/c = {mg1.c_local(r_test, 0)/C:.6f}')
print(f' Verwacht: {math.sqrt(1 - 1/10):.6f}')
Eén massa: 10 M_zon, rs = 29541 m
Bij r = 10·rs:
v_grav/c = 0.316228
|v_netto|/c = 0.316228
isotroop/c = 0.000001
c_local/c = 0.948683
Verwacht: 0.948683
L = 30 * rs
fig = plot_decomposition(mg1, 'Eén massa (10 M☉): $|v_{netto}| = v_{grav}$',
(-L, L), (-L, L), N=150)
plt.show()
3.2 Twee gelijke massa’s#
Op het middenpunt cancelen de velden: \(\vec{g} = 0\), dus \(\vec{v}_{netto} = 0\). Alle gravitatie-energie is isotroop — als water onder druk vanuit twee kanten.
Dicht bij één massa domineert die massa en \(|\vec{v}_{netto}| \approx v_{grav}\).
M_each = 5 * M_SUN
rs_each = 2 * G * M_each / C**2
sep = 100 * rs_each
mg2 = MultiGravity([M_each, M_each], [(-sep, 0), (sep, 0)])
print(f'Twee gelijke massa\'s: elk {M_each/M_SUN:.0f} M_zon')
print(f'Afstand: {2*sep/rs_each:.0f} rs')
# Middenpunt
d_mid = mg2.decompose(0, 0)
print(f'\nMiddenpunt:')
print(f' v_grav/c = {math.sqrt(d_mid["v_grav2"])/C:.6f}')
print(f' |v_netto|/c = {math.sqrt(d_mid["v_stroom2"])/C:.6f} (verwacht: 0)')
print(f' isotroop/c = {math.sqrt(d_mid["v_iso2"])/C:.6f}')
# Dicht bij massa 1
r_near = -sep + 10*rs_each
d_near = mg2.decompose(r_near, 0)
print(f'\nDicht bij massa 1:')
print(f' v_grav/c = {math.sqrt(d_near["v_grav2"])/C:.6f}')
print(f' |v_netto|/c = {math.sqrt(d_near["v_stroom2"])/C:.6f}')
print(f' isotroop/c = {math.sqrt(d_near["v_iso2"])/C:.6f}')
Twee gelijke massa's: elk 5 M_zon
Afstand: 200 rs
Middenpunt:
v_grav/c = 0.141421
|v_netto|/c = 0.000000 (verwacht: 0)
isotroop/c = 0.141421
Dicht bij massa 1:
v_grav/c = 0.324443
|v_netto|/c = 0.315329
isotroop/c = 0.076360
L = 2.5 * sep
fig = plot_decomposition(mg2, 'Twee gelijke massa\'s (elk 5 M☉)',
(-L, L), (-L, L), N=200)
plt.show()
fig = plot_profile(mg2, 'Profiel: twee gelijke massa\'s langs de x-as',
(-2.5*sep, 2.5*sep), N=2000, exclude_radius=5*rs_each)
plt.show()
3.3 Twee ongelijke massa’s#
Met ongelijke massa’s verschuift het punt \(\vec{g} = 0\) (het L1 Lagrange-punt) naar de lichtere massa. Op dat punt is \(\vec{v}_{netto} = 0\) en alle gravitatie-energie is isotroop.
M1, M2 = 8 * M_SUN, 2 * M_SUN
rs1 = 2 * G * M1 / C**2
sep3 = 100 * rs1
mg3 = MultiGravity([M1, M2], [(-sep3, 0), (sep3, 0)])
print(f'Massa 1: {M1/M_SUN:.0f} M_zon (links), Massa 2: {M2/M_SUN:.0f} M_zon (rechts)')
Massa 1: 8 M_zon (links), Massa 2: 2 M_zon (rechts)
L = 2.5 * sep3
fig = plot_decomposition(mg3, 'Twee ongelijke massa\'s (8 + 2 M☉)',
(-L, L), (-L, L), N=200)
plt.show()
fig = plot_profile(mg3, 'Profiel: twee ongelijke massa\'s langs de x-as',
(-2.5*sep3, 2.5*sep3), N=2000, exclude_radius=5*rs1)
plt.show()
3.4 Drie massa’s in een driehoek#
Drie massa’s in een gelijkzijdige driehoek. In het zwaartepunt cancelen de velden bij gelijke massa’s (\(\vec{g} = 0\), \(\vec{v}_{netto} = 0\)). Bij ongelijke massa’s verschuift het punt \(\vec{g} = 0\) naar de zwaardere massa’s.
# Drie gelijke massa's in een gelijkzijdige driehoek
M3 = 5 * M_SUN
rs3 = 2 * G * M3 / C**2
a = 100 * rs3 # zijde van de driehoek
h = a * math.sqrt(3) / 2 # hoogte
# Posities: gelijkzijdige driehoek gecentreerd op oorsprong
p1 = (0, 2*h/3)
p2 = (-a/2, -h/3)
p3 = (a/2, -h/3)
mg_tri = MultiGravity([M3, M3, M3], [p1, p2, p3])
print(f'Drie gelijke massa\'s: elk {M3/M_SUN:.0f} M_zon')
print(f'Driehoekszijde: {a/rs3:.0f} rs')
# Zwaartepunt (0, 0)
d_center = mg_tri.decompose(0, 0)
print(f'\nZwaartepunt (0, 0):')
print(f' v_grav/c = {math.sqrt(d_center["v_grav2"])/C:.6f}')
print(f' |v_netto|/c = {math.sqrt(d_center["v_stroom2"])/C:.6f} (verwacht: ~0)')
print(f' g_mag = {d_center["g_mag"]:.4e} m/s²')
Drie gelijke massa's: elk 5 M_zon
Driehoekszijde: 100 rs
Zwaartepunt (0, 0):
v_grav/c = 0.227951
|v_netto|/c = 0.000000 (verwacht: ~0)
g_mag = 4.7684e-07 m/s²
L = 2 * a
fig = plot_decomposition(mg_tri, 'Drie gelijke massa\'s (elk 5 M☉) in driehoek',
(-L, L), (-L, L), N=200)
plt.show()
# Drie ongelijke massa's: 7, 2, 1 zonsmassa's
M_a, M_b, M_c = 7 * M_SUN, 2 * M_SUN, 1 * M_SUN
mg_tri2 = MultiGravity([M_a, M_b, M_c], [p1, p2, p3])
print(f'Drie ongelijke massa\'s: {M_a/M_SUN:.0f}, {M_b/M_SUN:.0f}, {M_c/M_SUN:.0f} M_zon')
fig = plot_decomposition(mg_tri2, 'Drie ongelijke massa\'s (7 + 2 + 1 M☉) in driehoek',
(-L, L), (-L, L), N=200)
plt.show()
Drie ongelijke massa's: 7, 2, 1 M_zon
3.5 Bolschil#
Een dunne bolschil van massa \(M\) bij straal \(R\). Binnen de schil: \(\vec{g} = 0\), dus \(\vec{v}_{netto} = 0\) — puur isotroop. Buiten: identiek aan een puntmassa.
from ort_core import SphericalMassDistribution
# Dunne bolschil: massa 10 M_zon bij straal R
M_shell = 10 * M_SUN
R_shell = 100 * (2 * G * M_shell / C**2) # 100 rs
thickness = R_shell * 0.02 # 2% dikte
def rho_shell(r):
"""Dunne schil bij R_shell met dikte."""
if abs(r - R_shell) < thickness / 2:
return M_shell / (4 * math.pi * R_shell**2 * thickness)
return 0.0
shell = SphericalMassDistribution(rho_shell, R_shell * 3, N_shells=2000)
print(f'Bolschil: M = {M_shell/M_SUN:.0f} M_zon, R = {R_shell:.3e} m')
print(f'M_totaal (check): {shell.M_total/M_SUN:.2f} M_zon')
# Profiel
r_arr = np.linspace(R_shell * 0.01, R_shell * 2.5, 500)
v_grav = [shell.decompose(r)['v_grav'] / C for r in r_arr]
v_netto = [shell.decompose(r)['v_netto'] / C for r in r_arr]
v_iso = [shell.decompose(r)['v_iso'] / C for r in r_arr]
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(r_arr / R_shell, v_grav, 'r-', lw=2, label='$v_{grav}/c$ (scalair)')
ax.plot(r_arr / R_shell, v_netto, 'b-', lw=2, label='$|v_{netto}|/c$ (gericht)')
ax.plot(r_arr / R_shell, v_iso, 'g--', lw=2, label='isotroop deel / $c$')
ax.axvline(1.0, color='gray', ls='--', alpha=0.5, label='schilpositie')
ax.set_xlabel('$r / R_{schil}$')
ax.set_ylabel('$v / c$')
ax.set_title('Bolschil: $v_{grav}$, $|v_{netto}|$ en isotroop deel')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print('\nBinnen de schil (r = 0.5 R):')
d_in = shell.decompose(0.5 * R_shell)
print(f' v_grav/c = {d_in["v_grav"]/C:.6f}')
print(f' |v_netto|/c = {d_in["v_netto"]/C:.6f} (verwacht: ~0)')
print(f' isotroop/c = {d_in["v_iso"]/C:.6f}')
print('\nBuiten de schil (r = 2 R):')
d_out = shell.decompose(2.0 * R_shell)
print(f' v_grav/c = {d_out["v_grav"]/C:.6f}')
print(f' |v_netto|/c = {d_out["v_netto"]/C:.6f}')
print(f' isotroop/c = {d_out["v_iso"]/C:.6f} (verwacht: ~0)')
Bolschil: M = 10 M_zon, R = 2.954e+06 m
M_totaal (check): 9.75 M_zon
Binnen de schil (r = 0.5 R):
v_grav/c = 0.098730
|v_netto|/c = 0.000000 (verwacht: ~0)
isotroop/c = 0.098730
Buiten de schil (r = 2 R):
v_grav/c = 0.069805
|v_netto|/c = 0.069805
isotroop/c = 0.000000 (verwacht: ~0)
# 2D doorsnede bolschil
N2d = 200
L = 2.5
xn = np.linspace(-L, L, N2d)
yn = np.linspace(-L, L, N2d)
Xn, Yn = np.meshgrid(xn, yn)
Rn = np.sqrt(Xn**2 + Yn**2)
v_grav_2d = np.zeros_like(Rn)
v_netto_2d = np.zeros_like(Rn)
v_iso_2d = np.zeros_like(Rn)
c_local_2d = np.zeros_like(Rn)
for j in range(N2d):
for i in range(N2d):
r = max(Rn[j, i], 0.02) * R_shell
d = shell.decompose(r)
v_grav_2d[j, i] = d['v_grav'] / C
v_netto_2d[j, i] = d['v_netto'] / C
v_iso_2d[j, i] = d['v_iso'] / C
c_local_2d[j, i] = d['c_local'] / C
extent = [-L, L, -L, L]
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
fig.suptitle('Bolschil (10 M☉): 2D doorsnede', fontsize=16, fontweight='bold')
vmax = 0.15
ax = axes[0, 0]
im = ax.imshow(v_grav_2d, origin='lower', extent=extent, cmap='hot', aspect='equal', vmin=0, vmax=vmax)
plt.colorbar(im, ax=ax, label='{grav}/c$')
ax.set_title('{grav}/c$ (scalair)')
ax = axes[0, 1]
im = ax.imshow(v_netto_2d, origin='lower', extent=extent, cmap='Blues', aspect='equal', vmin=0, vmax=vmax)
plt.colorbar(im, ax=ax, label='$|v_{netto}|/c$')
ax.set_title('$|v_{netto}|/c$ (gericht — 0 binnen schil)')
skip = max(1, N2d // 15)
for j in range(skip//2, N2d, skip):
for i in range(skip//2, N2d, skip):
rn = Rn[j, i]
vn = v_netto_2d[j, i]
if rn > 0.15 and vn > 0.005:
ux = Xn[j,i] / rn
uy = Yn[j,i] / rn
scale = vn / vmax * 0.25
ax.arrow(Xn[j,i], Yn[j,i], ux*scale, uy*scale,
head_width=0.06, head_length=0.03, fc='navy', ec='navy', alpha=0.5)
ax = axes[1, 0]
im = ax.imshow(v_iso_2d, origin='lower', extent=extent, cmap='Greens', aspect='equal', vmin=0, vmax=vmax)
plt.colorbar(im, ax=ax, label='isotroop/$')
ax.set_title('isotroop/$ (alleen binnen schil)')
cl_min = np.min(c_local_2d)
ax = axes[1, 1]
im = ax.imshow(c_local_2d, origin='lower', extent=extent, cmap='coolwarm', aspect='equal',
vmin=cl_min - 0.001, vmax=1.0)
plt.colorbar(im, ax=ax, label='{local}/c$')
ax.set_title('{local}/c$')
for ax in axes.flat:
circle = Circle((0, 0), 1.0, fill=False, color='white' if ax == axes[0,0] else 'gray',
ls='--', lw=2)
ax.add_patch(circle)
ax.set_xlabel(' / R_{schil}$')
ax.set_ylabel(' / R_{schil}$')
plt.tight_layout()
plt.show()
3.6 Gasbol met dichtheidsprofiel#
Een uniforme gasbol: constante dichtheid \(\rho\) tot straal \(R\). Binnen de bol neemt \(M_{in}\) toe met \(r^3\), dus \(g \propto r\) (lineair). \(v_{netto}\) groeit naar buiten; \(v_{iso}\) komt van de massa buiten \(r\).
# Uniforme gasbol: constante dichtheid tot R
M_gas = 10 * M_SUN
R_gas = 100 * (2 * G * M_gas / C**2)
rho_0 = M_gas / (4/3 * math.pi * R_gas**3)
def rho_uniform(r):
return rho_0 if r <= R_gas else 0.0
gasbol = SphericalMassDistribution(rho_uniform, R_gas * 3, N_shells=2000)
print(f'Uniforme gasbol: M = {M_gas/M_SUN:.0f} M_zon, R = {R_gas:.3e} m')
print(f'Dichtheid: {rho_0:.2e} kg/m³')
print(f'M_totaal (check): {gasbol.M_total/M_SUN:.2f} M_zon')
# Profiel
r_arr = np.linspace(R_gas * 0.01, R_gas * 2.5, 500)
v_grav = [gasbol.decompose(r)['v_grav'] / C for r in r_arr]
v_netto = [gasbol.decompose(r)['v_netto'] / C for r in r_arr]
v_iso = [gasbol.decompose(r)['v_iso'] / C for r in r_arr]
c_local = [gasbol.decompose(r)['c_local'] / C for r in r_arr]
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10), sharex=True)
ax1.plot(r_arr / R_gas, v_grav, 'r-', lw=2, label='$v_{grav}/c$ (scalair)')
ax1.plot(r_arr / R_gas, v_netto, 'b-', lw=2, label='$|v_{netto}|/c$ (gericht)')
ax1.plot(r_arr / R_gas, v_iso, 'g--', lw=2, label='isotroop deel / $c$')
ax1.axvline(1.0, color='gray', ls='--', alpha=0.5, label='rand gasbol')
ax1.set_ylabel('$v / c$')
ax1.set_title('Uniforme gasbol: decompositie')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)
ax2.plot(r_arr / R_gas, c_local, 'purple', lw=2)
ax2.axvline(1.0, color='gray', ls='--', alpha=0.5)
ax2.set_xlabel('$r / R_{gasbol}$')
ax2.set_ylabel('$c_{local} / c$')
ax2.set_title('$c_{local}$ profiel')
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Waarden op sleutelposities
for label, r in [('Centrum', R_gas*0.01), ('r = 0.5R', R_gas*0.5),
('Rand', R_gas*0.99), ('r = 2R', R_gas*2)]:
d = gasbol.decompose(r)
print(f'{label:>10}: v_grav={d["v_grav"]/C:.4f}c '
f'|v_netto|={d["v_netto"]/C:.4f}c '
f'iso={d["v_iso"]/C:.4f}c '
f'c_local={d["c_local"]/C:.6f}c')
Uniforme gasbol: M = 10 M_zon, R = 2.954e+06 m
Dichtheid: 1.84e+11 kg/m³
M_totaal (check): 10.02 M_zon
Centrum: v_grav=0.1225c |v_netto|=0.0011c iso=0.1225c c_local=0.992464c
r = 0.5R: v_grav=0.1173c |v_netto|=0.0499c iso=0.1062c c_local=0.993094c
Rand: v_grav=0.1006c |v_netto|=0.0990c iso=0.0177c c_local=0.994930c
r = 2R: v_grav=0.0708c |v_netto|=0.0708c iso=0.0000c c_local=0.997493c
# 2D doorsnede gasbol
N2d = 200
L = 2.5
xn = np.linspace(-L, L, N2d)
yn = np.linspace(-L, L, N2d)
Xn, Yn = np.meshgrid(xn, yn)
Rn = np.sqrt(Xn**2 + Yn**2)
v_grav_2d = np.zeros_like(Rn)
v_netto_2d = np.zeros_like(Rn)
v_iso_2d = np.zeros_like(Rn)
c_local_2d = np.zeros_like(Rn)
for j in range(N2d):
for i in range(N2d):
r = max(Rn[j, i], 0.02) * R_gas
d = gasbol.decompose(r)
v_grav_2d[j, i] = d['v_grav'] / C
v_netto_2d[j, i] = d['v_netto'] / C
v_iso_2d[j, i] = d['v_iso'] / C
c_local_2d[j, i] = d['c_local'] / C
extent = [-L, L, -L, L]
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
fig.suptitle('Uniforme gasbol (10 M☉): 2D doorsnede', fontsize=16, fontweight='bold')
vmax = 0.15
ax = axes[0, 0]
im = ax.imshow(v_grav_2d, origin='lower', extent=extent, cmap='hot', aspect='equal', vmin=0, vmax=vmax)
plt.colorbar(im, ax=ax, label='{grav}/c$')
ax.set_title('{grav}/c$ (scalair)')
ax = axes[0, 1]
im = ax.imshow(v_netto_2d, origin='lower', extent=extent, cmap='Blues', aspect='equal', vmin=0, vmax=vmax)
plt.colorbar(im, ax=ax, label='$|v_{netto}|/c$')
ax.set_title('$|v_{netto}|/c$ (gericht)')
skip = max(1, N2d // 15)
for j in range(skip//2, N2d, skip):
for i in range(skip//2, N2d, skip):
rn = Rn[j, i]
vn = v_netto_2d[j, i]
if rn > 0.15 and vn > 0.003:
ux = Xn[j,i] / rn
uy = Yn[j,i] / rn
scale = vn / vmax * 0.25
ax.arrow(Xn[j,i], Yn[j,i], ux*scale, uy*scale,
head_width=0.06, head_length=0.03, fc='navy', ec='navy', alpha=0.5)
ax = axes[1, 0]
im = ax.imshow(v_iso_2d, origin='lower', extent=extent, cmap='Greens', aspect='equal', vmin=0, vmax=vmax)
plt.colorbar(im, ax=ax, label='isotroop/$')
ax.set_title('isotroop/$')
cl_min = np.min(c_local_2d)
ax = axes[1, 1]
im = ax.imshow(c_local_2d, origin='lower', extent=extent, cmap='coolwarm', aspect='equal',
vmin=cl_min - 0.001, vmax=1.0)
plt.colorbar(im, ax=ax, label='{local}/c$')
ax.set_title('{local}/c$')
for ax in axes.flat:
circle = Circle((0, 0), 1.0, fill=False, color='gray', ls='--', lw=2)
ax.add_patch(circle)
ax.set_xlabel(' / R_{gasbol}$')
ax.set_ylabel(' / R_{gasbol}$')
plt.tight_layout()
plt.show()
4 Samenvatting#
De twee grootheden#
\(v_{grav}\) (scalair) |
\(\vec{v}_{netto}\) (vector) |
|
|---|---|---|
Definitie |
\(\sqrt{2\Phi}\) |
richting en grootte uit \(\vec{g}\) |
Berekening |
\(\sum 2GM_i/r_i\) |
\(\vec{g} = \sum GM_i/r_i^2 \cdot \hat{r}_i\) |
Bolschil binnen |
\(> 0\) |
\(= 0\) |
Bolschil buiten |
$= |
\vec{v}_{netto} |
Eén massa |
$= |
\vec{v}_{netto} |
De formule#
De grootte van \(\vec{v}_{netto}\)#
waar \(d|\vec{g}|/ds\) de richtingsafgeleide van \(|\vec{g}|\) langs \(\vec{g}\) is.
Speciale gevallen#
Eén massa: \(|\vec{v}_{netto}| = v_{grav}\) → \(c^2 = v_{tijd}^2 + |\vec{v}_{ruimte} + \vec{v}_{grav}|^2\)
Binnen bolschil: \(\vec{v}_{netto} = 0\) → \(c^2 = v_{tijd}^2 + |\vec{v}_{ruimte}|^2 + v_{grav}^2\) (isotroop)
Buiten bolschil: identiek aan puntmassa in het centrum