GitHub

Geodesic integration strategies

Second-Order

The motivation behind the second-order methods is to permit the computation of geodesics in generic spacetimes, via the geodesic equation:

Gradus.compute_geodesic_equation — Function
compute_geodesic_equation(ginv, j1, j2, v)

Using the inverse metric ginv, the Jacobian of the metric for $r$ and $\theta$, j1 and j2 respectively, and velocity four-vector v, calculates the four-acceleration via the geodesic equation.

Returns the components of $\frac{\text{d}^2 u^\mu}{\text{d} \lambda^2}$ via

\[\frac{\text{d}^2 u^\mu}{\text{d} \lambda^2} + \Gamma^{\mu}_{\phantom{\mu}\nu\sigma} \frac{\text{d}u^\nu}{\text{d} \lambda} \frac{\text{d}u^\sigma}{\text{d} \lambda} = 0,\]

where $x^\mu$ is a position four-vector, $\Gamma^{\mu}_{\phantom{\mu}\nu\sigma}$ are the Christoffel symbols of the second kind, and $\lambda$ the affine parameter describing the curve.

The Christoffel symbols $\Gamma^{\mu}_{\phantom{\mu}\nu\sigma}$ are defined

\[\Gamma^{\mu}_{\phantom{\mu}\nu\sigma} = \frac{1}{2} g^{\mu\rho} \left( \partial_{\nu}g_{\rho \sigma} + \partial_{\sigma}g_{\rho \nu} - \partial_{\rho}g_{\sigma \nu} \right).\]

Limitations:

  • currenly pre-supposes static, axis-symmetric metric.
source

The above can be solved as a second-order ODE, subject to an initial position and initial velocity

\[u^\mu = \left(t, r, \theta, \phi \right), \quad \text{and} \quad \dot{u}^\mu = \left( \dot{t}, \dot{r}, \dot{\theta}, \dot{\phi} \right),\]

where the dot refers to the derivative with respect to $\lambda$. In general, the spatial components of the initial velocity are known a priori, and the time-component is determined via the constraint:

Gradus.constrain_time — Function
constrain_time(g_comp, v, μ = 0.0, positive::Bool = true)

Constrains the time component of the four-velocity v, given metric components g_comp and effective mass μ.

\[g_{\sigma\nu} \dot{u}^\sigma \dot{u}^\nu = -\mu^2,\]

for $v^t$. The argument positive allows the sign of $\mu$ to be changed. true corresponds to time-like geodesics, false to space-like.

This function should rarely be directly called, and instead is invoked by _constrain.

Limitations:

  • currenly pre-supposes static, axis-symmetric metric.
source

Methods

Gradus.TransferFunctionMethod — Type
TransferFunctionMethod()

Computing the underlying relativistic effects using Cunningham transfer functions. This method has a number of limitations related to the assumptions:

  • The disc must be axis-symmetric.
  • The accretion disc must have a velocity structure that depends only on radius.
  • The domain of integration must be in or above the equatorial plane (no false images).

The transfer function approach is often faster and converges to higher numerical accuracy than the other methods. They can also be pre-computed for use in spectral models for fitting.

source
Gradus.BinningMethod — Type
BinningMethod()

Compute the underlying relativistic effects by 'binning' the observer's plane into a number of regions, and tracing a single geodesic for each photon. Conceptually, this is like tracing a single photon for each pixel, however the regions need not be equi-rectangular, and may take on other shapes, see AbstractImagePlane.

This method is slow and computationally expensive, but has the benefit that it makes no assumptions about the model being computed.

It is internally used for testing conceptually more complex integration methods.

source

Image planes

Gradus.PolarPlane — Type
function PolarPlane(
    grid::AbstractImpactParameterGrid;
    Nr = 400,
    Nθ = 100,
    r_min = 1.0,
    r_max = 250.0,
    θ_min = 0.0,
    θ_max = 2π,
)

Divide the image plane into a polar grid centered at $\alpha = 0$ and $\beta = 0$.

source
Gradus.CartesianPlane — Type
function CartesianPlane(
    grid::AbstractImpactParameterGrid;
    Nx = 150,
    Ny = 150,
    x_min = 0.0,
    x_max = 150.0,
    y_min = 0.0,
    y_max = 150.0,
)

Represent the image plane as equi-rectangular regions.

source
Gradus.image_plane — Function
image_plane(plane::AbstractImagePlane)
image_plane(plane::AbstractImagePlane, x::SVector{4})

Return two vectors, representing the $\alpha$ and $\beta$ impact parameters that parameterise the image plane. Each pair of $\alpha$ and $\beta$ must correspond to an unnormalized_areas.

source
Gradus.trajectory_count — Function
trajectory_count(plane::AbstractImagePlane)

Return an integer that counts how many unique geodesics need to be calculated to map the plane. This should be equal to length(first(image_plane(plane))), and is used to pre-allocate buffers.

source
Gradus.unnormalized_areas — Function
unnormalized_areas(plane::AbstractImagePlane)

Return a vector where each element is the area (number) of a given geodesic element on the image plane. For a pixel image plane, each area will be a constant 1. These are used to weight the contributions of each region when calculating observational results.

source