Почему мой фотон не движется по орбите с правильным радиусом? ⇐ Python

[Anonymous]

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


Однако, если изначально проецировать из x = 3GM/(c^2), это действительно происходит на орбиту правильного радиуса! (Со временем ускользает, так как уязвим к небольшим неточностям)


Все уменьшено на 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)
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

Я ожидаю, что мой оператор нулевой геодезической печати внутри цикла for будет либо 0, либо относительно небольшим целым числом. Например, первый оператор нулевой геодезической печати может иметь значение 0, 16, -8 и т. д. из-за небольшой неточности при большом добавлении чисел с плавающей запятой (величины e17).
У меня есть попытался выполнить отладку, заменив оператор «else» на «elif j == 1» и, посмотрев на следующую итерацию, можно увидеть, что нулевая геодезическая печатает гораздо большее число с плавающей запятой.
Я верю, что подумаю Нулевая геодезическая ошибка покажет, почему мои траектории неверны. Я несколько раз проверял символы и реализацию Кристоффеля и не нашел ничего очевидного.

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