As usual in raytracing, each sample (pixel) in the image corresponds to a path starting from the position of the camera, initially travelling to a direction 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 takes the form (cf. [1, §4.2.1])
Without loss of generality, one can fix and, even though not very obviously (see, e.g.,  and ), this gives the ordinary differential equation
for . To solve it in the raytracer, I use the initial conditions
span the plane in which the ray travels.
Equation (1) is integrated using the Leapfrog method until either a maximum number of steps is reached or , which signals that the ray escapes at an angle . Given and , the corresponding spatial coordinates can be computed from (2) as
The Schwarzschild time coordinate is resolved by amending the numerical integration with
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 is varied using an ad hoc formula, in addition to using the classical approximation at large distances from the black hole.
When the observer is moving with velocity w.r.t. the Schwarzschild coordinate system, the ray directions 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 , , and
The observer (unless stationary) and the planet both move on circular orbits with angular velocities given by 
On such orbits, the time dilation of the observer relative to the planet is computed by updating the -coordinate as
where is the elapsed wall clock time between consecutive animation frames and .
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 , emitted by a source moving at velocity is [1, §1.7]
where 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
where and are the initial and final (at intersection with the light source) photon directions on the ray path. The velocities and 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 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 and intensities by . 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 ) or blueshift ().
The perceived (RGB) colors corresponding to any fixed spectrum 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 . The spectra used for the colors in this simulation are the black-body spectrum
and hydrogen-: .
 Jeremy Goodman, Topics in High-Energy Astrophysics http://www.astro.princeton.edu/~jeremy/heap.pdf
 Riccardo Antonelli, Visualizing a Black Hole http://spiro.fisica.unipd.it/~antonell/schwarzschild/
 Wikipedia, Schwarzschild Geodesics https://en.wikipedia.org/wiki/Schwarzschild_geodesics
 Munsell Color Science Laboratory, CIE data http://www.cis.rit.edu/research/mcsl2/online/cie.php
All web pages retrieved in December 2015.