Bei der Projektion aus der Ferne kreist mein Photon nicht mit 3GM/(c^2). Als Referenz: Die innere Kugel ist das Schwarze Loch mit dem Radius 2GM/(c^2). Auf der rechten Seite der Schleife tritt das Photon innerhalb des Umlaufbahnradius ein und bewegt sich nicht spiralförmig zum Ereignishorizont!

Bei der anfänglichen Projektion von x = 3GM/(c^2) funktioniert es jedoch tatsächlich in die Umlaufbahn Radius! (Irgendwann entkommt es, da es anfällig für leichte Ungenauigkeiten ist)

Alles um 1e6 verkleinert, um es in Blender darzustellen
Code: Select all
import bpy
from math import sin, cos, pi, sqrt, atan2
from mathutils import Vector
# Constants
M = 9e32 # Mass of the black hole in kg
c = 299792458 # Speed of light in m/s
G = 6.6743e-11 # Gravitational constant in N m^2/kg^2
# Simulation parameters
curve_num = 4000 # Resolution of the curve
#test_point = Vector((0, 4, 0)) # Initial position
test_point = Vector((3*G*M/(1e6 * c**2), 0, 0))
alpha = 4 * pi/8 # initial projection angle (at x = 0, y = 4 try 2.55 * pi/8)
d_lambda = 1e-4 # Step size for integration
# Initialize variables
r = test_point.length * 1e6 # Convert to meters
p = atan2(test_point.y, test_point.x) # p = phi = Initial planar angle
# initial radial and angular velocities
v_r = -c * cos(alpha) # Radial component of velocity
v_p = c * sin(alpha) / r # Angular component of velocity
# Schwarzschild metric components (2D)
g_tt = -(1 - 2 * G * M / (r * c**2)) # time
g_rr = 1 / (1 - 2 * G * M / (r * c**2)) # radial
g_pp = r**2 # phi
# Initialize velocity in the time direction
v_t = sqrt(-(g_rr * v_r**2 + g_pp * v_p**2) / g_tt)
# Debug: Initial null geodesic condition
print(g_tt * v_t**2 + g_rr * v_r**2 + g_pp * v_p**2)
# Plot black hole sphere mesh
horizon = 2 * G * M / c**2 # Event horizon radius
bpy.ops.mesh.primitive_uv_sphere_add()
black_hole = bpy.context.object
black_hole.scale = (horizon / 1e6, horizon / 1e6, horizon / 1e6) # Scale for Blender units
# Create a 2D trajectory curve
bpy.ops.mesh.primitive_circle_add(radius=1, vertices=curve_num + 1)
curve = bpy.context.object
verts = curve.data.vertices
# Delete unnecessary vertex to make the trajectory a 2D curve
bpy.ops.object.mode_set(mode='EDIT')
bpy.ops.mesh.select_all(action='DESELECT')
bpy.ops.object.mode_set(mode='OBJECT')
verts[-1].select = True
bpy.ops.object.mode_set(mode='EDIT')
bpy.ops.mesh.delete(type='VERT')
bpy.ops.object.mode_set(mode='OBJECT')
# Main simulation loop
for i, v in enumerate(verts):
if i == 0:
v.co = test_point
else:
if r > horizon + 1000:
# Relevant Christoffel symbols
Γ_r_tt = (G * M / (r**2 * c**2)) * (1 - 2 * G * M / (r * c**2))
Γ_r_rr = -(G * M / (r**2 * c**2)) / (1 - 2 * G * M / (r * c**2))
Γ_r_pp = -r * (1 - 2 * G * M / (r * c**2))
Γ_p_rp = 1 / r
# Get accelerations
dv_r = -Γ_r_tt * v_t**2 - Γ_r_rr * v_r**2 - Γ_r_pp * v_p**2
dv_p = -Γ_p_rp * v_r * v_p
# Update velocities
v_r += dv_r * d_lambda
v_p += dv_p * d_lambda
# Update positions
r += v_r * d_lambda
p += v_p * d_lambda
# Recalculate metric coefficients
g_tt = -(1 - 2 * G * M / (r * c**2))
g_rr = 1 / (1 - 2 * G * M / (r * c**2))
g_pp = r**2
# Debug: Print null geodesic condition
print(g_tt * v_t**2 + g_rr * v_r**2 + g_pp * v_p**2)
# Convert spherical to Cartesian coordinates (2D in equatorial plane)
x = (r / 1e6) * cos(p)
y = (r / 1e6) * sin(p)
# Update vertex position
v.co = Vector((x, y, 0))
# Finalize the curve
bpy.context.view_layer.objects.active = curve
bpy.ops.object.convert(target='CURVE')
curve = bpy.context.object
curve.data.bevel_depth = 0.01
curve.data.fill_mode = 'FULL'
curve.data.bevel_resolution = 12
Das habe ich Ich habe versucht zu debuggen, indem ich meine „else“-Anweisung durch „elif j == 1“ ersetzt habe, und als ich mir die nächste Iteration ansah, konnte ich sehen, dass die Null-Geodäte einen viel größeren Float druckt.
Ich glaube, das herauszufinden Ein geodätischer Nullfehler zeigt, warum meine Flugbahnen falsch sind. Ich habe meine Christoffel-Symbole und die Implementierung mehrmals überprüft und kann nichts Auffälliges erkennen.