1
$\begingroup$

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?

$\endgroup$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.