I am trying to raytrace an image of a black hole, but given my current difficulties of doing so with the formalism of relativity, I've decided to do my raytracing using Newton's law of gravitation instead:
a = -GM / r^2
This is my code so far:
import numpy as np
from PIL import Image
M = 500.0
max_iterations = 100
dt = 0.1
img_width = 160
img_height = 90
resolution = (img_width, img_height)
sensor_width = 0.035
focal_length = 0.025
camera_origin = np.array([0, 7 * M, 0])
def norm(x):
return (x.dot(x)) ** (1/2)
def normalize(x):
return x / norm(x)
class Camera:
def __init__(self, pos=camera_origin, resolution=resolution, sensor_width=sensor_width, focal_length=focal_length, z_rotation=np.pi / 4):
rot_z = np.array([
[np.cos(z_rotation), -np.sin(z_rotation), 0],
[np.sin(z_rotation), np.cos(z_rotation), 0],
[0, 0, 1]
])
self.rot_matrix = rot_z
self.pos = np.matmul(rot_z, pos)
self.resolution = resolution
self.sensor_width = sensor_width
self.focal_length = focal_length
def ray_vel(self, i, j):
w, h = self.resolution
aspect_ratio = w / h
# The u and v vectors are
# the vectors of the camera
# viewport
u = self.sensor_width
v = -self.sensor_width / aspect_ratio
du = u / w
dv = v / h
# Find the directional vector
# of a ray that passes through
# a given pixel
x = -u/2 + (i + 1) * du
y = self.focal_length
z = -v/2 + (j + 1) * dv
pixel_dir = np.array([x, y, z])
# Rotate the directional vector about
# the z-axis with given amount
pixel_dir = np.matmul(self.rot_matrix, pixel_dir)
# Make these unit vectors
# as speed of light in geometrized
# units is 1
return normalize(pixel_dir)
class Ray:
def __init__(self, origin, direction):
self.origin = origin
self.direction = direction
def compute(self, t):
pos = self.origin
vels = self.direction
# Use Euler's method to
# numerically solve
# the differential equation
# a = -GM/r^2
for i in range(t):
r = norm(pos)
vels -= M / (r ** 2) * dt
pos += vels * dt
return vels
def trace(self, tmax=100):
ray_dir = self.compute(tmax)
v_x, v_y, v_z = normalize(ray_dir)
color_r = int(abs(v_x) * 255)
color_g = int(abs(v_y) * 255)
color_b = int(abs(v_z) * 255)
return (color_r, color_g, color_b)
bh_img = Image.new("RGB", (img_width, img_height))
counter = 0
camera = Camera()
print("Beginning raytracing...")
for i in range(0, img_width):
for j in range(0, img_height):
ray_origin = camera.pos
ray_dir = camera.ray_vel(i, j)
ray = Ray(ray_origin, ray_dir)
color = ray.trace(max_iterations)
bh_img.putpixel((i, j), color)
if counter % 1000 == 0:
print(f"Computed ray {counter} of {img_width * img_height}")
counter += 1
bh_img.save("render_output.jpg")
My code essentially emits rays from the camera, computes their trajectories as they orbit around the black hole, and then finds their final direction, which I have outputted as an RGB color for debugging purposes. What is bizarre to me is that regardless of the values of M I use, the raytraced image looked exactly as if I just outputted the initial ray direction. There is seemingly no change in ray direction as the rays are traced around the black hole. Is there any reason why this is?
