[Python] Why does my photon not orbit at the correct radius?

In Blender, I am path tracing some 2D photons bending due to a non-rotating standard black hole defined by Schwarzschild metric. I set my initial velocities in terms of r (radius), phi (angle), and t (time) with respect to the affine parameter lambda and then iteratively update the space-time vector based on the Christoffel symbols at the respective point. Simply paste the code into a Blender python script and run it.

When projected from afar, my photon does not orbit at 3GM/(c^2). For reference, the inner sphere is the black hole of radius 2GM/(c^2). At the right side of the loop, the photon enters within the orbit radius and does not spiral to the event horizon!



However, when initially projected from x = 3GM/(c^2), it does indeed go into orbit at the correct radius! (Eventually escaping as it's vulnerable to slight inaccuracies)



Everything scaled down by 1e6 to plot in 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
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


I expect my null geodesic print statement inside the for loop to be either 0 or a relatively small integer. For example, the first null geodesic print statement may be 0, 16, -8, etc. due to a small amount of imprecision with the large float addition (e17 magnitudes).

I have tried to debug by replacing my "else" statement with "elif j == 1" and looking at the very next iteration, it can be seen the null geodesic prints a much larger float.

I believe figuring out the null geodesic error will reveal why my trajectories are incorrect. I have checked my Christoffel symbols and implementation multiple times and can't catch anything obvious.

Continue reading...