Верна ли моя двухмерная симуляция черной дыры/света? ⇐ Python

[Anonymous]

В Blender я отслеживаю путь некоторых двумерных фотонов, изгибающихся из-за невращающейся стандартной черной дыры, определенной метрикой Шварцшильда. Я устанавливаю свои начальные скорости в терминах r (радиус), фи (угол) и t (время) относительно аффинного параметра лямбда, а затем итеративно обновляю вектор пространства-времени на основе символов Кристоффеля в соответствующей точке. Просто вставьте код в скрипт Python Blender и запустите его.
При проецировании издалека мой фотон не вращается по орбите 3GM/(rc^2). Внутренняя сфера представляет собой черную дыру радиуса
2GM/(rc^2)


Однако при первоначальной проекции из x = 3GM/(rc^2) действительно выходит на орбиту в нужном месте радиус!
[img]https://i.sstatic.net /rUKnZfAk.png[/img]

Все уменьшено на 1e6 для печати в Blender

Код: Выделить всё

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 (near photon sphere)
d_lambda = 1e-4  # Step size for integration

# Black hole properties
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')

# Initialize variables
r = test_point.length * 1e6  # Convert to meters
phi = atan2(test_point.y, test_point.x)  # Planar angle

alpha = 2.55 * pi/8

v_r = -c * cos(alpha)  # Radial component of velocity
v_phi = c * sin(alpha) / r  # Angular component of velocity

# Schwarzschild metric components (2D)
g_tt = -(1 - 2 * G * M / (r * c**2))
g_rr = 1 / (1 - 2 * G * M / (r * c**2))
g_pp = r**2

# Initialize velocity in the time direction (null geodesic condition)
v_t = sqrt(-(g_rr * v_r**2 + g_pp * v_phi**2) / g_tt)

# Debug: Initial null geodesic condition
print(g_tt * v_t**2 + g_rr * v_r**2 + g_pp * v_phi**2)

# 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_phi**2
dv_phi = -Γ_p_rp * v_r * v_phi

# Update velocities
v_r += dv_r * d_lambda
v_phi += dv_phi * d_lambda

# Update positions
r += v_r * d_lambda
phi += v_phi * 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_phi**2)

# Convert spherical to Cartesian coordinates (2D in equatorial plane)
x = (r / 1e6) * cos(phi)
y = (r / 1e6) * sin(phi)

# 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

Я хочу знать, почему мой оператор нулевой геодезической печати внутри цикла for не печатает что-то относительно маленькое. Первый оператор печати может иметь значения 0, 16, -8 и т. д., поскольку при добавлении большого количества чисел с плавающей запятой может возникнуть небольшая неточность.
Я попытался выполнить отладку, заменив оператор "else" на "elif j == 1" и, глядя на следующую итерацию, можно увидеть, что нулевая геодезическая печатает гораздо большее число с плавающей запятой.
Помимо нулевой геодезической, я хочу знать, есть ли у меня орбиты/траектории света верны. Я несколько раз проверял символы и реализацию Кристоффеля и не могу понять, правильно ли это.

Подробнее здесь: https://stackoverflow.com/questions/793 ... on-correct