API

Base.:+(m1::AbstractModel, m2::AbstractModel)

Combine two models to create a composite AddModel. This adds two models pointwise, i.e.

julia> visibility(m1 + m2, 1.0, 1.0) == visibility(m1, 1.0, 1.0) + visibility(m2, 1.0, 1.0)
m1 = Gaussian()

added(m1::AbstractModel, m2::AbstractModel)

Combine two models to create a composite AddModel. This adds two models pointwise, i.e.

julia> visibility(added(m1, m2), 1.0, 1.0) == visibility(m1, 1.0, 1.0) + visibility(m2, 1.0, 1.0)
m1 = Gaussian()

convolved(m1::AbstractModel, m2::AbstractModel)

Convolve two models to create a composite ConvolvedModel.

julia> convolved(m1, m2)
m1 = Ring()

components(m::AbstractModel)

Returns the model components for a composite model. This will return a Tuple with all the models you have constructed.

Example

julia> components(m)
m = Gaussian() + Disk()

smoothed(m::AbstractModel, σ::Number)

Smooths a model m with a Gaussian kernel with standard deviation σ.

Notes

This uses convolved to created the model, i.e.

julia> convolved(m1, m2) == smoothed(m1, 1.0)
m1 = Disk()

abstract type CompositeModel{M1, M2} <: ComradeBase.AbstractModel

Abstract type that denotes a composite model. Where we have combined two models together.

Implementation

Any implementation of a composite type must define the following methods:

• visibility_point
• uv_combinator
• imanalytic
• visanalytic
• ComradeBase.intensity_point if model intensity is IsAnalytic
• intensitymap_numeric! if model intensity is NotAnalytic
• intensitymap_numeric if model intensity is NotAnalytic
• flux
• radialextent

In addition there are additional optional methods a person can define if needed:

• intensitymap_analytic! if model intensity is IsAnalytic (optional)
• intensitymap_analytic if model intensity is IsAnalytic (optional)
• visibilitymap_analytic if visanalytic is IsAnalytic (optional)
• visibilitymap_numeric if visanalytic is Not Analytic (optional)
• visibilitymap_analytic! if visanalytic is IsAnalytic (optional)
• visibilitymap_numeric! if visanalytic is Not Analytic (optional)

struct AddModel{T1, T2} <: VLBISkyModels.CompositeModel{T1, T2}

Pointwise addition of two models in the image and visibility domain. An end user should instead call added or Base.+ when constructing a model

Example

julia> m1 = Disk() + Gaussian()
julia> m2 = added(Disk(), Gaussian()) + Ring()

struct ConvolvedModel{M1, M2} <: VLBISkyModels.CompositeModel{M1, M2}

Pointwise addition of two models in the image and visibility domain. An end user should instead call convolved. Also see smoothed(m, σ) for a simplified function that convolves a model m with a Gaussian with standard deviation σ.

abstract type GeometricModel{T} <: ComradeBase.AbstractModel

A type that defines it is a geometric model. These are usually primitive models, and are usually analytic in Fourier and the image domain. As a result a user only needs to implement the following methods

• visibility_point
• intensity_point
• radialextent

Note that if the geometric model isn't analytic then the usual methods listed in ComradeBase.AbstractModel for non-analytic models need to be implemented.

struct ConcordanceCrescent{T} <: VLBISkyModels.GeometricModel{T}

Creates the ConcordanceCrescent model, i.e. a flat-top crescent with a displacment and a slash and shadow depth. Note this creates a crescent with unit flux. If you want a different flux please use the renomed modifier.

Fields

• router: Outer radius of the crescent

• rinner: Inner radius of the crescent (i.e. inside this radius there is a hole)

• shift: Displacment of the inner disk radius

• slash: Strength of the linear slash. Note that s∈[0.0,1.0] to ensure positivity in the image.

Notes

Unlike the Gaussian and Disk models this does not create the unit version. In fact, this model could have been created using the Disk and primitives by using VLBISkyModels's model composition functionality.

Creates a Kamruddin and Dexter crescent model. This works by composing two disk models together.

Arguments

• router: The radius of the outer disk
• rinner: The radius of the inner disk
• shift: How much the inner disk radius is shifted (positive is to the right)
• floor: The floor of the inner disk 0 means the inner intensity is zero and 1 means it is a large disk.

Disk{T}() where {T}

Tophat disk geometrical model, i.e. the intensity profile

\[ I(x,y) = \begin{cases} \pi^{-1} & x^2+y^2 < 1 \\ 0 & x^2+y^2 \geq 0 \end{cases}\]

i.e. a unit radius and unit flux disk.

By default if T isn't given, Disk defaults to Float64

SlashedDisk{T}(slash::T) where {T}

Tophat disk geometrical model, i.e. the intensity profile

\[ I(x,y) = \begin{cases} \pi^{-1} & x^2+y^2 < 1 \\ 0 & x^2+y^2 \geq 0 \end{cases}\]

i.e. a unit radius and unit flux disk.

By default if T isn't given, Disk defaults to Float64

struct ExtendedRing{T} <: VLBISkyModels.GeometricModel{T}

A symmetric extended ring whose radial profile follows an inverse gamma distributions.

The formula in the image domain is given by

I(r,θ) = βᵅrᵅ⁻²exp(-β/r)/2πΓ(α)

where α = shape and β = shape+1

Note

We mainly use this as an example of a non-analytic Fourier transform (although it has a complicated expression)

Fields

• shape: shape of the radial distribution

Note that if T isn't specified at construction then it defaults to Float64.

struct Gaussian{T} <: VLBISkyModels.GeometricModel{T}

Gaussian with unit standard deviation and flux.

By default if T isn't given, Gaussian defaults to Float64

struct MRing{T, V<:Union{AbstractArray{T, 1}, Tuple{Vararg{T, N}} where {N, T}}} <: VLBISkyModels.GeometricModel{T}

m-ring geometric model. This is a infinitely thin unit flux delta ring whose angular structure is given by a Fourier expansion. That is,

I(r,θ) = (2π)⁻¹δ(r-1)∑ₙ(αₙcos(nθ) - βₙsin(nθ))

The N in the type defines the order of the Fourier expansion.

Fields

• α: Real Fourier mode coefficients

• β: Imaginary Fourier mode coefficients

struct Ring{T} <: VLBISkyModels.GeometricModel{T}

A infinitely thin ring model, whose expression in the image domain is I(r,θ) = δ(r - 1)/2π i.e. a unit radius and flux delta ring.

By default if T isn't given, Gaussian defaults to Float64

struct ParabolicSegment{T} <: VLBISkyModels.GeometricModel{T}

A infinitely thin parabolic segment in the image domain. The segment is centered at zero, with roots ±1 and a yintercept of 1.

Note that if T isn't specified at construction then it defaults to Float64.

struct ZeroModel{T} <: ComradeBase.AbstractModel

Defines a model that is empty that is it return zero for everything.

Notes

This returns 0 by using FillArrays so everything should be non-allocating

MultiComponentModel(beam::AbstractModel, fluxes::AbstractVector, x::AbstractVector, y::AbstractVector)

Build a model with a base model type beam where fluxes, x, y corresond to the flux, and positions of the components. This can be used to construct clean like models.

struct PolarizedModel{TI, TQ, TU, TV} <: ComradeBase.AbstractPolarizedModel

Wrapped model for a polarized model. This uses the stokes representation of the image.

Fields

• I: Stokes I model

• Q: Stokes Q Model

• U: Stokes U Model

• V: Stokes V Model

abstract type AbstractImageTemplate <: ComradeBase.AbstractModel

An abstract ComradeBase.AbstractModel that serves as the parent type of all the VIDA templates implemented in this repo.

By default this model assumes the following ComradeBase traits

ComradeBase.visanalytic(::Type{<:AbstractImageTemplate}) = NoAnalytic()
ComradeBase.imanalytic(::Type{<:AbstractImageTemplate})  = IsAnalytic()
ComradeBase.ispolarized(::Type{<:AbstractImageTemplate}) = NotPolarized()

As a result if a user wishes the implement their own subtype (e.g., MyTemplate) of AbstratImageTemplate they will need to implement the following methods

• ComradeBase.intensity_point(m::MyTemplate, p): which computes the potentially unormalized brightness of the template at the point p.
• ComradeBase.radialextent(m::MyTemplate): which computes the rough radial extent of the model m.

For more information about the total interface see VLBISkyModels.jl

struct RingTemplate{R<:VLBISkyModels.AbstractRadial, A<:VLBISkyModels.AbstractAzimuthal} <: VLBISkyModels.AbstractImageTemplate

A flexible ring template that forms a ring by taking the product of a radial and azimuthal brightness profile.

A list of radial profiles is given by subtypes(AbstractRadial)

A list of azimuthal profiles is given by subtypes(AbstractAzimuthal)

Examples

julia> rad = RadialGaussian(0.1)
julia> azi = AzimuthalUniform()
julia> ring = modify(RingTemplate(rad, azi), Stretch(10.0), Shift(1.0, 2.0))

Fields

• radial: Radial profile of the ring

• azimuthal: Azimuthal profile of the ring

RadialGaussian(σ)

Create a profile that is a radial gaussian ring with radius unity and standard deviation σ.

Notes

This is usually couple with a azimuthal profile to create a general ring template

julia> rad = RadialGaussian(0.1)
julia> azi = AzimuthalUniform()
julia> t = RingTemplate(rad, azi)

Arguments

• σ: The standard deviation for the Gaussian ring.

RadialDblPower(αinner, αouter)

Radial profile that given by a double power-law with the function form

    r^αinner*inv(1+ r^(αinner + αouter + 1))

Notes

This is usually couple with a azimuthal profile to create a general ring template

julia> rad = RadialDblPower(3.0, 3.0)
julia> azi = AzimuthalUniform()
julia> t = RingTemplate(rad, azi)

Arguments

• αinner the power law index for r<1.
• αouter the power law index for r≥1.

RadialTruncExp(σ)

Radial profile that has a exponential profile up to a cutoff radius of 1 where for r<1 the profile is identically zero.

This evalutes to

    exp(-(r-1)/σ)

Notes

This is usually couple with a azimuthal profile to create a general ring template

julia> rad = RadialTruncExp(2.0)
julia> azi = AzimuthalUniform()
julia> t = RingTemplate(rad, azi)

Arguments

• σ: Exponential inverse fall off parameter.

RadialJohnsonSU(αinner, αouter)

Radial profile that given by positive unit radius Johnson SU distributions with the functional form

    exp[-1/2(γ + arcsinh((r - μ)/σ)]^2*inv(1 + (r - 1)² + σ²)

where controls σ ≥ 0 the width, and γ the asymmetry.

Notes

This is usually couple with a azimuthal profile to create a general ring template

julia> rad = RadialJohnsonSU(1/2, -3/2)
julia> azi = AzimuthalUniform()
julia> t = RingTemplate(rad, azi)

Arguments

• σ : the width of the ring
• γ : the asymmetry of the ring

AzimuthalUniform()
AzimuthalUniform{T}()

A azimuthal profile that is uniform for all angles.

Notes

This is usually couple with a radial profile to create a general ring template

julia> rad = RadialDblPower(3.0, 3.0)
julia> azi = AzimuthalUniform() # Defaults to Float64
julia> t = RingTemplate(rad, azi)

AzimuthalCosine(s::NTuple{N,T}, ξ::NTuple{N, T}) where {N, T}

A azimuthal profile that is given by a cosine expansion of order N. The expansion is given by

    1 - ∑ₙ sₙcos(nϕ - ξₙ)

Notes

This is usually couple with a radial profile to create a general ring template

julia> rad = RadialDblPower(3.0, 3.0)
julia> azi = AzimuthalCosine((0.5, 0.2), (0.0, π/4)) # Defaults to Float64
julia> t = RingTemplate(rad, azi)

Arguments

• s : amplitudes of the N order cosine expansion of the azimuthal brightness
• ξs: phase of the N order cosine expansion of the azimuthal brightness

GaussianRing(σ)

Implements the gaussian ring template.

Notes

This is a convienence constructor that uses RingTemplate under the hood. To create this function your self do

RingTemplate(RadialGaussian(σ), AzimuthalUniform())

Arguments

• σ : standard deviation of the Gaussian ring

GaussianRing(r0, σ, x0, y0)

Implements the gaussian ring template.

Notes

This is a convienence constructor that uses RingTemplate under the hood. To create this function your self do

modify(GaussianRing(σ/r0), Stretch(r0), Shift(x0, y0))

Arguments

• r0: radius of the ring
• σ : standard deviation of the Gaussian ring
• x0: location of the ring center horizontally
• y0: location of the ring center vertically

SlashedGaussianRing(σ, s)

Implements the slashed gaussian ring template, that uses a cosine to implement the slash and has a brightness position angle of 0 degrees east of north.

Arguments

• σ : standard deviation of the Gaussian ring
• s : Slash amplitude. 0 means no slash, and 1 is maximal.

Notes

This is a convienence constructor that uses RingTemplate under the hood. To create this function your self do

RingTemplate(RadialGaussian(σ), AzimuthalCosine((s,), (zero(s),)))

SlashedGaussianRing(r0, σ, s, ξ, x0, y0)

Implements the slashed gaussian ring template, that uses a cosine to implement the slash.

Arguments

• r0: radius of the ring
• σ : standard deviation of the Gaussian ring
• s : Slash amplitude. 0 means no slash, and 1 is maximal.
• ξs: azimuthal peak brightness P.A. measured east of north
• x0: location of the ring center horizontally
• y0: location of the ring center vertically

Notes

This is a convienence constructor that uses RingTemplate under the hood. To create this function your self do

modify(SlashedGaussianRing(σ/r0, s), Stretch(r0), Rotate(ξ), Shift(x0, y0))

EllipticalGaussianRing(r0, σ, τ, ξτ, x0, y0)

Implements the elliptical gaussian ring template with semi-minor axis a and semi-major axis b.

Arguments

• r0: geometric mean radius (√ab) of the ring
• σ : standard deviation of the Gaussian ring
• τ : asymmetry of the ring τ = 1-b/a
• ξτ: semi-major axis measured east of north
• x0: location of the ring center horizontally
• y0: location of the ring center vertically

Notes

This is a convienence constructor that uses RingTemplate under the hood. To create this function your self do

julia modify(GaussianRing(σ/r0), Stretch(r0*sqrt(1-τ), r0/sqrt(1-τ)), Rotate(ξτ), Shift(x0, y0))

EllipticalSlashedGaussianRing(r0, σ, τ, ξτ, s, ξs, x0, y0)

Combination of the elliptical and slashed gaussian ring.

Details

Adds ellipticity to the ring. The radius r0 of the ring is now defined as the geometric mean of the semi-major (a) and semi-minor (b) axis lengths i.e.

r0 = √(a*b).

The ellipticity τ is given by τ = 1-b/a.

The brightness asymmetry uses a cosine to implement the slash.

Arguments

• r0: geometric mean radius (√ab) of the Gaussian ring
• σ : standard deviation of the Gaussian ring
• τ : asymmetry of the Gaussian ring τ = 1-b/a
• ξτ: semi-major axis measured east of north
• s : Slash amplitude. 0 means no slash, and 1 is maximal.
• ξs: azimuthal peak brightness P.A. measured east of north
• x0: location of the ring center horizontally
• y0: location of the ring center vertically

Notes

This is a convienence constructor that uses RingTemplate under the hood. To create this function your self we define

julia modify(SlashedGaussianRing(σ/r0, s), Rotate(ξs-ξτ), Stretch(r0*sqrt(1-τ), r0/sqrt(1-τ)), Rotate(ξτ), Shift(x0, y0) )

struct CosineRing{N, M, T} <: VLBISkyModels.AbstractImageTemplate

A flexible symmetric ring model. The thickness is modeled as a cosine expansion with N terms and the slash by a expansion with M terms.

Details

The ring is forced to be symmetric for a significant speed boost over CosineRing. The thickness of the ring is modeled by a cosine expansion in azimuthal angle. N specifies the number of cosine modes to fit, where the first mode is the constant thickness portion and so has no corresponding angle. The slash is modeled as a separate cosine expansion, with M terms. Here the zero order term is forced to be unity, so M defines the M additional terms.

CosineRingwFloor(r0, σ0, σ, ξσ, s, ξs, floor, x0, y0)

A Cosine ring with a GaussDisk at the center of the ring that matches the r0 and σ0 of the CosineRing.

This is a convienence constructor of the more basic VLBISkyModel image constructors. It is equivalent to

julia> CosineRing(r0, σ0, σ, ξσ, s, ξs, x0, y0) + floor*modify(GaussDisk(σ0/r0), Stretch(r0), Shift(x0, y0))

Details

Adds ellipticity to the ring. The radius r0 of the ring is now defined as the geometric mean of the semi-major (a) and semi-minor (b) axis lengths i.e.

r0 = √(a*b).

The ellipticity τ is given by τ = 1-b/a.

Arguments

• r0: geometric mean radius (√ab) of the ring
• σ0: standard deviation of the Gaussian ring
• σ : amplitudes of the N order cosine expansion of the ring thickness
• ξσ: phase of the N order cosine expansion of the ring thickness
• s : amplitudes of the M order cosine expansion of the azimuthal brightness
• ξs: phase of the M order cosine expansion of the azimuthal brightness
• floor: The flux of the center Gaussian. This is relative to the CosineRing.
• x0: location of the ring center horizontally
• y0: location of the ring center vertically

CosineRingwGFloor(r0, σ0, σ, ξσ, s, ξs, floor, σG, x0, y0)

A Cosine ring with a Gaussian blob at the center of the ring. This is a convienence constructor of the more basic VLBISkyModel image constructors. It is equivalent to

julia> CosineRing(r0, σ0, σ, ξσ, s, ξs, x0, y0) + floor*modify(Gaussian(), Stretch(σG), Shift(x0, y0))

Details

Adds ellipticity to the ring. The radius r0 of the ring is now defined as the geometric mean of the semi-major (a) and semi-minor (b) axis lengths i.e.

r0 = √(a*b).

The ellipticity τ is given by τ = 1-b/a.

Arguments

• r0: geometric mean radius (√ab) of the ring
• σ0: standard deviation of the Gaussian ring
• σ : amplitudes of the N order cosine expansion of the ring thickness
• ξσ: phase of the N order cosine expansion of the ring thickness
• s : amplitudes of the M order cosine expansion of the azimuthal brightness
• ξs: phase of the M order cosine expansion of the azimuthal brightness
• floor: The flux of the center Gaussian. This is relative to the CosineRing.
• σG: The size of the central Gaussian.
• x0: location of the ring center horizontally
• y0: location of the ring center vertically

EllipticalCosineRing(r0, σ0, σ, ξσ, τ, ξτ, s, ξs, x0, y0)

An Elliptical Cosine ring. This is a convienence constructor of the more basic VLBISkyModel image constructors. It is equivalent to

julia> modify(CosineRing(σ0/r0, σ/r0, ξσ, s, ξs .- ξτ),
            Stretch(r0*sqrt(1-τ), r0/sqrt(1-τ)),
            Rotate(ξτ),
            Shift(x0, y0)
        )

Details

Adds ellipticity to the ring. The radius r0 of the ring is now defined as the geometric mean of the semi-major (a) and semi-minor (b) axis lengths i.e.

r0 = √(a*b).

The ellipticity τ is given by τ = 1-b/a.

Arguments

• r0: geometric mean radius (√ab) of the ring
• σ0: standard deviation of the Gaussian ring
• σ : amplitudes of the N order cosine expansion of the ring thickness
• ξσ: phase of the N order cosine expansion of the ring thickness
• τ : asymmetry of the ring τ = 1-b/a
• ξτ: semi-major axis measured east of north
• s : amplitudes of the M order cosine expansion of the azimuthal brightness
• ξs: phase of the M order cosine expansion of the azimuthal brightness
• x0: horizontal location of the ring center
• y0: vertical location of the ring center vertically

struct LogSpiral{T<:Real} <: VLBISkyModels.AbstractImageTemplate

Template type for a logarithmic spiral segment

Fields

• κ: Unit curvature of the logarithmic spiral

• σ: thickness of the Gaussian spiral arm

• δϕ: Azimuthal extent of the spiral arm

struct Constant{T} <: VLBISkyModels.AbstractImageTemplate

An constant template.

Details

Defines an image that just has constant flux. This is very useful for soaking up low levels of flux in image reconstructions that can bias the results.

Since images or normalized to unity, this means the Constant template has no additional parameters.

Fields

• scale: Sets the angular scale of the image. This is usually equal to the image FOV

struct GaussDisk{T} <: VLBISkyModels.AbstractImageTemplate

A smoothed disk model

Details

Defines a template for an image that has a smoothed disk model.

Fields

• α: Disk edge standard deviation

ContinuousImage{A<:IntensityMap, P} <: AbstractModel
ContinuousImage(img::Intensitymap, kernel)

The basic continuous image model for VLBISkyModels. This expects a IntensityMap style object as its imag as well as a image kernel or pulse that allows you to evaluate the image at any image and visibility location. The image model is

I(x,y) = ∑ᵢ Iᵢⱼ κ(x-xᵢ, y-yᵢ)

where Iᵢⱼ are the flux densities of the image img and κ is the intensity function for the kernel.

Pulse

Pixel response function for a radio image model. This makes a discrete sampling continuous by picking a certain smoothing kernel for the image.

Notes

To see the implemented Pulses please use the subtypes function i.e. subtypes(Pulse)

struct DeltaPulse{T} <: VLBISkyModels.Pulse

A dirac comb pulse function. This means the image is just the dicrete Fourier transform

$(TYPEDEF)

Uses the basis spline (BSpline) kernel of order N. These are the kernel that come from recursively convolving the tophat kernel

\[ B_0(x) = \begin{cases} 1 & |x| < 1 \\ 0 & otherwise \end{cases}\]

N times.

Notes

BSpline kernels have a number of nice properties:

1. Simple frequency response $\sinc(u/2)^N$
2. preserve total intensity

For N>1 these kernels aren't actually interpolation kernels however, this doesn't matter for us.

Currently only the 0,1,3 order kernels are implemented.

struct RaisedCosinePulse{T} <: VLBISkyModels.Pulse

RaisedCosinePulse()
RaisedCosinePulse(rolloff)

Raised cosine pulse function. This tends to be a very flat response, where the roll off controls the speed of decay. By default we set rolloff = 0.5.

struct BicubicPulse{T} <: VLBISkyModels.Pulse

BicubicPulse(b = 0.5)

The bicubic pulse for imaging. This pulse tends to have a flat spectrum but for most values of b can produce negative intensities in an image. This is the pulse used in Broderick et al. 2020.

Butterworth{N}()
Butterworth{N, T}()

Construct a model that corresponds to the Butterworth filter of order N. The type of the output is given by T and if not given defaults to Float64

FourierDualDomain(imgdomain::AbstractSingleDomain, visdomain::AbstractSingleDomain, algorithm)

Constructs a set of grids that live in the image and visibility domains. The transformation between the grids is specified by the algorithm which is a subtype of VLBISkyModels.FourierTransform.

Arguments

• imgdomain: The image domain grid
• visdomain: The visibility domain grid
• algorithm: The Fourier transform algorithm to use see subtypes(VLBISkyModels.FourierTransform) for a list

FourierDualDomain(imgdomain::AbstractRectiGrid, alg::FFTAlg)

Constructs a FourierDualDomain that uses the FFT algorithm to compute the transformation. For this no visibilty domain is specified since we assume it is the default grid from the FFT with padding specified in FFTAlg.

Arguments

• imgdomain: The image domain grid
• alg: The FFT algorithm to use

xygrid(grid::AbstractRectiGrid)

Converts from a visi1bility domain recti-grid with spatial dimensions (U,V) to a image domain grid with spatial dimensions (X,Y). Note the other dimensions are not changed.

For the inverse see uvgrid

uvgrid(grid::AbstractRectiGrid)

Converts from a image domain recti-grid with spatial dimension (X,Y) to a Fourier or visibility domain grid with spatial dimensions (U,V). Note the other dimensions are not changed.

For the inverse see xygrid

uviterator(nx, dx, ny dy)

Construct the u,v iterators for the Fourier transform of the image with pixel sizes dx, dy and number of pixels nx, ny

If you are extending Fourier transform stuff please use these functions to ensure that the centroid is being properly computed.

DFTAlg

Uses a discrete fourier transform. This is not very efficient for larger images. In those cases NFFTAlg or FFTAlg are more reasonable. For small images this is a reasonable choice especially since it's easy to define derivatives.

struct FFTAlg{T} <: VLBISkyModels.FourierTransform

Use an FFT to compute the approximate numerical visibilitymap of a model. For a DTFT see DFTAlg or for an NFFT NFFTAlg

Fields

• padfac: The amount to pad the image by. Note we actually round up to use small prime factor.

• flags: FFTW flags or wisdom for the transformation. The default is FFTW.ESTIMATE, but you can use FFTW.MEASURE for better performance if you plan on evaluating the sample FFT multiple times.

NFFTAlg

Uses a non-uniform FFT to compute the visibilitymap. You can optionally pass uv which are the uv positions you will compute the NFFT at. This can allow for the NFFT plan to be cached improving performance

Fields

• padfac: Amount to pad the image

• reltol: relative tolerance of the NFFT

• precompute: NFFT interpolation algorithm.

• fftflags: Flag passed to inner AbstractFFT. The fastest FFTW is FFTW.MEASURE but takes the longest to precompute

InterpolatedModel(model, grid; algorithm)

Computes an representation of the model in the Fourier domain using interpolations and FFTs. This is useful to construct models that aren't directly representable in the Fourier domain.

Note

This is mostly used for testing and debugging purposes. In general people should use the FourierDualDomain functionality to compute the Fourier transform of a model.

modify(m::AbstractModel, transforms...)

Modify a given model using the set of transforms. This is the most general function that allows you to apply a sequence of model transformation for example

modify(Gaussian(), Stretch(2.0, 1.0), Rotate(π/4), Shift(1.0, 2.0), Renorm(2.0))

will create a asymmetric Gaussian with position angle π/4 shifted to the position (1.0, 2.0) with a flux of 2 Jy. This is similar to Flux's chain.

modify(img::IntensityMap, transforms...)

This modifies the img by applying the transforms... returning a transformed IntensityMap

Unlike when modify is applied to a <:AbstractModel this returns an already modified image.

basemodel(model::ModifiedModel)

Returns the ModifiedModel with the last transformation stripped.

Example

julia> basemodel(stretched(Disk(), 1.0, 2.0)) == Disk()
true

unmodified(model::ModifiedModel)

Returns the un-modified model

Example

julia> m = stretched(rotated(Gaussian(), π/4), 2.0, 1.0)
julia> umodified(m) == Gaussian()
true

renormed(model, f)

Renormalizes the model m to have total flux f*flux(m). This can also be done directly by calling Base.:* i.e.,

julia> renormed(m, f) == f*M
true

rotated(model, ξ)

Rotates the model by an amount ξ in radians in the clockwise direction.

posangle(model)

Returns the rotation angle of the rotated model

shifted(model, Δx, Δy)

Shifts the model m in the image domain by an amount Δx,Δy in the x and y directions respectively.

stretched(model, α, β)

Stretches the model m according to the formula Iₛ(x,y) = 1/(αβ) I(x/α, y/β), where were renormalize the intensity to preserve the models flux.

struct ModifiedModel{M, MT<:Tuple} <: ComradeBase.AbstractModel

Container type for models that have been transformed in some way. For a list of potential modifiers or transforms see subtypes(ModelModifiers).

Fields

• model: base model

• transform: model transforms

abstract type ModelModifier{T}

General type for a model modifier. These transform any model using simple Fourier transform properties. To modify a model you can use the ModifiedModel constructor or the modify function.

julia> visanalytic(stretched(Disk(), 2.0, 2.0)) == visanalytic(Disk())
true

To implement a model transform you need to specify the following methods:

• transform_uv
• transform_image
• scale_uv
• scale_image
• radialextent

See thee docstrings of those methods for guidance on implementation details.

Additionally these methods assume the modifiers are of the form

I(x,y) -> fᵢ(x,y)I(gᵢ(x,y)) V(u,v) -> fᵥ(u,v)V(gᵥ(u,v))

where g are the transformimage/uv functions and f are the scaleimage/uv function.

Stretch(α, β)
Stretch(r)

Stretched the model in the x and y directions, i.e. the new intensity is Iₛ(x,y) = 1/(αβ) I(x/α, y/β), where were renormalize the intensity to preserve the models flux.

Example

julia> modify(Gaussian(), Stretch(2.0)) == stretched(Gaussian(), 2.0, 1.0)
true

If only a single argument is given it assumes the same stretch is applied in both direction.

julia> Stretch(2.0) == Stretch(2.0, 2.0)
true

struct Shift{T} <: VLBISkyModels.ModelModifier{T}

Shifts the model by Δx units in the x-direction and Δy units in the y-direction.

Example

julia> modify(Gaussian(), Shift(2.0, 1.0)) == shifted(Gaussian(), 2.0, 1.0)
true

Rotate(ξ)

Type for the rotated model. This is more fine grained constrol of rotated model.

Example

julia> modify(Gaussian(), Rotate(2.0)) == rotated(Gaussian(), 2.0)
true

struct Renormalize{T} <: VLBISkyModels.ModelModifier{T}

Renormalizes the flux of the model to the new value scale*flux(model). We have also overloaded the Base.:* operator as syntactic sugar although I may get rid of this.

Example

julia> modify(Gaussian(), Renormalize(2.0)) == 2.0*Gaussian()
true

imageviz(img::IntensityMap; scale_length = fieldofview(img.X/4), kwargs...)

A default image visualization for a IntensityMap.

If eltype(img) <: Real i.e. an image of a single stokes parameter this will plot the image with a colorbar in units of Jy/μas². The plot will accept any kwargs that are a supported by the Makie.Image type an can be queried by typing ?image in the REPL

If eltype(img) <: StokesParams i.e. full polarized image this will use polimage. The plot will accept any kwargs that are a supported by the PolImage type an can be queried by typing ?polimage in the REPL.

polimage(img::IntensityMap{<:StokesParams};
            colormap = :bone,
            colorrange = Makie.automatic,
            pcolorrange=Makie.automatic,
            pcolormap=Reverse(:jet1),
            nvec = 30,
            min_frac = 0.1,
            min_pol_frac=0.2,
            length_norm=1.0,
            plot_total=true)

Plot a polarized intensity map using the image img.

The plot follows the conventions from EHTC M87 Paper VII.

The stokes I image will be plotted with the attributes

• colormap
• colorrange
• alpha
• colorscale

The polarized image will consist of a set. The attribute plot_total changes what polarized quantities are considered.

If plot_total = true

• The total polarization will be considered and the markers will be given by ellipses.
• The orientation of the ellipse is equal to the EVPA.
• The area of the ellipse is proportional to |V|².
• The semi-major axis is related to the total polarized intensity times the length_norm.
• The color of the ellipse is given by the fractional total polarization times the

sign of Stokes V.

If plot_total = false

• Only the linear polarization is considered and the markers will be ticks.
• The orientation of the ticks is equal to the EVPA.
• The length of the ticks is equal to the total linear polarized intensity, i.e. √(Q² + U²) times the length_norm.
• The color of the tick is given by the fractional linear polarization.

Attributes

• colormap: The colormap of the stokes I image. The default is :bone.
• colorrange: The color range of the stokes I image. The default is (0, maximum(stokes(img, :I)))
• pcolorrange: The color range for the polarized image
• pcolormap: The colormap used for fractional total/linear polarization markers.
• nvec: The number of polarization vectors to plot
• min_frac: Any markers with I < min_frac*maximum(I)) will not be plotted
• min_pol_frac: Any markers with P < min_frac*maximum(P)) where P is the total/linear polarization flux will not be plotted.
• length_norm: Specifies the normalization used for the ticks. The default is that the pixel with the largest polarization intensity will have a tick length = 10x the pixel separation. For an image with a maximum polarized intensity of 10Jy/μas² and pixel spacing of 1μas the marker length will be equal to 10μas.
• plot_total: If true plot the total polarization. If false only plot the linear polarization.

center_image(img::SpatialIntensityMap)

centers the img such that the centroid of the image is approximately at the origin.

convolve(img::IntensityMap, m::AbstractModel)

Convolves an img with a given analytic model m. This is useful for blurring the image with some model. For instance to convolve a image with a Gaussian you would do

convolve(img, Gaussian())

For the inplace version of the function see convolve!

Notes

This method does not automatically pad your image. If there is substantial flux at the boundaries you will start to see artifacts.

convolve!(img::IntensityMap, m::AbstractModel)

Convolves an img with a given analytic model m. This is useful for blurring the image with some model. For instance to convolve a image with a Gaussian you would do

convolve!(img, Gaussian())

Notes

This method does not automatically pad your image. If there is substantial flux at the boundaries you will start to see artifacts.

smooth(img::SpatialIntensityMap)

Smooths the img using a symmetric Gaussian with σ standard deviation.

For more flexible convolution please see convolve.

PoincareSphere2Map(I, p, X, grid)
PoincareSphere2Map(I::IntensityMap, p, X)
PoincareSphere2Map(I, p, X, cache::AbstractCache)

Constructs an polarized intensity map model using the Poincare parameterization. The arguments are:

• I is a grid of fluxes for each pixel.
• p is a grid of numbers between 0, 1 and the represent the total fractional polarization
• X is a grid, where each element is 3 numbers that represents the point on the Poincare sphere that is, X[1,1] is a NTuple{3} such that ||X[1,1]|| == 1.
• grid is the dimensional grid that gives the pixels locations of the intensity map.

PolExp2Map(a, b, c, d, grid::AbstractRectiGrid)

Constructs an polarized intensity map model using the matrix exponential representation from Arras 2021 (Thesis).

Each Stokes parameter is parameterized as

I = exp(a)cosh(p)
Q = exp(a)sinh(p)b/p
U = exp(a)sinh(p)c/p
V = exp(a)sinh(p)d/p

where a,b,c,d are real numbers with no conditions, and p=√(a² + b² + c²).

regrid(img, g)

Regrids the spatial parts of an image img on the new domain g

rad2μas(x)

Converts a number from radians to micro-arcseconds (μas)

μas2rad(x)

Converts a number from micro-arcseconds (μas) to rad

load_clean_components(fname::AbstractString, beam=nothing)

Load a clean component model from a file. The file can be a FITS file or a .mod file. If the beam argument is not given it will try to extract the beam from the FITS file.

load_fits(fitsfile::String, IntensityMap)

This loads in a fits file that is more robust to the various imaging algorithms in the EHT, i.e. is works with clean, smili, eht-imaging. The function returns an tuple with an intensitymap and a second named tuple with ancillary information about the image, like the source name, location, mjd, and radio frequency.

save_fits(file::String, img::IntensityMap, obs)

Saves an image to a fits file. You can optionally pass an EHTObservation so that ancillary information will be added.

abstract type AbstractFourierDualDomain <: ComradeBase.AbstractDomain

This defines an abstract cache that can be used to hold or precompute some computations.

scale_image(model::AbstractModifier, u, u)

Returns a number on how to scale the image visibility at u v for an modified model

scale_image(model::AbstractModifier, x, y)

Returns a number of how to to scale the image intensity at x y for an modified model

transform_uv(model::AbstractModifier, u, u)

Returns a transformed u and v according to the model modifier

transform_image(model::AbstractModifier, x, y)

Returns a transformed x and y according to the model modifier

unpack(m::AbstractModel)

Unpacks all elements of a struct into a named tuple. Note that this may include elements that aren't direcltly model parameters.

stokes_intensitymap(I, Q, U, V)

Constructs an IntensityMap from four maps for I, Q, U, V.

struct FFTPlan{A<:FFTAlg, P} <: VLBISkyModels.AbstractPlan

The cache used when the FFT algorithm is used to compute visibilties. This is an internal type and is not part of the public API

struct NUFTPlan{A, P, M, I, T} <: VLBISkyModels.AbstractNUFTPlan

Internal type used to store the cache for a non-uniform Fourier transform (NUFT).

The user should instead create this using the FourierDualDomain function.