Ray tracing black holes

Lately I’ve been studying up on ray tracing, and one of my goals has been to build a nonlinear ray tracer — that is, a ray tracer that works in curved space, for example space that is curved by a nearby black hole. (See the finished source code!)

In order to do this, the path of each ray must be calculated in a stepwise fashion, since we can no longer rely on the geometry of straight lines in our world. With each step taken by the ray, the velocity vector of the ray is updated based on an equation of motion determined by a “force field” present in our space.

This idea has certainly been explored in the past, notably by Riccardo Antonelli, who derived a very clever and simple equation for the force field that guides the motion of the ray in the vicinity of a black hole, namely

$$\vec F(r) = – \frac{3}{2} h^2 \frac{\hat r}{r^5}$$

I decided to use the above equation in my own ray tracer because it’s very efficient computationally (and because I’m not nearly familiar enough with the mathematics of GR to have derived it myself). The equation models a simple Schwarzschild black hole (non-rotating, non-charged) at the origin of our coordinate system. The simplicity of the equation has the tradeoff that the resulting images will be mostly unphysical, meaning that they’re not exactly what a real observer would “see” in the vicinity of the black hole. Instead, the images must be interpreted as instantaneous snapshots of how the light bends around the black hole, with no regard for redshifting or distortions relative to the observer’s motion.

Nevertheless, this kind of ray tracing provides some powerful visualizations that help us understand the behavior of light around black holes, and help demystify at least some of the properties of these exotic objects.

My goal is to build on this existing work, and create a ray tracer that is more fully featured, with support for other types of objects in addition to the black hole. I also want it to be more extensible, with the ability to plug in different equations of motion, as well as to build more complex scenes, or even to build scenes algorithmically. So, now that my work on this ray tracer has reached a semi-publishable state, let’s dive into all the things it lets us do.

Accretion disk

The ray tracer supports an accretion disk that is either textured or plain-colored. It also supports multiple disks, at arbitrary radii from the event horizon, albeit restricted to the horizontal plane around the black hole. The collision point of the ray with the disk is calculated by performing a binary search for the exact intersection. If we don’t search for the precise point of intersection, we would see artifacts due to the “resolution” of the steps taken by each ray (notice the jagged edges at the bottom of the disk):

Once the intersection search is implemented, the lines and borders become nice and crisp:

We can also apply different colors to the top and bottom of the disk. Observe that the black hole distorts the disk in a way that makes the bottom (colored in green) appear around the lower semicircle of the photon sphere, even though we’re looking at the disk from above:

Note that the dark black circle is not the event horizon, but is actually the photon sphere. This is because photons that cross into the photon sphere from the outside cannot escape. (Only photons that are emitted outward from inside the photon sphere can be seen by an outside observer.)

If we zoom in on the right edge of the photon sphere, we can see higher-order images of the disk appear around the sphere (second- and even third-order images are visible). These are rays of light that have circled around the photon sphere one or more times, and eventually escaped back to the observer.

And here is the same image with a more realistic-looking accretion disk:

Great! Now that we have the basics out of the way, it’s time to get a little more crazy with ray tracing arbitrary materials around the black hole.

Additional spheres

The ray tracer allows adding an unlimited number of spheres, positioned anywhere (outside the event horizon, that is!) and either textured or plain-colored. Here is a scene with one hundred “stars” randomly positioned in an “orbit” around the black hole (click to view larger versions of the images):

Notice once again how we can see second- and third-order images of the spheres as we get closer to the photon sphere. By the way, here is a similar image of stars around the black hole, but with the curvature effects turned off (as if the black hole did not curve the surrounding space):

And here is a video, generated using the ray tracer, that shows the observer circling around the black hole with stars in its vicinity. Once again, this is not a completely realistic physical picture, since the stars are not really “orbiting” around the black hole, but rather it’s a series of snapshots taken at different angles:

Notice how the spherical stars are distorted around the Einstein ring, as well as how the background sky is affected by the curvature.

Reflective spheres

And finally, the ray tracer supports adding spheres that are perfectly reflective:

All that’s necessary for doing this is to calculate the exact point of impact by the ray on the sphere (again using a binary intersection search) and get the corresponding reflected velocity vector based on the normal vector on the sphere at that point. Here is a similar image, but with a textured accretion disk:

Future work

Eventually I’d like to incorporate more algorithms for different equations of motion for the rays. For example, someone else has encoded a similar algorithm for a Kerr black hole (i.e. a black hole with angular momentum), and there is even a port of it to C# already, which I was able to integrate into my ray tracer easily:

A couple more ideas:

  • There’s no reason the ray tracer couldn’t support different types of shapes besides spheres, or even arbitrary mesh models (e.g. STL files).
  • I’d also like to use this ray tracer to create some more animations or videos, but that will have to be the subject of a future post.
  • Make it run on CUDA?

Pi is wrong! Long live Tau!

At one point or another, we’ve all had a feeling that something is not quite right in the world. It’s a huge relief, therefore, to discover someone else who shares your suspicion. (I’m also surprised that it’s taken me this long to stumble on this!)

It has always baffled me why we define \(\pi\) to be the ratio of the circumference of a circle to its diameter, when it should clearly be the ratio of the circumference to its radius. This would make \(\pi\) become the constant 6.2831853…, or 2 times the current definition of \(\pi\).

Why should we do this? And what effect would this have?

Well, for starters, this would remove an unnecessary factor of 2 from a vast number of equations in modern physics and engineering.

Most importantly, however, this would greatly improve the intuitive significance of \(\pi\) for students of math and physics. \(\pi\) is supposed to be the “circle constant,” a constant that embodies a very deep relationship between angles, radii, arc lengths, and periodic functions.

The definition of a circle is the set of points in a plane that are a certain distance (the radius) from the center. The circumference of the circle is the arc length that these points trace out. The circle constant, therefore, should be the ratio of the circumference to the radius.

To avoid confusion, we’ll use the symbol tau (\(\tau\)) to be our new circle constant (as advocated by Michael Hartl, from the Greek τόρνος, meaning “turn”), and make it equal to 6.283…, or \(2\pi\).

In high school trigonometry class, students are required to make the painful transition from degrees to radians. And what’s the definition of a radian? It’s the ratio of the length of an arc (a partial circumference) to its radius! Our intuition should tell us that the ratio of a full circumference to the radius should be the circle constant.

Instead, students are taught that a full rotation is \(2\pi\) radians, and that the sine and cosine functions have a period of \(2\pi\). This is intuitively clunky and fails to illustrate the true beauty of the circle constant that \(\pi\) is supposed to be. This is surely part of the reason that so many students fail to grasp these relationships and end up hating mathematics. A full rotation should be τ radians! The period of the sine and cosine functions should be \(\tau\)!

But… wouldn’t we have to rewrite all of our textbooks and scientific papers that make use of \(\pi\)?

Yes, we would. And, in doing so, we would make them much easier to understand! You can read the Tau Manifesto website to see examples of the beautiful simplifications that \(\tau\) would bring to mathematics, so I won’t repeat them here. You can also read the original opinion piece by Bob Palais that explores this subject.

It’s not particularly surprising that the ancient Greeks used the diameter of a circle (instead of the radius) in their definition of \(\pi\), since the diameter is easier to measure, and also because they couldn’t have foreseen the ubiquity of this constant in virtually all sciences.

However, it’s a little unfortunate that someone like Euler, Leibniz, or Bernoulli didn’t pave the way for redefining \(\pi\) to be 6.283…, thus missing the opportunity to simplify mathematics for generations to come.

Aside from all the aesthetic improvements this would bring, considering how vitally important it is for more of our high school students (and beyond) to understand and appreciate mathematics, we need all the “optimizations” we can get to make mathematics more palatable for them. This surely has to be an optimization to consider seriously!

From now on, I’m a firm believer in tauism! Are you?

The Math Book: Get it Now!

Cliff Pickover, the prolific author of more than forty popular science and mathematics books, has outdone himself with his latest compilation: The Math Book. This is a collection of 250 “milestones” of mathematics throughout history, complete with breathtaking glossy color illustrations for each entry (a first for his books), as well as insightful descriptions that explain the history and the significance of each of these marvels of mathematics.

This book is especially significant in one other way: it contains my artwork! The book’s entry on Knight’s Tours (p. 186) familiarizes the reader with the history of this problem, dating all the way back to Euler in 1759. And, alongside the article, Pickover displays a 30×30 knight’s tour that was solved by my neural network knight’s tour implementation. For the picture in the book, I used a modified version of the program that generated a sufficiently hi-res image. That particular knight’s tour took about 3 days for my computer to generate.

I’m deeply grateful to have one of my creations published in a book by someone as influential as Cliff Pickover. Of course, it’s all of the 250 entries in the book that make it an incredibly fascinating stroll through the history of mathematics. As mentioned elsewhere, this book definitely has bestseller potential, and could easily be one of Pickover’s best works. Buy the book now!

Binomial Coefficients and Stirling Numbers in C#

As long as I’m on a roll with my posts on number theory in C#, I thought I’d briefly discuss how to generate binomial coefficients, and then move on to Stirling numbers of the second kind, both of which are extremely useful in combinatorics and finite calculus.

Some prerequisites for calculating binomial coefficients include a standard factorial function, nothing special:

Note that I use longs whenever possible because, despite the performance hit from 64-bit operations, it’s worth it to be able to work with numbers that are double the magnitude of ints.

We will also need a function for calculating a falling power of a number, which is defined as: $$x^{\underline{n}} = x(x-1)(x-2)\ldots(x-(n-1))$$You’ll see why in a moment.

Notice that the above function handles the case $$n^{\underline{0}}$$ returning 1. However, it does not handle the case of p > n, which would give an incorrect result, but this is not necessary for our purposes.

Binomial Coefficients

Recall that the definition for the binomial coefficient is $${n \choose k} = \frac{n!}{k!(n-k)!}$$ However, using this exact formula to compute binomial coefficients is a bit naive. If we use falling powers (sometimes called falling factorials), the above formula easily reduces to: $$\frac{n^{\underline{k}}}{k!}$$

We can improve the algorithm a bit more by adding the condition:

$${n \choose k} =
\begin{cases}
\frac{n^{\underline{k}}}{k!} \quad \mbox{if } k \leq \lfloor n/2 \rfloor,\\
\frac{n^{\underline{n-k}}}{(n-k)!} \quad \mbox{if } k > \lfloor n/2 \rfloor.
\end{cases}
$$

Not only is this algorithm faster, but it can also handle larger coefficients than the original formula, since neither the falling power nor the factorial ever gets larger than n/2.
The code for this is straightforward:

However, this is still not as optimal as it can be. The most optimal approach would be to accumulate the falling power while dividing by each factor of the factorial in place. This would minimize the chance of overflow errors, and allow for even larger coefficients to be calculated. The disadvantage of this algorithm is the necessary use of floating-point math:

Stirling Numbers

Stirling numbers (of the second kind) are useful for, among other things, enumerating the coefficients of the falling-power expansion of a regular power. For example, how would we express x3 in terms of falling powers of x? That is, how do we arrive at an equation of the form $$x^3 = ax^{\underline 3} + bx^{\underline 2} + cx^{\underline 1}$$ Well, we could just solve the equation directly, but this would get unwieldy for higher powers. A neat way of doing this involves the use of Stirling numbers of the second kind, $$\left\{\begin{matrix} n \\ k \end{matrix}\right\}$$ A useful theorem for computing these numbers is

$$ \left\{\begin{matrix} n \\ k \end{matrix}\right\}=\frac{1}{k!} \sum_{i=0}^{k}(-1)^i{k \choose i}(k-i)^n $$

With the Stirling numbers in hand, we can now obtain the coefficients for falling power expansions:

$$x^m = \sum_{k=0}^{m}\left\{\begin{matrix} m \\ k \end{matrix}\right\}x^{\underline k}$$

Project Euler — Problem 197

I never thought that solving random math problems can be addictive, but Project Euler does exactly that. Not only does it make you flex your math muscles, it also challenges you to take your programming language of choice to its limits. So it’s a total win-win: you brush up on your math skills, and broaden your programming repertoire at the same time.

Once I discovered Project Euler, I couldn’t pull myself away from the computer until I solved as many problems as I could. As of this writing I’ve solved 187 of their 211 problems.

Problem 197 has to do with finding the nth term of a particular recursively defined sequence: $$u_{n+1} = f(u_n)$$ with $$u_0 = -1, f(x) = \lfloor 2^{30.403243784-x^2}\rfloor \cdot 10^{-9}$$

Of course, as with most Project Euler problems, the value for n is set ridiculously high, presumably to eliminate the possibility of brute-forcing the problem (within the lifetime of the universe).

Fortunately, it takes little more than a superficial examination to see that this problem is actually quite simple in disguise. If we look at the first few terms of the sequence, we can already guess that this sequence appears to be “converging” to a function that oscillates between two values, approximately 0.681 and 1.029 (I won’t give precise numbers, since that would give away the solution).

This means that all we need to do is go far enough into the sequence that the deviation of the oscillations is less than the desired precision asked by the problem (10-9). And it so happens that we don’t need to go out far at all. The sequence actually settles on its two oscillatory values as early as the 1000th term (probably even earlier)! Therefore, the sum of the 1000th and 1001st term will be equivalent to the sum of the 1012th and (1012+1)st term, which is what the problem asks for!

The code to do this is elementary. I accomplished it with a mere 5 lines of C# code. Can you do better?