Geodesic integration strategies

velocity_to_sky_angles(
    m::AbstractMetric,
    x_source::SVector{4},
    v_source::SVector{4},
    v::SVector{4},
)

Map the velocity vector v to angles on the local sky $(\theta, \phi)$ of a source at position x_source moving with velocity v_source. Returns a tuple with (θ, ϕ).

See sky_angles_to_velocity for the inverse operation.

sky_angles_to_velocity(
    m::AbstractMetric,
    x_source::SVector{4},
    v_source::SVector{4},
    θ,
    ϕ
)

Convert angles on the local sky $(\theta, \phi)$ at some point x_source to a velocity vector. Also applies the Lorentz transformation due to the velocity v_source of x_source. Returns a four-velocity vector (SVector{4}) that is unconstrained. Use constrain_all or similar to normlise the four-velocity.

See velocity_to_sky_angles for the inverse operation.

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.

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.

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.

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.

AbstractImagePlane

An plane abstraction used to represent the observer's image plane. These are particularly relevant for BinningMethod computations.

Concerete implementations include:

• PolarPlane
• CartesianPlane

An AbstractImagePlane must implement the following functions:

• image_plane
• trajectory_count
• unnormalized_areas

function PolarPlane(
    grid::Abstract2DGrid;
    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$.

function CartesianPlane(
    grid::Abstract2DGrid;
    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.

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.

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.

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.