Почему нулевая геодезическая черная дыра не печатает ноль? Верны ли мои траектории? ⇐ Python

[Anonymous]

Используя matplotlib, я отслеживаю путь некоторых двумерных фотонов, изгибающихся из-за невращающейся стандартной черной дыры, определенной метрикой Шварцшильда. Я устанавливаю свои начальные скорости в терминах r (радиус), фи (угол) и t (время) относительно аффинного параметра лямбда, а затем итеративно обновляю вектор пространства-времени на основе символов Кристоффеля в соответствующей точке. /p>
Основная проблема заключается в том, почему нулевые геодезические операторы не выводят что-то близкое к нулю.
Все масштабируется на 1e6 для построения графика
>

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

import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos, sqrt, atan2, pi

# 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 = 2000  # Resolution of the curve
d_lambda = 1e-4  # Step size for integration

# Black hole parameters
horizon = 2 * G * M / c**2  # Event horizon radius
photon_sphere = 3 * G * M / c**2  # Photon sphere radius

# Test point parameters
num = 300  # Number of photons
x_i = 7.5
y_i = 2.5
y_f = 4

# Initial velocity in Cartesian coordinates
v_x = -c  # Speed of light in negative x-direction
v_y = 0   # No velocity in the y-direction

# Prepare for matplotlib plotting
trajectories = []

for i in range(num):
trajectory = []

# Define initial test point
if num > 1:
test_point = np.array([x_i, y_i + i / (num - 1) * (y_f - y_i), 0])
else:
test_point = np.array([x_i, y_i, 0])

# Convert to polar coordinates
r = np.linalg.norm(test_point) * 1e6  # Radius in meters
p = atan2(test_point[1], test_point[0])  # Initial planar angle

# 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

# null geodesic
print(g_tt * v_t**2 + g_rr * v_r**2 + g_pp * v_p**2)

# Initial velocities
v_r = (v_x * cos(p) + v_y * sin(p))  # Radial velocity
v_p = (-v_x * sin(p) + v_y * cos(p)) / r  # Angular velocity
v_t = sqrt(-(g_rr * v_r**2 + g_pp * v_p**2) / g_tt)

# Integrate geodesics
for j in range(curve_num):
if r > horizon + 10000:
# 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
Γ_t_rt = (G * M / (r**2 * c**2)) / (1 - 2 * G * M / (r * c**2))

# Update change in velocities
dv_r = -Γ_r_tt * v_t**2 - Γ_r_rr * v_r**2 - Γ_r_pp * v_p**2
dv_p = -2 * Γ_p_rp * v_r * v_p
dv_t = -2 * Γ_t_rt * v_r * v_t

# Update velocities
v_r += dv_r * d_lambda
v_p += dv_p * d_lambda
v_t += dv_t * d_lambda

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

# null geodesic
print(g_tt * v_t**2 + g_rr * v_r**2 + g_pp * v_p**2)

# Update positions
r += v_r * d_lambda
p += v_p * d_lambda

# Store Cartesian coordinates
x = (r / 1e6) * cos(p)
y = (r / 1e6) * sin(p)

# Only store points within the -10 to 10 range
if -10 

Подробнее здесь: [url]https://stackoverflow.com/questions/79357401/why-is-black-hole-null-geodesic-not-printing-zero-are-my-trajectories-correct[/url]