This is a slightly polished version the article that I previously published in Russian [here](https://habr.com/ru/articles/732926/).
In this article, I'll show you how to render a physically accurate black hole, in full compliance with General Relativity.
In order to do that, we'll derive motion equations for light rays, implement Runge-Kutta integrator in GLSL, and finally combine the two into a fragment shader that traces the light rays emitted from the camera into the past.
# Preparation: Metrics, Connection, Geodesics
![[colored_background.png|Gravitational lensing]]
Light travels in a straight line. But what does straight really mean? Something described by a linear equation? Well, no. If you take a polar coordinate system, the equation for a straight line is anything but linear. And if the space itself isn't flat (like this torus), there's no choice of coordinates where straight lines are linear.
![[thorus_geo.png|Straight line on a torus]]
So, how do we consider a curve "straight"? Well, let's take a vector that points along the curve at some point, and then move it following the curve. Will it still point along it? If yes for any point, such a curve is called _geodesic_, and this is the closest it gets to the straight line on some arbitrary manifold.
Deriving equations for geodesics on a 4-dimensional manifold is quite a laborious task. Fortunately, we don't have to do it by hand. Let the machines do laborious tasks. Take the python package [sympy](https://docs.sympy.org/latest/index.html).
```python
from sympy import *
from sympy.diffgeom import *
import numpy as np
```
The [sympy.diffgeom](https://docs.sympy.org/latest/modules/diffgeom.html) module is somewhat inconvenient and clunky, but it does the job of deriving equations on manifolds described by an [atlas](https://en.wikipedia.org/wiki/Atlas_(topology)) reasonably well. Our atlas will have just one patch, but with a couple of coordinate frames: one spherical (for convenience when working with the Schwarzschild metric), and one Cartesian - for our future camera.
Let's define the coordinate symbols:
```python
t_, x_, y_, z_ = symbols("t x y z")
rho_, theta_, phi_ = symbols("\\rho \\theta \\phi")
```
(I use the `_` suffix as a convention for terminal symbols)
Let's construct the equations for transitioning from spherical coordinates to Cartesian ones:
```python
x = rho_ * sin(theta_) * sin(phi_)
y = rho_ * sin(theta_) * cos(phi_)
z = rho_ * cos(theta_)
print(latex(x_) + " = " + latex(x))
print(latex(y_) + " = " + latex(y))
print(latex(z_) + " = " + latex(z))
```
$
\begin{array}l
x = \rho \sin{\left(\phi \right)} \sin{\left(\theta \right)} \\
y = \rho \sin{\left(\theta \right)} \cos{\left(\phi \right)} \\
z = \rho \cos{\left(\theta \right)} \\
\end{array}
$
Let's create a new 4-dimensional manifold `spacetime`, add a `patch` to its atlas, and attach the two coordinate systems to that patch:
```python
spacetime = Manifold("spacetime", 4)
patch = Patch("patch", spacetime)
relation_dict = {
("spherical", "cart"): [(x_, y_, z_, t_), (x, y, z, t_)]
}
coord_sh = CoordSystem("spherical", patch, (t_, rho_, theta_, phi_), relation_dict)
coord_cart = CoordSystem("cart", patch, (x_, y_, z_, t_), relation_dict)
```
One inconvenience of `sympy.diffgeom` is that one can't simply use the coordinate symbols fed into the `CoordSystem` function, and instead has to take the new ones from there:
```python
c_t, c_rho, c_theta, c_phi = coord_sh.base_scalars()
print(",".join(map(latex, (c_t, c_rho, c_theta, c_phi))))
```
$
\mathbf{t},\mathbf{\rho},\mathbf{\theta},\mathbf{\phi}
$
As you can see, the symbols are the same, but only the new set of variables can be used in `sympy.diffgeom`.
Now let's define the Schwarzschild metric:
```python
Rs = symbols("r_s")
d_t, d_rho, d_theta, d_phi = coord_sh.base_oneforms()
TP = TensorProduct
metric = (
+ (1 - Rs/c_rho) * TP(d_t, d_t)
- 1 / (1 - Rs / c_rho) * TP(d_rho, d_rho)
- (c_rho**2) * TP(d_theta, d_theta)
- (c_rho**2) * (sin(c_theta)**2) * TP(d_phi, d_phi))
print(latex(metric))
```
$
\left(- \frac{r_{s}}{\mathbf{\rho}} + 1\right) \operatorname{d}t \otimes \operatorname{d}t - \sin^{2}{\left(\mathbf{\theta} \right)} \mathbf{\rho}^{2} \operatorname{d}\phi \otimes \operatorname{d}\phi - \mathbf{\rho}^{2} \operatorname{d}\theta \otimes \operatorname{d}\theta - \frac{\operatorname{d}\rho \otimes \operatorname{d}\rho}{- \frac{r_{s}}{\mathbf{\rho}} + 1}
$
Being forced to use tensor products of one-forms here is another inconvenience. One would assume having a 2d-array of $g_{ij}$ components would be easier, but alas. Differential forms it is.
Either way, once we've got the metric in this form, we can use it to work out the components of the [Levi-Civita connection](https://en.wikipedia.org/wiki/Levi-Civita_connection) (a.k.a the Christoffel symbols of the 2nd kind). A connection on a manifold defines how vectors in the tangent space at one point transform when continuously moved to a tangent space at a different point. And the Levi-Civita connection is a special kind of connection that preserves the metric tensor under such movements. Meaning that the [covariant derivative](https://en.wikipedia.org/wiki/Covariant_derivative) of the metric tensor is zero (all the while its components may change).
`sympy.diffgeom` provides us with a function `metric_to_Christoffel_2nd` that derives the Christoffel symbols for us, based on that metrics tensor preservation condition:
```python
christoffel = simplify(metric_to_Christoffel_2nd(metric))
print("\\begin{cases}")
# v v v intentional reordering here, to sort by lower indices first
for (j,k,i) in np.ndindex(*np.shape(christoffel)):
if not christoffel[i,j,k].is_zero:
print(f"\\Gamma^{{{i}}}_{{{j}{k}}} =" + latex(christoffel[i,j,k]) + " \\\\")
print("\\end{cases}")
```
$
\begin{cases}
\Gamma^{1}_{00} =\frac{r_{s} \left(- r_{s} + \mathbf{\rho}\right)}{2 \mathbf{\rho}^{3}} \\
\Gamma^{0}_{01} =\frac{r_{s}}{2 \left(- r_{s} + \mathbf{\rho}\right) \mathbf{\rho}} \\
\Gamma^{0}_{10} =\frac{r_{s}}{2 \left(- r_{s} + \mathbf{\rho}\right) \mathbf{\rho}} \\
\Gamma^{1}_{11} =\frac{r_{s}}{2 \left(r_{s} - \mathbf{\rho}\right) \mathbf{\rho}} \\
\Gamma^{2}_{12} =\frac{1}{\mathbf{\rho}} \\
\Gamma^{3}_{13} =\frac{1}{\mathbf{\rho}} \\
\Gamma^{2}_{21} =\frac{1}{\mathbf{\rho}} \\
\Gamma^{1}_{22} =r_{s} - \mathbf{\rho} \\
\Gamma^{3}_{23} =\frac{1}{\tan{\left(\mathbf{\theta} \right)}} \\
\Gamma^{3}_{31} =\frac{1}{\mathbf{\rho}} \\
\Gamma^{3}_{32} =\frac{1}{\tan{\left(\mathbf{\theta} \right)}} \\
\Gamma^{1}_{33} =\left(r_{s} - \mathbf{\rho}\right) \sin^{2}{\left(\mathbf{\theta} \right)} \\
\Gamma^{2}_{33} =- \frac{\sin{\left(2 \mathbf{\theta} \right)}}{2} \\
\end{cases}
$
With the connection components at hand, one can derive the geodesics equation. And we'll have to do it manually. Fortunately, it's not that hard.
Let some curve be defined by a parametric equation $x^\mu(t)$. Then, its tangent vector has a form $\dot{x}^\mu(t) = \frac{d}{dt}x^\mu(t)$. By the definition of geodesics, the covariant derivative of this vector along itself should be zero:
$
\nabla_{\dot{x}} \dot{x} = 0
$
Writing it out with the Christoffel symbols:
$
\nabla_{\dot{x}} \dot{x}
=
\left( \frac{\partial}{\partial{x^\mu}} \dot{x}^\nu + \Gamma^{\nu}_{\kappa\mu}\dot{x}^\kappa \right) \dot{x}^\mu
=
\left( \frac{\partial}{\partial{x^\mu}} \frac{dx^\nu}{dt} \right) \frac{dx^\mu}{dt}
+
\Gamma^{\nu}_{\kappa\mu} \frac{dx^\kappa}{dt} \frac{dx^\mu}{dt}
$
In the first term, let's swap the order of the factors, and use the chain rule:
$
\left( \frac{\partial}{\partial{x^\mu}} \frac{dx^\nu}{dt} \right) \frac{dx^\mu}{dt}
=
\frac{dx^\mu}{dt}\frac{\partial}{\partial{x^\mu}} \left( \frac{dx^\nu}{dt} \right)
=
\frac{d}{dt}\left( \frac{dx^\nu}{dt} \right)
=
\frac{d^2 x^\nu}{dt^2}
$
Then, the equation comes down to
$
\frac{d^2 x^\nu}{dt^2}
+ \Gamma^{\nu}_{\kappa\mu} \frac{dx^\kappa}{dt} \frac{dx^\mu}{dt}
= 0
$
Or, simply,
$
\boxed{
\ddot{x}^\nu =
-\Gamma^{\nu}_{\kappa\mu} \dot{x}^\kappa \dot{x}^\mu
}
$
In python:
```python
def geo_second_derivatives(christoffel, derivatives):
'''
Parameters:
christoffel: tensor of rank 3 - analytical expression of christoffel symbol in given coordinates
derivatives: iterable consisting of sympy symbols corresponding to each first derivative of a coordinate variable
'''
d2u = np.zeros(len(derivatives), dtype='object')
for n,k,m in np.ndindex(*np.shape(christoffel)):
d2u[n] -= christoffel[n,k,m] * derivatives[k] * derivatives[m]
return d2u
```
```python
d2_t, d2_rho, d2_theta, d2_phi = geo_second_derivatives(christoffel, symbols("t' \\rho' \\theta' \\phi'"))
print("\\begin{array}{l}")
print(f"{{{latex(c_t)}}}'' = {latex(d2_t)} \\\\")
print(f"{{{latex(c_rho)}}}'' = {latex(d2_rho)} \\\\")
print(f"{{{latex(c_theta)}}}'' = {latex(d2_theta)} \\\\")
print(f"{{{latex(c_phi)}}}'' = {latex(d2_phi)} \\\\")
print("\\end{array}")
```
$
\begin{array}{l}
{\mathbf{t}}'' = - \frac{\rho' r_{s} t'}{\left(- r_{s} + \mathbf{\rho}\right) \mathbf{\rho}} \\
{\mathbf{\rho}}'' = - \phi'^{2} \left(r_{s} - \mathbf{\rho}\right) \sin^{2}{\left(\mathbf{\theta} \right)} - \frac{\rho'^{2} r_{s}}{2 \left(r_{s} - \mathbf{\rho}\right) \mathbf{\rho}} - \theta'^{2} \left(r_{s} - \mathbf{\rho}\right) - \frac{r_{s} t'^{2} \left(- r_{s} + \mathbf{\rho}\right)}{2 \mathbf{\rho}^{3}} \\
{\mathbf{\theta}}'' = \frac{\phi'^{2} \sin{\left(2 \mathbf{\theta} \right)}}{2} - \frac{2 \rho' \theta'}{\mathbf{\rho}} \\
{\mathbf{\phi}}'' = - \frac{2 \phi' \rho'}{\mathbf{\rho}} - \frac{2 \phi' \theta'}{\tan{\left(\mathbf{\theta} \right)}} \\
\end{array}
$
Thus, we have arrived at the system of differential equations for the coordinates of geodesic lines.
# Raytracing
![[sphere_pole.png|Duplication of the sphere pole image]]
So, the light rays (or any other objects for that matter) follow "straight" lines defined by second-order differential equations. One can easily compute these paths numerically with [scipy.integrate](https://docs.scipy.org/doc/scipy/reference/integrate.html) in python.
But wouldn't we want to **see** the black hole? One could send a ray through every pixel of a screen, and let them go around the black hole and hit its surroundings. And if the rays are sent back in time, it's as if one received the light emitted by the surrounding objects in the past!
Well, we would have to send a million of rays then. Not an easy task for a python script. If only we could run all of it in parallel! Well, that's what GPUs are for.
For running the same calculations on every pixel, there is a special kind of GPU program: a fragment shader. To make our own shader, let's take OpenGL and its shading language GLSL.
The plan is extremely simple: place a colorful spherical shell around the black hole, put a camera next to it, send light rays from the camera, and let them either hit the shell, giving colors to the corresponding pixels, or be lost somewhere on the event horizon.
## Geodesics equations in GLSL...
Nobody wants to retype the equations that we got from sympy into a different language. Fortunately, we can just use the [ccode](https://docs.sympy.org/latest/modules/codegen.html) function to produce the code, and then it only requires some minor cleanup. After grooming the code, we get the following function for our geodesics equation $p'' = f(p, p')$:
```c
vec4 d2p(vec4 p, vec4 dp)
{
float t = p[0];
float rho = p[1];
float theta = p[2];
float phi = p[3];
float D_t = dp[0];
float D_rho = dp[1];
float D_theta = dp[2];
float D_phi = dp[3];
// Geodesics equation:
float D2_t = - (D_rho * rs * D_t) / (rho * (rho - rs));
float D2_rho = (rho - rs) * (
+ pow(D_phi * sin(theta), 2)
+ pow(D_theta, 2)
- (rs * pow(D_t, 2)) / (2 * pow(rho, 3))
)
+ (rs * pow(D_rho, 2)) / (2 * rho * (rho - rs));
float D2_theta = pow(D_phi, 2) * sin(theta) * cos(theta) - (2 * D_rho * D_theta) / rho;
float D2_phi = -(2 * D_phi * D_theta) / tan(theta) - (2 * D_rho * D_phi) / rho;
vec4 d2p = {D2_t, D2_rho, D2_theta, D2_phi};
return d2p;
}
```
This function takes the coordinates of a point on the light ray, and the velocity 4-vector at that point. It returns the second derivatives of the coordinates as they move along the curve.
Here one can immediately spot an issue: see the division by $tan(\theta)$? When this angle is zero, the second derivative value would be undefined. And even when $\theta$ is close to 0, it would cost us precision. Fortunately, it's not a big deal. We'll navigate around the precision issue by using an adaptive integration step. Alternatively, one could notice that for non-rotating black holes, the rays never leave their plane, and use that to avoid the zeros problem, but I ~~can't be bothered~~ will leave it as an exercise for the reader.
Now, we can compute the second derivatives for any geodesics passing through a point `p` in the direction given by a 4-velocity `dp`. This second derivative determines how that 4-velocity changes, and the 4-velocity defines what the next point is going to be. The classic numerical integration problem.
## ... and its numerical solution
In the quest to find balance between accuracy and performance, I tried out several different integration methods and eventually settled on the [adaptive Runge-Kutta method](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#Adaptive_Runge%E2%80%93Kutta_methods) of orders 3 and 4.
The idea behind adaptive Runge-Kutta methods is that two methods of different orders (in this case, 3rd and 4th) run in parallel. The difference between their results gives an estimate of the error, which is then used to adjust the integration step: if the error is too big, the step gets smaller; if it’s small enough, the step gets bigger.
What makes it doubly convenient is that intermediate results from the method of one order can be reused - with different coefficients - inside the method of the other order, to boost performance. Choosing those coefficients is a science in itself, and there are hundreds of research papers written on this topic.
Naturally, I didn’t read those hundreds of papers. I just took the classic RK4 as the main method, and then hacked together a third-order method to reuse as much of what's already computed as possible.
>[!info]- On the RK4 method
> RK4 is a method for solving the initial value problem: finding a curve $y(t)$ that satisfies the differential equation $y'(t) = f(t, y(t))$ and goes through the point $(t_0, y_{0})$. The easiest method to numerically solve it is the Euler method:
>$y_{i+1} = y_{i} + f(t_{i}, y_{i}) \Delta t$
> Unfortunately, the higher the second derivative value is, the worse this method works. And it works especially bad if the curve keeps bending in the same direction - and that's precisely our case.
> However, one might add a few corrections to it: compute the derivatives not only at the point $(t_{i}, y_{i})$, but around it too. Consider this:
> 1. $a = f(t_i,y_i)$ - the coefficient of the Euler method
> 2. $b=f\left( t_{i} + \frac{\Delta t}{2}, y_{i} + a \frac{\Delta t}{2} \right)$ - go half way to the resulting point of Euler method and compute the derivative value here.
> 3. Now use this $b$ value to move half-step from the initial point: $y_{b} = y_{i} + b \frac{\Delta t}{2}$
> 4. And compute the derivative value there:
> $c=f\left( t_{i} + \frac{\Delta t}{2}, y_{i} + b \frac{\Delta t}{2} \right)$
> 5. Finally, use $c$ to take a full step from the initial point and compute another derivative there:
$d = f(t_{i} + \Delta t, y_{i} + c \Delta t)$
> 6. With these 4 values, take a weighted average of $a,b,c,d$ to be the actual derivative value at $(t_{i}, y_{i})$:
$y_{i}'(t) = \frac{1}{6}(a + 2b + 2c + d)$
>
>And only with this value, the actual step is taken:
>$y_{i+1} = y_{i} + \frac{1}{6}(a + 2b + 2c + d) \Delta t$
> This makes the classic RK4. Naturally, one could pick different weights and steps, creating a huge family of other Runge-Kutta methods.
>
> For the 3-rd order method running along with the main RK4 I simply chose
>$y_{i+1} = y_{i} + \frac{1}{6}(a + 2b + 3c) \Delta t$
But since RK4 solves first-order equations, how is it going to solve a second-order one? Not a problem! Any equation of the form
$
p'' = f(\tau, p, p')
$
can be reduced to a first order equation by adding extra dimensions, as one can consider $p'$ to be an independent variable:
$
y =
\begin{bmatrix}
p \\ p'
\end{bmatrix}
$
Then,
$
y' =
\begin{bmatrix}
p \\ p'
\end{bmatrix}'
=
\begin{bmatrix}
p' \\ p''
\end{bmatrix}
=
\begin{bmatrix}
p' \\ f(\tau,p,p')
\end{bmatrix}
=g(\tau,y)
$
This equation is now of the first order, so we can solve it.
Let's make a wrapper for `d2p` function do be used in this way:
```c
void rk4_diff(vec4 y[2], out vec4 dy[2])
{
dy[0] = y[1];
dy[1] = d2p(y[0], y[1]);
}
```
This wrapper takes an 8-dimensional vector that contains the point position and velocity, and returns its derivative: velocity as the derivative of position, and the "acceleration" computed with `d2p`.
The resulting RK34 integrator isn't the prettiest, but here it is:
>[!example]- rk34_step
> ```c
> /**
> * Adaptive Runge-Kutta of orders 3 and 4, tailored for taking 2nd order derivative for a moving 4D point.
> * [in] epsilon - desired value for error estimate. Step is updated based on that value
> * [in] max_step - hard limit to how big a step can get
> * [in] dtau - integration step
> * [inout] p0 - point position
> * [inout] p1 - point velocity
> * [out] error_estimate - estimated error value (difference between RK3 and RK4 results)
> *
> * Return value:
> * true - error is not too big, can continue processing
> * false - error is too big, need to rerun with the updated dtau
> */
>
> bool rk34_step(
> float epsilon, float max_step, inout float dtau, inout vec4 p0, inout vec4 p1, out float error_estimate)
> {
> vec4 a[2], b[2], c[2], d[2];
>
> vec4 y0[2] = {p0, p1};
> rk4_diff(y0, a);
>
> vec4 y1[2] = {p0 + a[0] * dtau/2, p1 + a[1] * dtau/2};
> rk4_diff(y1, b);
>
> vec4 y2[2] = {p0 + b[0] * dtau/2, p1 + b[1] * dtau/2};
> rk4_diff(y2, c);
>
> vec4 y3[2] = {p0 + c[0] * dtau, p1 + c[1] * dtau};
> rk4_diff(y3, d);
>
> vec4 rk3_0 = (a[0] + 2 * b[0] + 3 * c[0]) * dtau / 6;
> vec4 rk3_1 = (a[1] + 2 * b[1] + 3 * c[1]) * dtau / 6;
>
> vec4 rk4_0 = (a[0] + 2 * b[0] + 2 * c[0] + d[0]) * dtau / 6;
> vec4 rk4_1 = (a[1] + 2 * b[1] + 2 * c[1] + d[1]) * dtau / 6;
>
> vec4 diff0 = rk3_0 - rk4_0;
> vec4 diff1 = rk3_1 - rk4_1;
> error_estimate = sqrt(abs(dot(diff0, diff0) + dot(diff1, diff1)));
>
> float new_dtau = dtau*sqrt(epsilon/error_estimate);
> if (new_dtau < dtau / 2) {
> dtau = new_dtau;
> return false;
> }
> dtau = new_dtau;
> if (dtau > max_step) {
> dtau = max_step;
> }
> p0 += rk4_0;
> p1 += rk4_1;
>
> return true;
> }
> ```
# Coordinates
![[sphere_with_back.png|Spherical shell around a black hole]]
So, we've figured out how to take the point position and velocity, and feed them into the `rk34_step` integrator. It takes a small step along the geodesics and adjusts the size of the next step to keep the error within bounds.
Now we have all the tools needed for sending rays of light wherever we want. But no tools to know where to actually send them.
And for that, like before, we need two coordinate systems: spherical for computations, and Cartesian for working with the camera. And a way to switch between them.
## Point coordinates transformation
We already know how to get from a spherical coordinate system to Cartesian one. Let's just turn it into GLSL code:
```c
vec4 sh_to_cart(vec4 sh)
{
float t = sh[0];
float r = sh[1];
float th = sh[2];
float phi = sh[3];
float x = r*sin(th)*sin(phi);
float y = r*sin(th)*cos(phi);
float z = r*cos(th);
vec4 ret = {x, y, z, t};
return ret;
}
```
For Cartesian coordinates, the time is chosen to be the last one so that one can access $x,y,z$ just like regular GLSL 4-vector components: `point.x`, `point.y`, `point.z`. $t$ can be accessed via `point.w`.
To determine where to send our light ray, we'll need an inverse transform as well: from Cartesian coordinates (of the pixels on the screen) to spherical ones.
This is something sympy struggles with. But it's a pretty easy geometry problem for a human:
```c
vec4 cart_to_sh(vec4 pos)
{
float r = sqrt(pos.x*pos.x + pos.y*pos.y + pos.z*pos.z);
float th = PI / 2. - atan(sqrt(pos.x*pos.x + pos.y*pos.y), pos.z);
float phi = atan(pos.y, pos.x);
float t = pos.w;
vec4 ret = {t, r, th, phi};
return ret;
}
```
Note how we are feeding 2 arguments to the `atan` function. In GLSL, that's a shortcut for finding the angle of a triangle given its legs (same as the `atan2` function in many other languages). Sympy can't really use this, and using a regular tangent would lose information about the signs, so the machine can't figure out the inverse transform.
## Vector coordinates transformation
To reiterate, geodesics (i.e. the light rays) are paths of vectors being transported along themselves. But in curvilinear coordinates, the same vector starting at different points would have different components, as the coordinate lines themselves change. So at each point there is a unique way to go between different coordinate systems in its tangent space. This transformation is captured by a [Jacobian matrix](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant).
And again, sympy comes to aid us:
```python
J = coord_sh.jacobian(coord_cart)
IJ = simplify(J.inv())
print(f"J_{{spherical \\rightarrow cart}} = {latex(J)}")
print(f"J_{{cart \\rightarrow spherical}} = {latex(IJ)}")
```
$
J_{spherical \rightarrow cart} = \left[\begin{matrix}0 & \sin{\left(\phi \right)} \sin{\left(\theta \right)} & \rho \sin{\left(\phi \right)} \cos{\left(\theta \right)} & \rho \sin{\left(\theta \right)} \cos{\left(\phi \right)}\\0 & \sin{\left(\theta \right)} \cos{\left(\phi \right)} & \rho \cos{\left(\phi \right)} \cos{\left(\theta \right)} & - \rho \sin{\left(\phi \right)} \sin{\left(\theta \right)}\\0 & \cos{\left(\theta \right)} & - \rho \sin{\left(\theta \right)} & 0\\1 & 0 & 0 & 0\end{matrix}\right]
$
$
J_{cart \rightarrow spherical} = \left[\begin{matrix}0 & 0 & 0 & 1\\\sin{\left(\phi \right)} \sin{\left(\theta \right)} & \sin{\left(\theta \right)} \cos{\left(\phi \right)} & \cos{\left(\theta \right)} & 0\\\frac{\sin{\left(\phi \right)} \cos{\left(\theta \right)}}{\rho} & \frac{\cos{\left(\phi \right)} \cos{\left(\theta \right)}}{\rho} & - \frac{\sin{\left(\theta \right)}}{\rho} & 0\\\frac{\cos{\left(\phi \right)}}{\rho \sin{\left(\theta \right)}} & - \frac{\sin{\left(\phi \right)}}{\rho \sin{\left(\theta \right)}} & 0 & 0\end{matrix}\right]
$
With this, for any point defined by its spherical coordinates, there is a matrix that determines how the vector coordinates are transformed.
Well, except for the points with $\theta = 0$, since there's division by zero. But let's not worry about the details.
For use in GLSL, these matrices have to be transposed, since they are stored in column-major format:
```c
// Value of Jacobi matrix at point sh_pos (point is given in Schwarzschild coordinates)
// cartesian_vector = J * sh_vector
mat4 jacobian_inner_to_cart_at(vec4 sh_pos)
{
float rho = sh_pos[1];
float theta = sh_pos[2];
float phi = sh_pos[3];
// OpenGL stores matrices in column-major format. So this is a TRANSPOSED jacobian matrix
mat4 ret = {
{0, 0, 0, 1},
{sin(phi)*sin(theta), sin(theta)*cos(phi), cos(theta), 0},
{rho*sin(phi)*cos(theta), rho*cos(phi)*cos(theta), -rho*sin(theta), 0},
{rho*sin(theta)*cos(phi), -rho*sin(phi)*sin(theta), 0, 0},
};
return ret;
}
// Value of inverse Jacobi matrix at point sh_pos (point is given in Schwarzschild coordinates)
// sh_vector = IJ * cartesian_vector
mat4 jacobian_cart_to_inner_at(vec4 sh_pos)
{
float rho = sh_pos[1];
float theta = sh_pos[2];
float phi = sh_pos[3];
// OpenGL stores matrices in column-major format. So this is a TRANSPOSED jacobian matrix
mat4 ret = {
{0, sin(theta)*sin(phi), sin(phi)*cos(theta)/rho, cos(phi)/(rho*sin(theta))},
{0, sin(theta)*cos(phi), cos(phi)*cos(theta)/rho, -sin(phi)/(rho*sin(theta))},
{0, cos(theta), -sin(theta)/rho, 0},
{1, 0, 0, 0},
};
return ret;
}
```
# Drawing the rest of the owl
![[outer_shell.png|Spherical shell just above the event horizon. Its whole surface is visible when looking from one side.]]
The fragment shader will run for every pixel of a framebuffer and will output the color for that pixel.
```c
layout(location = 0) out vec4 diffuseColor;
#define PI 3.1415926535897932384626433832795
uniform vec2 window_size;
uniform vec4 camera_pos;
uniform mat3 camera_rotation;
uniform float viewport_dist = 1.0;
uniform float rs = 0.25;
uniform float shell_radius = 2.0;
```
The inputs are provided as the following [uniform](https://www.khronos.org/opengl/wiki/Uniform_(GLSL)) variables:
- size of the framebuffer
- camera parameters: position defined in spherical coordinates, the camera rotation matrix in Cartesian ones, and the viewport distance
- Schwarzschild radius
- radius of the shell
Coordinates of the pixel itself are available in the `gl_FragCoord` variable.
The rest is pretty straightforward:
- Send a ray through each pixel, back in time
- In a loop, move it along a geodesic line using `rk34_step` function
- Upon hitting the shell, paint the pixel with its color. If the ray is leaving the shell, return.
Let's go!
>[!example]- void main()
> ```c
> void main()
> {
> // Position of the current pixel in the (-1, 1) range
> vec2 window_center = {window_size.x / 2, window_size.y / 2};
> vec2 pixel_pos = (vec2(gl_FragCoord) - window_center) / max(window_center.x, window_center.y);
>
> // Light ray.
> // Direction from the viewport into the screen pixel:
> vec3 direction3 = {viewport_dist, pixel_pos.x, pixel_pos.y};
> direction3 = normalize(direction3);
> // Adjusted for the camera pointing direction:
> direction3 = camera_rotation * direction3;
> vec4 direction4 = {direction3.x, direction3.y, direction3.z, -length(direction3)};
>
> // Current position and 4-velocity in Schwarzschild coordinates:
> vec4 pos_s = camera_pos;
> vec4 dp = jacobian_cart_to_inner_at(pos_s) * direction4;
>
> // Initial step
> float dtau = 0.01;
>
> vec3 color = {0,0,0};
> float prev_pho = pos_s[1];
> float max_step = 0.25;
> int i;
> for (i = 0; i < 2000; i++) {
> float r = pos_s[1];
>
> if (abs(r - shell_radius) < 1) {
> // increase precision near the shell
> max_step = 0.01 + abs(r - shell_radius)/20.;
> } else {
> max_step = 0.5;
> }
>
> float error_estimate;
> if (rk34_step(0.01, max_step, dtau, pos_s, dp, error_estimate) == false) {
> // Spend this iteration on adjusting the step size.
> continue;
> }
> if ((r > shell_radius * 2) && (r - prev_pho > 0)) {
> // Optimization for looking from outside.
> // This ray has missed the shell completely and keeps moving away from it.
> break;
> }
>
> // Optimization. If we are this close, not going away anyway
> if (r < rs*1.01) {
> break;
> }
>
> bool entering_shell = (prev_pho >= shell_radius) && (r < shell_radius);
> bool exiting_shell = (prev_pho < shell_radius) && (r >= shell_radius);
> if ( entering_shell || exiting_shell) {
> float lat = pos_s[2]/PI*180; // in degrees
> float lon = pos_s[3]/PI*180; // in degrees
> float line_freq = 15; // degrees between grid lines
>
> // Grid lines
> float shine = 1.0;
> float thicc = 0.5;
>
> float prev_lat_line = floor(lat/line_freq)*line_freq;
> float next_lat_line = prev_lat_line + line_freq;
> color.b += shine * exp(-pow((lat - prev_lat_line)/thicc, 2));
> color.b += shine * exp(-pow((lat - next_lat_line)/thicc, 2));
>
> float prev_lon_line = floor(lon/line_freq)*line_freq;
> float next_lon_line = prev_lon_line + line_freq;
> color.b += shine * exp(-pow((lon - prev_lon_line)/thicc, 2));
> color.b += shine * exp(-pow((lon - next_lon_line)/thicc, 2));
>
> color.g = color.b/2;
>
> if (exiting_shell)
> break;
> }
> prev_pho = r;
> }
>
> diffuseColor = vec4(color, 1.0);
> }
> ```
The results:
![[sphere_1.png|]]
This shader is actually sub-optimal. Division by near-zero values ruins performance near the poles, as it forces the integrator to choose very small steps to preserve precision. Just look at this evil red glow! It's proportional to the number of iterations spent on rendering each pixel.
![[sphere_with_iteration_count.png|]]
And here's the error estimate from the integrator:
![[sphere_with_max_error.png|]]
The rest of the pictures in this article (except for the torus) were made by just tweaking the shell size and colors.
I'd say there’s room to make this quite a few times faster. But I’ll leave that as an exercise for the reader.
<div class=custom-page-footer><hr/>
<div class=spaced-lr>
<a href="https://some.3b1b.co/">#SoME4</a>
June 2025
</div>
</div>