As usual in raytracing, each sample (pixel) in the image corresponds to a path starting from the position ${p}_{0}=\left({x}_{0},{y}_{0},{z}_{0}\right)$ of the camera, initially travelling to a direction ${d}_{0}\in {\mathbb{R}}^{3}$ specified by the camera coordinate system. The path represents photons arriving to the camera from that direction. Instead of traveling in straight lines until hitting an object, the paths solve the photonic geodesic equation in the Schwarzschild space-time, which in certain units and spherical coordinates $\left(t,r,\mathit{\theta},\varphi \right)$ takes the form (cf. [1, §4.2.1])

$$d{s}^{2}=\left(1-\frac{1}{r}\right)d{t}^{2}+{\left(1-\frac{1}{r}\right)}^{-1}d{r}^{2}-{r}^{2}\left(d{\mathit{\theta}}^{2}-{sin}^{2}\mathit{\theta}\cdot d{\varphi}^{2}\right)=0.$$

Without loss of generality, one can fix $\mathit{\theta}=\frac{\pi}{2}$ and, even though not very obviously (see, e.g., [3] and [2]), this gives the ordinary differential equation

$${u}^{\u2033}\left(\varphi \right)=-u\left(\varphi \right)\left(1-\frac{3}{2}{u}^{2}\left(\varphi \right)\right),$$ | (1) |

for $u:=1\u2215r$. To solve it in the raytracer, I use the initial conditions

$$u\left(0\right)=\frac{1}{\parallel {p}_{0}\parallel},\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}{u}^{\prime}\left(0\right)=-u\left(0\right)\frac{{d}_{0}\cdot n}{{d}_{0}\cdot t},$$where

span the plane in which the ray travels.

Equation (1) is integrated using the Leapfrog method until either a maximum number of steps $N$ is reached or ${u}_{j}+\Delta {u}_{j}<0$, which signals that the ray escapes at an angle $\varphi \in \left({\varphi}_{j},{\varphi}_{j}+\Delta {\varphi}_{j}\right)$. Given $\varphi $ and $u$, the corresponding spatial coordinates $p=\left(x,y,z\right)$ can be computed from (2) as

$$p=\left(ncos\varphi +tsin\varphi \right)\u2215u.$$The Schwarzschild time coordinate $t$ is resolved by amending the numerical integration with

$${t}^{\prime}\left(\varphi \right)=-\frac{\sqrt{{\left({u}^{\prime}\right)}^{2}+{u}^{2}\left(1-u\right)}}{{u}^{2}\left(1-u\right)},$$ which is needed to compute light travel times for simulating apparent positions of
moving objects. In order to do the integration in a visually accurate way, the step
size $\Delta {\varphi}_{j}$
is varied using an *ad hoc* formula, in addition to using the classical approximation
$d{t}^{2}\approx d{x}^{2}+d{y}^{2}+d{z}^{2}$ at large
distances $r$
from the black hole.

When the observer is moving with velocity
$v$
w.r.t. the Schwarzschild coordinate system, the ray directions
$d$ in
the coordinate system of the observer (the camera coordinate system) are
transformed by *relativistic aberration*, which is corresponds to the Lorentz
velocity transformation

where $v={v}^{0}v$, ${d}_{\parallel}={v}^{0}\cdot d$, ${d}_{\perp}=d-{d}_{\parallel}$ and

$$\gamma =\frac{1}{\sqrt{1-{v}^{2}}}.$$ | (3) |

The observer (unless stationary) and the planet both move on circular orbits with angular velocities given by [2]

$$\frac{d\varphi}{dt}=\frac{1}{r\sqrt{2\left(r-1\right)}}.$$ | (4) |

On such orbits, the time dilation of the observer relative to the planet is computed by updating the $t$-coordinate as

$$\Delta t=\sqrt{\frac{{\left(\Delta {t}_{0}\right)}^{2}\left(1-{v}^{2}\right)}{1-\frac{1}{r}}},$$where $\Delta {t}_{0}$ is the elapsed wall clock time between consecutive animation frames and $v=rd\varphi \u2215dt$.

When computing intersection between the planet and light paths, the motion of the planet has to be taken into account. Within the ODE integration steps, a ray-sphere intersection model can be used after mapping the ray to the moving coordinate system of the planet using a Lorentz velocity transformation.

The *Doppler factor* for light arriving to the observer from direction
$d$, emitted by a source
moving at velocity $v$
is [1, §1.7]

$$\delta =\gamma \left(1+d\cdot v\right),$$ | (5) |

where $\gamma $ is the Lorentz factor (3). The Doppler factor for a ray path connecting a moving light source and a moving observer can be computed as

$$\delta ={\delta}_{o}{\delta}_{s}={\gamma}_{o}{\gamma}_{s}\left(1-{d}_{o}\cdot {v}_{o}\right)\left(1+{d}_{s}\cdot {v}_{s}\right),$$where ${d}_{o}$ and ${d}_{s}$ are the initial and final (at intersection with the light source) photon directions on the ray path. The velocities ${v}_{o}$ and ${v}_{s}$ of the observer and the source can be determined by (4). All of these are given in the ‘‘stationary’’ Schwarzschild coordinate system. For the accretion disk, the velocity ${v}_{s}$ is computed assuming the disk consists of glowing matter moving in circular orbits given by (4).

The wavelengths of the emitted and observed light are related by
${\lambda}_{o}=\delta {\lambda}_{s}$ and
intensities by $I\left({\lambda}_{o}\right)={\delta}^{-3}I\left({\lambda}_{s}\right)$
[1]. The change in intensity is called *relativistic beaming* and the change in
wavelength is Doppler shift, which is also commonly referred to as *redshift* (when
$\delta >1$) or
*blueshift* ($\delta <1$).

The perceived (RGB) colors corresponding to any fixed spectrum $I$ can be precomputed in a one-dimensional lookup table (texture) as a function of the Doppler factor. This is done by integrating the shifted spectra with respect to the standard CIE color sensitivity measures [4]. The spectra used for the colors in this simulation are the black-body spectrum

$$I\left(\lambda \right)\propto \frac{2h{c}^{2}}{{\lambda}^{5}}\frac{1}{{e}^{\frac{hc}{\lambda {k}_{B}T}}-1}$$and hydrogen-$\alpha $: $I\left(\lambda \right)\propto {\delta}_{656.281nm}\left(\lambda \right)$.

[1] Jeremy Goodman, Topics in High-Energy Astrophysics http://www.astro.princeton.edu/~jeremy/heap.pdf

[2] Riccardo Antonelli, Visualizing a Black Hole http://spiro.fisica.unipd.it/~antonell/schwarzschild/

[3] Wikipedia, Schwarzschild Geodesics https://en.wikipedia.org/wiki/Schwarzschild_geodesics

[4] Munsell Color Science Laboratory, CIE data http://www.cis.rit.edu/research/mcsl2/online/cie.php

All web pages retrieved in December 2015.