Repository

Point functions

Point functions operate on AbstractGeodesicPoint, and are used to represent calculations that can be applied to any single geodesic. They usually map to some sort of physical value (an energy shift, a proper time, and so on), but can also be used to filter geodesics based on certain predicates.

PointFunctions are mainly indented for use with rendergeodesics, as short-hands for expressing how a pixel in the image should be calculated.

There are many already defined in Gradus.jl (e.g. ConstPointFunctions.redshift) but it's easy to define your own.

For example, a point function which calculates the angle from the zenith in the local coordinates at the endpoint of a geodesic. First we'll use a pre-made point function which filters those geodesics that hit some accretion geometry:

using Gradus, Makie, CairoMakie

m = KerrMetric(1.0, 0.998)
d = ThinDisc(0.0, 50.0)
x_obs = SVector{4}(0.0, 1e4, deg2rad(80), 0.0)

# use a pre-made point function
pf = ConstPointFunctions.filter_intersected()

# apply the point function to colour the image
a, b, img = rendergeodesics(m, x_obs, d, 2x_obs[2]; verbose = true, pf = pf)

fig = Figure()
ax = Axis(fig[1,1])
heatmap!(ax, a, b, img')
fig

Now we'll write a point function that calculates the angle on the local sky. All point functions are given three things: the metric, the GeodesicPoint, and the value of the affine parameter.

Note

The affine parameter argument is now redundant, since that information is included in the GeodesicPoint but exists for backwards compatibility.

We can define our function like this:

source_velocity = SVector(1.0, 0.0, 0.0, 0.0)
angle_pf = PointFunction(
    (m, gp, lambda) ->
        # Only need the first element, as we're only interested in theta
        Gradus.velocity_to_sky_angles(m, gp.x, source_velocity, gp.v)[1]
)

This can be composed together with the filter point function to first apply the filter and then our custom point function:

new_pf = angle_pf ∘ ConstPointFunctions.filter_intersected()
# Apply the point function to colour the image
a, b, img = rendergeodesics(m, x_obs, d, 2x_obs[2]; verbose = true, pf = new_pf)

fig = Figure()
ax = Axis(fig[1,1])
heatmap!(ax, a, b, img')
fig

We can include the effects of the motion of the disc also. Assuming the disc has Keplerian velocities, the point function could be amened to use the local four-velocity of the disc instead of a stationary point for the angle transformation:

function to_local_angle(m, gp, lambda)
    v_disc = CircularOrbits.fourvelocity(m, gp.x[2])
    Gradus.velocity_to_sky_angles(m, gp.x, v_disc, gp.v)[1]
end

# Wrap our function as before
comoving_angle_pf = PointFunction(to_local_angle)

# And compose and trace!
comoving_pf = comoving_angle_pf ∘ ConstPointFunctions.filter_intersected()
a, b, img = rendergeodesics(m, x_obs, d, 2x_obs[2]; verbose = true, pf = comoving_pf)

fig = Figure()
ax = Axis(fig[1,1])
heatmap!(ax, a, b, img')
fig

The effects of special relativistic beaming are here clear to see!

Point function documentation

Gradus.AbstractGeodesicPoint — Type
abstract type AbstractGeodesicPoint{T}

Supertype for geodesic points, used to store information about specific points along geodesic trajectories.

Note

Currently limited to storing the start and endpoint of any given trajectory. To keep the full geodesic path, it is encouraged to use the SciMLBase.AbstractODESolution directly.

Must minimally have the same fields as GeodesicPoint. Examples include Gradus.FirstOrderGeodesicPoint.

source
Gradus.GeodesicPoint — Type
struct GeodesicPoint <: AbstractGeodesicPoint

A structure that represents the start and end point along a geodesic. To obtain a GeodesicPoint, use unpack_solution or unpack_solution_full on the result of a tracegeodesics.

Fields:

  • status::StatusCodes.T: A StatusCodes code of the integrator for this geodesic.
  • λ_min::T: Start affine time
  • λ_max::T: End affine time
  • x_init::SVector{4,T}: Start four position
  • x::SVector{4,T}: End four position
  • v_init::SVector{4,T}: Start four velocity
  • v::SVector{4,T}: End four velocity
  • aux::A: Auxillary values (polarisation, intensities, etc.)

Example

sol = tracegeodesics(m, x, v, 100.0)
gp = unpack_solution(sol)
gp.status
source
Gradus.PointFunction — Type
struct PointFunction <: AbstractPointFunction
PointFunction(func)

Point functions are functions that are used to calculate physical parameters from geodesic integrations, and to compose more complex models. A number of default and utility PointFunction are defined in Gradus.ConstPointFunctions.

Principally, point functions return a single value per geodesic, and are used to fill rendered images with values, e.g. colouring redshift.

Point functions may be instantiated by wrapping a function with the following signature

function func(m::AbstractMetric{T}, gp::AbstractGeodesicPoint, max_time::T; kwargs...)::T where {T}
    # ...
end

pf = PointFunction(func)
  • The AbstractMetric argument may be used to dispatch for different metrics.
  • gp is an AbstractGeodesicPoint corresponding to a given geodesic.
  • The max_time parameter is the maximum integration time used to integrate the geodesics. This may be useful when trying to determine whether a geodesic terminated early or not.

They may be invoked by invoking the instance

result = pf(m, gp, max_time)
Note

As of version 0.1.0, the kwargs parameter is reserved only for passing optional results when chaining multiple point functions (see below). This is subject to revision and breaking changes in future versions.

Multiple AbstractPointFunction may be chained together using the ∘ operator, and are evaluated from right to left

pf3 = pf2 ∘ pf1

This may be useful for constructing filters using FilterPointFunction. When used with two PointFunction objects, the output of the previous PointFunction is passed to the next via the value keyword argument.

source
Gradus.FilterPointFunction — Type
struct FilterPointFunction <: AbstractPointFunction
FilterPointFunction(func, default_value)

Point functions used to filter geodesics. They may be constructed with

function func(m::AbstractMetric{T}, gp::AbstractGeodesicPoint, max_time::T; kwargs...)::Bool where {T}
    # ... return Bool
end

fpf = FilterPointFunction(func, NaN64)

The second argument to the constructor is the default value, given to the pixel if the boolean condition of func is false.

Example

A filter for geodesics within a certain radius, used to only calculate redshift within 10 $\r_\text{g}$

func(m, gp, max_time) = gp.u[2] < 10.0
pf = ConstPointFunctions.redshift(m, u) ∘ FilterPointFunction(func, NaN64)
source

Pre-defined point functions

Gradus.ConstPointFunctions.shadow — Function
shadow(m::AbstractMetric, gp::AbstractGeodesicPoint, max_time)

A PointFunction which colours the shadow of the black hole for any disc-less render. Equivalent to ConstPointFunctions.affine_time ∘ ConstPointFunctions.filter_early_term.

source
Gradus.ConstPointFunctions.redshift — Function
redshift(m::AbstractMetric, x_obs)

Returns a PointFunction that calculates the redshift of a medium as observed by an observer at x_obs.

Calculate the analytic redshift at a given geodesic point, assuming equatorial, geometrically thin accretion disc. Implementation depends on the metric type, with non-Kerr metrics using an interpolation scheme within the plunging region via interpolate_redshift.

Calls redshift_function to dispatch different implementations.

source