GitHub

Lens collections

The hierarchy Gravity.SingleLens – Gravity.LensPlane – Gravity.LensSystem, as explained before, makes it possible to define complex lensing models efficiently.

Lens Plane

A lens plane is a collection of single lenses at the same redshift. We compute the potential, the deflection, and the Jacobian of a lens plane as the sum of the corresponding quantities for the single lenses belonging to the lens plane. Special care is taken to make the code as much efficient as possible: in particular

  1. all methods acting on a lens plane make use of type-stable and concrete variables only, so that they are highly optimized by the compiler, and

  2. all methods acting on a lens plane do not make any allocation.

These results are obtained using a special data-structure for all lenses in a lens plane, which are organized by the lens type.

Gravity.LensPlane — Type
LensPlane

Structure representing a single lens plane (i.e., a vector of SingleLens'es) at a given redshift z.

Members

  • z: plane's redshift (common to all lenses)
  • lenses: a tuple of vectors of subtypes of SingleLens'es, containing all lenses located at redshift z

Constructors

LensPlane(z, ([l1a, l1b...], [l2a, l2b...]...))
LensPlane(l1, l2, ...)
LensPlane([l1, l2, ...])

where each l1, l2, ... is a concrete SingleLens.

Methods

A LensPlane, being a collection of single lenses at the same redshift, has the same lensing methods of a single lens, namely convergence, potential, deflection, jacobian, deflectionjacobian, and Gravity.redshift.

Note

The particular structure of the lenses field is a design choice to make computations fast. A natural choice would have been to set lenses just as a vector of concrete SingleLens'es. Unfortunately, this way we would be unable to perform type inference on the various methods (potential, deflection...) as this would require a precise knowledge of the type of each element of lenses (unavailable to the compiler). An alternative would have been to define lenses as a simple tuple of concrete SingleLens'es, but this technique would be very expensive computationally in case the number of lenses is more than a handful. The solution is to use for lenses a tuple of vectors, with each vector composed of just one specific kind of concrete SingleLens: since there are a limited number SingleLens types, and since it is anyway unlikely that one would want to use more than a handful different kind of SingleLens'es in a lensing model, this approach works nicely.

source

Lens System

A lens system is an ordered collection of different lens planes at increasing redshifts. All methods acting on a lens system generally require the use of a SourceCoeff, so that the lens equation is targeted to the specific source.

Similarly to a lens plane, methods acting on a lens system are highly optimized, again avoiding non-concrete variables and allocations. However, note that in the current version of Julia (1.8.0 at the time of this writing), the methods use concrete types only if the number of planes in a lens system does not exceed 3: this is generally not a severe limitation, since one would hardly model a lens system using more than one or two planes.

Gravity.LensSystem — Type
LensSystem{T, V, P}

A set of lens planes at different redshifts within a cosmology.

The three type parameters identify the base type of the LensSystem (T, typically, a Float64), the vector type (V), and the type of the tuple for lensplanes (P).

Members

  • lensplanes: a tuple of LensPlane's, with increasing redshifts
  • λ: a LensCoeff structure associated to the lens planes
  • cosmo: the assumed cosmology
  • hubblefactor: a dimensionless real representing the factor cosmology.H₀ / H₀. It is used to scale the timedelay appropriately, in case the Hubble constant reported in the cosmology strucure is different from the true one.

Constructors

LensSystem([plane1, plane2...] [, cosmology])
LensSystem(plane1, plane2...[; cosmo=cosmology])

Build a LensSystem from the given lens planes; optionally one can specify a non-standard cosmology (by default the standard Gravity.default_cosmology is assumed). The last two constructors allow the direct use of SingleLens'es, which are then sorted for their redshifts and put into LensPlane's.

Warning

The constructors are not type-stable, because all of them check the redshifts of the various planes and reorder them if necessary. The only exception is the case of a single plane, as in LensSystem(plane; [cosmo=cosmology]) Alternatively, one can use the basic constructor, by just passing all fields in turn: LensSystem((plane1, plane2...), λ, cosmology).

Note

The lensing coefficients, which are linked to the cosmology, can have a different type than the lensing planes. This is to avoid derivatives of the lensing coefficients with automatic differentiation, with all the related issues.

source
Gravity.timedelay — Function
timedelay(lenssys, ς, γ)
timedelay(lenssys, ς, γ₁, γ₂...)

Compute the time delay associated to a given path γ.

The parameter ς is a SourceCoeff structure holding information on the source redshift and lensing strength; it must be created for the specific lenssys.

Warning

This computation of the time delay requires the full path information. This is either obtained from the fullpath. The procedure used here is less efficient of a direct computation of the timedelay and of the path; it is also less numerically stable.

source
timedelay(lenssys, ς, θ)

Compute the time delay associated to a given point θ.

The parameter ς is a SourceCoeff structure holding information on the source redshift and lensing strength; it must be created for the specific lenssys.

source
Gravity.lensmap — Function
lensmap(lenssys, ς, θ)
lensmap(lenssys, ς, θs; use_threads=false)

Compute the source position β corresponding to an image position θ.

This method applies the multi-plane the lens equation (ray-tracing). lenssys must be a LensSystem structure, with associated SourceCoeff information in ς.

When applied to a vector of points θs, the method returns a vector of source positions βs. This is equivalent to a broadcast call lensmap.(lenssys, ς, θs), but uses a different approach based on lensmap!, which generally results in a slightly faster computation (by approximately a 30% factor), at the espense of more memory allocation. The optional parameter use_threads allows to parallelize the computation using threads (see lensmap! for more details).

Warning

For efficiency reasons, if ς.i == 0 (i.e., if the sources are in front of all lenses), the method for vectors points returns directly θs, without making any copy. Make sure you understand this point if you apply updates to the result of this function.

source
lensmap(lenssys, ς, image)

Compute the source corresponding to an image.

This method applies the multi-plane the lens equation (ray-tracing). lenssys must be a LensSystem structure, with associated SourceCoeff information in ς. image can be any Gravity.SingleSource, such as a Gravity.PointSource or a Gravity.LuminousSource, or any Gravity.ImageFamily.

The returned output will be of the same type of image.

source
lensmap(lenssys, image_plane)

Map an image-plane into a source plane.

This method calls the appropriate lens-mapping routines for all sources in the image-plane and collects the results into a source-plane.

Depending on the value of Gravity.use_threads, the lensing operations for the various sources will be performed using threads.

source
Gravity.lensmap! — Function
lensmap!(sys, ς, xs [, δs]; use_threads=false)

Perform a lens mapping operation in-place for an array of SPoints xs.

The parameter ς is a SourceCoeff structure holding information on the source redshift and lensing strength; it must be created for the specific lens system sys.

The optional parameter δs is an array of SPoints, which will be used to store the deflection vectors. If not provided, δs will be allocated with the same size as xs.

If use_threads is true and length(xs) > 1000, the computation will be parallelized using threads: no threads are used for less than 1000 points, as the overhead of thread creation is generally not worth it in this case.

The method returns the array xs with the source positions.

See also: lensmap.

source
Gravity.timedelaylensmap — Function
timedelaylensmap(lenssys, ς, θ)

Compute together both the time delay and the lensmap at a given point θ.

The parameter ς is a SourceCoeff structure holding information on the source redshift and lensing strength; it must be created for the specific lenssys.

This is generally faster than performing the separate computation of the lensing path (using fullpath) and of the timedelay.

source
Gravity.fullpath — Function
fullpath(lenssys, ς, θ)

Compute the entire path associated to a image position θ.

This procedure, given a LensSystem lenssys, a SourceCoeff ς and an image position θ, returns an array of tuples in the format (xᵢ, fᵢ) which can be used to compute the source-plane position for any source distance. In particular, calling d the source comoving transverse distance, one has

β = xᵢ + fᵢ / d

Example

```julia-repl julia> ls = rand(LensSystem);

julia> sc = SourceCoeff(0.78, ls)

source
Gravity.fullpath! — Function
fullpath!(γ, lenssys, ς, θ)

In-place version of fullpath.

The parameter γ must be a pre-allocated vector of points with size at least ς.i + 1.

source
Gravity.xfullpath — Function
xfullpath(lenssys, x)

Compute the entire path associated to a image position x.

This procedure, given a LensSystem lenssys and an image position x, returns a tuples ((x₁, δ₁), (x₂, δ₂), ..., (xₙ, δₙ)) which can be used to compute the source-plane position for any source distance. In particular, one has

β = x[ς.i+1] - δ[ς.i+1] * ς.ρᵢ

where ς is the SourceCoeff associated to the source redshift. This operation can be also performed with the help of xfullpath_map.

See also: xfulldistortion.

source
xfullpath(lenssys, xs)

Compute the entire path associated to an array of image positions xs.

This procedure, given a LensSystem lenssys and an array of image position xs, returns a tuple (ys, δs), where both ys and δs are two arrays with dimensions higher (by unity) than xs. The last axes of ys and δs encode the position of the light ray at the various lens planes of lenssys.

In particular, one can compute the source positions for xs as

βs = xs[..., ς.i+1] - δs[..., ς.i+1] * ς.ρᵢ

where ς is the SourceCoeff associated to the source redshift. This operation can be also performed with the help of xfullpath_map.

See also: xfullpath!, xfulldistortion.

source
Gravity.xfullpath! — Function
xfullpath!(xp, δp, lenssys, xs)

In-place version of xfullpath.

The two arrays xp and δp must be pre-allocated and must have one dimension more than xs to store all the intercepts at the lens planes. For example, they can be allocated using the expression

xp = similar(xs, (axes(xs)..., length(lenssys.lensplanes) + 1))
δp = similar(xp)
source
Gravity.xfullpath_map — Function
xfullpath_map(path, ς)
xfullpath_map(xp, δp, ς)
xfullpath_map((xp, δp), ς)

Compute the source position associated to a given extended full path.

This function takes a path obtained with xfullpath, or alternatively the two arrays xp and δp obtained with the same function (with array arguments), togehter with a SourceCoeff ς, and computes the lensmap. Therefore, one has always

@assert xfullpath_map(xfullpath(sys, x), ς) ≈ lensmap(sys, ς, x)
source
Gravity.lensmapdistortion — Function
lensmapdistortion(lenssys, ς, θ)

Compute the image position and its associated distortion for an image position θ.

lenssys must be a LensSystem structure, with associated SourceCoeff information in ς. The result will be a SPointDiff, i.e. a static 6-vector: the first two elements represent the source position β, and the last four the four components of its distortion 𝒜 = ∂β/∂θ.

source
Gravity.xfulldistortion — Function
xfulldistortion(lenssys, x)

Compute the entire path and distortion associated to a image position x.

This procedure, given a LensSystem lenssys and an image position x, returns the tuple ((x₁, δ₁, 𝒜₁, 𝒟₁), (x₂, δ₂, 𝒜ₙ, 𝒟ₙ), ..., (xₙ, δₙ, 𝒜ₙ, 𝒟ₙ)). This can be used to compute the source-plane position and distortion for any source distance. In particular, one has

β = x[ς.i+1] - δ[ς.i+1] * ς.ρᵢ
𝒜 = 𝒜[ς.i+1] - 𝒟[ς.i+1] * ς.ρᵢ

where ς is the SourceCoeff associated to the source.

See also: xfullpath.

source
Gravity.magnification — Function
magnification(lenssys, ς, θ)
magnification(lenssys, z, θ)
magnification(lenssys, ς)
magnification(lenssys, z)

Compute the inverse of the magnification of a lens system at the position θ.

In last two versions, return a closure that, when called with a point θ as an argument, will compute the magnification. This form is useful as an argument to Gravity.refine_grid_isocontour!.

The inverse magnification is calculated as det(distortion(lenssys, ς, θ)), with ς = SourceCoeff(z, lenssys).

source