SIM5 - a general purpose library for GR raytracing and radiation transport

Introduction

SIM5 is a C library with a Python interface that contains a collection of routines for relativistic raytracing and radiation transfer in GR. It has a special focus on raytracing from accretion disks, tori, hot spots or any other 3D configuration of matter in Kerr geometry, but it can be used with any other metric as well. It can handle both optically thick and thin sources as well as transport of polarization of the radiation and helps to calculate the propagation of light rays from the source to an observer through a curved spacetime. It supports parallelization and runs on GPUs.

Features

The library contains

• a selection of physical and unit conversion constants related to black hole disk accretion
• an implementation of relativistic thin disk model
• basic routines related to GR and Kerr geometry in particular
• routines for computing null geodesics and GR raytracing with polarization

The library is thread-safe and supports parallelization of the calculations at both CPUs (OpenMP/MPI) and GPUs (CUDA/OpenCL).

Installation

Acquire the source code by cloning the git repository

git clone https://github.com/mbursa/sim5

The C code does not have any external library dependencies except the standard system ones. Typically, it is enough to run

make

The compilation process creates a lib folder and puts two files there, sim5lib.c, sim5lib.h and sim5lib.o, which contain a complete C source code merged by the compiler process from all the individual source files in src folder. Use those files to statically link SIM5 in your code (most typical usage in C/C++ codes), see an example bellow.

Compilation of the Python interface

The Python interface does have a few external dependecies and thus it does not compile by default when calling make. Before compiling the Python interface, make sure you have the following pre-requisities installed:

    # python headers
    apt install python3-dev
    # SWIG - a tool that connects C/C++ programs with high-level programming languages
    apt install swig

Having those do

    make python

<!-- The third file is a compiled library used for run-time linking. It is used by the Python interface, but you may use it for dynamical linking of SIM5 to your code too. -->

Usage

The library can be directly used in C/C++ projects and with a proper interface it can be used also in other languages (Fortran, ...). It also works as a Python package making most of its functions accessible from Python scripts with no extra effort.

Note: only Python 3 is supported

Using SIM5 in your C/C++ project (static linking)

The most typical use of the code in a C/C++ project is to use static linking. Having a sample code that uses SIM5

#include <stdlib.h>
#include "sim5lib.h"

int main(int argc, char *argv[])
//***************************************************
{
    double a, rms, rbh;

    // for BH spin values from 0 to 1 print the radius of BH horizon and 
    // the location of the marginally stable orbit
    for (a=0.0; a<1.0; a+=0.01) {
        rbh = r_bh(a);  // radius of the event horizon
        rms = r_ms(a);  // radius of the marginally stable orbit
        fprintf(stderr, "%.4e  %e  %e\n", a, rbh, rms);
    }
    
    return 0;
}

we compile the code with (assuming SIM5 is located in subfolder sim5 of the current directory and assuming make has been run on SIM5)

gcc kerr-orbits.c -o kerr-orbits -I./sim5/lib -L./sim5/lib ./sim5/lib/sim5lib.o -lm

Refer to examples for further inspiration.

Using SIM5 in your Python project

SIM5 can be used as a Python module, which allows its routines to be called from a Python script. It is enough to import sim5 in your script and call the member functions. Some care has to be taken when working with arrays and objects that are passed to SIM5 routines or returned from them. The SIM5 Python interface uses SWIG tool to make an interface between C and Python and you may refer to SWIG documentation to learn about how to do that.

A basic example of a Python script that calls SIM5 may look like this

import sys
import math
import sim5

bh_spin = 0.99    # BH spin parameter
bh_mass = 10.0    # BH mass in units of Solar mass
mdot    = 0.1     # disk accreation rate in units of Eddington accreation rate
alpha   = 0.1     # alpha parameter for viscosity
rout    = 500     # outer disk boudary in units of gravitational radius

# setup a relativistic thin disk model
sim5.disk_nt_setup(bh_mass, bh_spin, mdot, alpha, 0)

print '# Flux from the accretion disk surface as a function of radius for Novikov-Thorne disk model'
print '# mdot =', sim5.disk_nt_mdot(), '[Mdot_Edd]; ', sim5.disk_nt_mdot()*bh_mass*sim5.Mdot_Edd/1e18, 'g/s'
print '# format: radius [GM/c2] local flux [erg/s/cm2]'

r = sim5.disk_nt_r_min()
while (r < rout):
    print r, sim5.disk_nt_flux(r)
    r = r * 1.05
#end of while 

This example assumes that sim5 subfolder does exist in the same folder along with the script. If that is not the case, put a line specifying the path to SIM5 library before the import statement.

sys.path.append('/path/to/sim5')
import sim5
...

Python calls to SIM5 routines are handled by the compiled C code in lib/sim5lib.so library. Although there is some overhead connected with the call to SIM5 routines, the actual evaluation of SIM5 functions is as fast as C code can be. This makes complex tasks (like raytracing) reasonably fast even if they are coded in Python (see examples in demo folder).

Note: only Python 3 is supported

Raytracing with SIM5

SIM5 library implements two essential ways of raytracing photons through curved spacetimes:

• by evaluating positional elliptic integrals along a defined trajectory
• by a stepwise integration of the geodesic equation for null geodesics

Both approaches have their advatages and disadvantages and the choice of one or the other depends mainly on the specific task that has to be carried out. In general, the use of the stepwise integration of the geodesic equation is necessary when the spacetime is not described by Kerr metric, otherwise the required accuracy, performance and computational complexity have to be considered to make a choice.

Raytracing by evaluating positional elliptic integrals

In Kerr spacetime, photon motion in r-θ plane is governed by equation
where R(r) and Θ(θ) are polynomials, and a solution can be expressed in terms of Jacobi elliptic functions. Any photon geodesic is determined by it two constants of motion (traditionally denoted λ and q) and a position along the geodesic can be uniquely expressed using the value of the above positional integral. In this way, once the trajectory is fixed by its motion constants, it is possible to follow the photon path by changing the value of the positional integral and by inverting it to obtain spacetime coordinates:

// initialize the geodesic with an initial position `x` and direction `k`
geodesic_init_src(bh_spin, x[1], x[2], k, k[1]<0?1:0, &gd, &error);
// get the initial value of the positional integral
P = geodesic_P_int(&gd, x, 0);

r = x[1];   // radial coordinate
m = x[2];   // poloidal cordinate, m=cos(theta)
while (r < R_MAX) {
    // make a step along the geodesic; both `P`, `r` and `m` are updated
    geodesic_follow(&gd, deltaP, &P, &r, &m, &status);   
    if (!status) break;
    ...	 // do some stuff at the new place
}

Refer to the following examples:

• [tbd]

Equatorial plane

An important special case of the above approach is raytracing of photons that originate from the equatorial plane, as this is a commonly adopted approximation e.g. in case of dealing with accretion disks. Integrating over the observer's sky, we can directly determine the radius from which the photons originate by an analytic inversion of an elliptic positional integral. Thus, by using

geodesic_init_inf(inclination, spin, alpha, beta, &gd, &error);
P = geodesic_find_midplane_crossing(&gd, 0);
r = geodesic_position_rad(&gd, P);

one can quickly and efficiently find the original radius of a photon coming out of the disk surface that has impact parameters alpha and beta on the observer's image plane.

Refer to the following examples:

• Example 04 - accretion disk image
• Example 05 - accretion disk spectrum

Raytracing by stepwise geodesic equation integration

Having an initial position x^μ and direction of propagation dx^μ/dλ of the photon, its subsequent trajectory through the spacetime can be computed by solving a second-order ordinary differential equation (so called geodesic equation): SIM5 uses Verlet integration method to numericaly integrate geodesic equations because it offers good performance while keeping reasonable precission (see Dolence at al. 2009). The basic scheme of the integration is the following (in Python):

# initialize trajectory from initial position `x` and direction `k`
sim5.raytrace_prep