Simulation of light rays under the influence of a large gravitational source is of great importance in physics and astronomy. The standard techniques to achieve this simulation are either fast but inaccurate, or accurate but slow. We have implemented each of these techniques, as well as an experimental technique which could provide a balance between speed and accuracy.

Method 1: Bent Ray

The first technique we used involves bending the ray at a single point. The bend point is chosen to be the point along the initial ray which is closest to the point mass (the source of gravitational influence), at time \(t_c\). The ray bends towards the point mass, with angle of deflection \(\theta = \frac{4Gm}{c^2 r}\), where \(r\) is the distance to the point mass.

There were not many technical difficulties with this approach. The main issue had to do with discrepancies in time values. The second ray-segment needed to start at time "t = 0" in its own reference frame, but at time \(t_c\) in the global frame. The solution involved simply checking intersections for each ray-segment independently, and then handling the various cases accordingly.

Method 2: Numerical Integration

The second technique we implemented was the most physically accurate. The acceleration due to gravity is given by \(\frac{d^2x}{dt^2} = \frac{Gm}{r^2}\). By integrating the acceleration over time, we can find the velocity function, and integrating the velocity function over time gives us the position function. To estimate these integrals, the position and direction of the photon are updated at discrete time steps using a numerical integration technique such as Euler's method, Euler-Cromer, or Runge-Kutta. One thing we had to keep in mind was that, since a photon cannot "speed up" or "slow down," we only want to concern ourselves with the component of force perpendicular to the current velocity. We also enforced this "constant speed" rule every time we updated the velocity. We implemented a total of four integration techniques, listed below.

One drawback to this approach is that intersections with objects must be checked iteratively.

Method 3: Hyperbolic Ray

The third and final technique we implemented was intended to be a hybrid between the first two methods, achieving a curved path but in constant time. We want a hyperbola defined by \(\frac{y^2}{a^2} - \frac{x^2}{b^2} = 1\), but embedded in a plane with its own coordinate system. We start by finding the same bent ray from the first technique. Given an initial ray with origin \(\vec{o}\) and direction \(\vec{d}\), let \(\vec{c}\) be the point along the initial ray which is closest to \(\vec{p}\), the position of the point mass. Let \(\vec{d_b}\) be the direction of the bend, determined by angle \(\alpha\) as before. Let \(\vec{d_m}\) be the average of \(\vec{d_b}\) and \(-\vec{d}\). Let \(\vec{c}'\) be the point along the initial ray such that \(\vec{d_m}\) points directly towards \(\vec{p}\) when starting at \(\vec{c}'\). This point \(\vec{c}'\) will be the origin of our coordinate system. The x-axis is defined by the unit vector \(\vec{x} = \frac{-\vec{d_b}-\vec{d}}{||\vec{d_b}+\vec{d}||}\), and the y-axis is defined as \(-\vec{d_m}\). These definitions are illustrated in diagrams 1 and 2.

Finally, the parameters \(a\) and \(b\) defining the shape of the hyperbola are calculated. The asymptotes should have slope \(\pm \frac{a}{b}\), so we know \(tan(\alpha) = \frac{a}{b}\). We define \(a = D cos(\theta)\), where \(D\) is the distance from \(\vec{p}\) to \(\vec{c}'\), and \(\theta\) is the angle between \(d\) and \(\vec{p} - \vec{o}\). This way, as the initial ray points more and more towards the point mass, the hyperbola deviates further from this initial ray.

The position of the photon at time \(t\) is given by \(\vec{r}(t) = \vec{c}' - b\cdot sinh(t) \vec{x} - a\cdot cosh(t) \vec{y}\).

Intersecting hyperbolic rays with primitives is tricky, as this requires solving the intersection in implicit form, rather than parametric, and then calculating the time of intersection using inverse cosh or sinh. To intersect with triangles, the approach we took was to first intersect the planes that the triangle and hyperbola are each embedded in. This intersection gives a line, which we then intersect with the hyperbola in its own coordinate system. This simply involves solving a degree 2 polynomial, so it can return at most two solutions. We then convert these points back to world coordinates, and check whether either one lies inside the triangle.

To intersect with bounding boxes, we first check whether the vertices of the box all lie to one side of the hyperbolas native plane. If they do, we return false immediately. Otherwise, we form 12 triangles representing the surface of the box, and then check intersection with each of these using the same process as above. In order to avoid copy pasting the code for this, we created a "Basic Triangle" struct, defined only by its vertices, then placed the intersection testing code there, so that bounding boxes and mesh triangles can each access the method for their own use.

To intersect with spheres, we first check whether the sphere intersects the plane. If it does, we determine the circle of intersection, and find the representation of this circle in the private coordinate system of the hyperbola. We then check intersection between the hyperbola and this circle, which involves solving a degree 4 polynomial.

One thing we discovered was that it made more sense to store the point of intersection directly in the "intersection" object, rather than recalculating it using the time of intersection.

Comparisons

Here is a comparison of the run-times for these images:

• Bentray: 15.3380 seconds
• Euler: 1405.0714 seconds
• Euler-Cromer: 1333.2845 seconds
• 2nd Order Runge-Kutta: 1422.1506 seconds
• 4th Order Runge-Kutta: 1616.7614 seconds
• Hyperbolic Ray: 1065.3778 seconds

The hyperbolic method clearly overestimates the size of the "Einstein-ring" around the black hole, and also underestimates the radius of the black hole itself. However, it does manage to be a little bit faster than any of the numerical integration methods.

As an aside, we also considered the possibility of using a feed-forward neural network to solve the geodesic equation. We did not have enough time for this, and neither of us has taken a class on machine learning. Still, this paper suggests that this technique could provide speed and accuracy. To achieve this, it would probably be necessary to "memoize" a set of network states so that it does not have to retrain for each new camera ray. The output would be a piece-wise linear solution, similar to that provided by the numerical integration techniques, so intersections with objects would still be computationally intensive.

Project Video

Paper

https://drive.google.com/a/berkeley.edu/file/d/0B-aLZs_0_nk7MEdhUE9MVTh5c1U/view?usp=sharing

References

• Similar project from previous semester:
http://grtracer.herokuapp.com/
• Gonsalves, R.J., 2011.
Solar system orbits and runge-kutta algorithms.
• Lagaris, I.E., Likas, A., and Fotiadis, D.~I. 1997.
Artificial neural networks for solving ordinary and partial differential equations.
• Sandvik, A.W., 2016.
Numerical solutions of classical equations of motion.
• Ericson, C., 2004.
Real-Time Collision Detection
• Eigen
http://eigen.tuxfamily.org/dox/

Contributions

Charles wrote the implementations for Bent Ray and the Numerical Integration techniques. Tyler wrote the implementation for Hyperbolic Ray. We both were involved in debugging the code and deciding how we would extend the existing program structure.