The Frankenstein Python API

Geometry classes

Given a set of visibilities, together these classes: (1) optionally fit for the source geometry and (2) deproject the visibilities by the given or fitted geometry.

class frank.geometry.FixedGeometry(inc, PA, dRA=0.0, dDec=0.0)

Disc Geometry class using pre-determined parameters.

Centre and deproject the source to ensure axisymmetry

Parameters:

inc (float, unit = deg) – Disc inclination
PA (float, unit = deg) – Disc positition angle.
dRA (float, default = 0, unit = arcsec) – Phase centre offset in right ascension
dDec (float, default = 0, unit = arcsec) – Phase centre offset in declination

Notes

The phase centre offsets, dRA and dDec, refer to the distance to the source from the phase centre.

class frank.geometry.FitGeometryGaussian(inc_pa=None, phase_centre=None, guess=None)

Determine the disc geometry by fitting a Gaussian in Fourier space.

Centre and deproject the source to ensure axisymmetry

Parameters:

inc_pa (tuple = (inclination, position angle) or None (default), unit = deg) – Determine whether to fit for the source’s inclination and position angle. If inc_pa = None, the inclination and PA are fit for. Else inc_pa should be provided as a tuple
phase_centre (tuple = (dRA, dDec) or None (default), unit = arcsec) – Determine whether to fit for the source’s phase centre. If phase_centre = None, the phase centre is fit for. Else the phase centre should be provided as a tuple
guess (list of len(4), default = None) – Initial guess for the source’s inclination [deg], position angle [deg], right ascension offset [arcsec], declination offset [arcsec].

Notes

The phase centre offsets, dRA and dDec, refer to the distance to the source from the phase centre.

class frank.geometry.FitGeometryFourierBessel(Rmax, N, inc_pa=None, phase_centre=None, guess=None, verbose=False)

Determine the disc geometry by fitting a non-parametric brightness profile in visibility space.

The best fit is obtained by finding the geometry that minimizes the weighted chi^2 of the visibility fit.

The brightness profile is modelled using the FourierBesselFitter, which is equivalent to a FrankFitter fit without the Gaussian Process prior. For this reason, a small number of bins is recommended for fit stability.

Parameters:

Rmax (float, unit = arcsec) –

Radius of support for the functions to transform, i.e.,
f(r) = 0 for R >= Rmax
N (int) – Number of collocation points
inc_pa (tuple = (inclination, position angle) or None (default), unit = deg) – Determine whether to fit for the source’s inclination and position angle. If inc_pa = None, the inclination and PA are fit for. Else inc_pa should be provided as a tuple
phase_centre (tuple = (dRA, dDec) or None (default), unit = arcsec) – Determine whether to fit for the source’s phase centre. If phase_centre = None, the phase centre is fit for. Else the phase centre should be provided as a tuple
guess (list of len(4), default = None) – Initial guess for the source’s inclination [deg], position angle [deg], right ascension offset [arcsec], and declination offset [arcsec]
verbose (bool, default=False) – Determines whether to print the iteration progress.

class frank.geometry.SourceGeometry(inc=None, PA=None, dRA=None, dDec=None)

Base class for geometry corrections.

Centre and deproject the source to ensure axisymmetry

Parameters:

inc (float, unit = deg) – Inclination of the disc
PA (float, unit = deg) – Position angle of the disc
dRA (float, unit = arcsec) – Phase centre offset in right ascension.
dDec (float, units = arcsec) – Phase centre offset in declination.

Notes

The phase centre offsets, dRA and dDec, refer to the distance to the source from the phase centre.

apply_correction(u, v, V, use3D=False)

Correct the phase centre and deproject the visibilities

Parameters:

u (array of real, size = N, unit = \(\lambda\)) – u-points of the visibilities
v (array of real, size = N, unit = \(\lambda\)) – v-points of the visibilities
V (array of real, size = N, units = Jy) – Complex visibilites
use3D (bool, default=False) – If True, also return the 3rd compoent of the de-projected visibilities, wp.

Returns:

up (array of real, size = N, unit = \(\lambda\)) – Corrected u-points of the visibilities
vp (array of real, size = N, unit = \(\lambda\)) – Corrected v-points of the visibilities
wp (array of real, size = N, unit = \(\lambda\)) – [Optional] Corrected w-points of the visibilities
Vp (array of real, size = N, unit = Jy) – Corrected complex visibilites

deproject(u, v, use3D=False)

Convert uv-points from sky-plane to deprojected space (u,v)

Parameters:

u (array of real, size = N, unit = \(\lambda\)) – u-points of the visibilities
v (array of real, size = N, unit = \(\lambda\)) – v-points of the visibilities
use3D (bool, default=False) – If True, also return the 3rd compoent of the de-projected visibilities, wp.

Returns:

up (array of real, size = N, unit = \(\lambda\)) – Corrected u-points of the visibilities
vp (array of real, size = N, unit = \(\lambda\)) – Corrected v-points of the visibilities
wp (array of real, size = N, unit = \(\lambda\)) – [Optional] Corrected w-points of the visibilities

fit(u, v, V, weights)

Determine geometry using the provided uv-data

Parameters:

u (array of real, size = N, unit = \(\lambda\)) – u-points of the visibilities
v (array of real, size = N, unit = \(\lambda\)) – v-points of the visibilities
V (array of complex, size = N, unit = Jy) – Complex visibilites
weights (array of real, size = N, unit = Jy) – Weights on the visibilities

reproject(u, v)

Convert uv-points from deprojected space to sky-plane

Parameters:

u (array of real, size = N, unit = \(\lambda\)) – u-points of the visibilities
v (array of real, size = N, unit = \(\lambda\)) – v-points of the visibilities

Returns:

up (array of real, size = N, unit = \(\lambda\)) – Corrected u-points of the visibilities
vp (array of real, size = N, unit = \(\lambda\)) – Corrected v-points of the visibilities

undo_correction(u, v, V)

Undo the phase centre correction and deprojection

Parameters:

u (array of real, size = N, unit = \(\lambda\)) – u-points of the visibilities
v (array of real, size = N, unit = \(\lambda\)) – v-points of the visibilities
V (array of real, size = N, unit = Jy) – Complex visibilites

Returns:

up (array of real, size = N, unit = \(\lambda\)) – Corrected u-points of the visibilities
vp (array of real, size = N, unit = \(\lambda\)) – Corrected v-points of the visibilities
Vp (array of real, size = N, unit = Jy) – Corrected complex visibilites

Fitting classes

Together these classes reconstruct the 1D radial brightness profile of a source by fitting the deprojected visibilities.

class frank.radial_fitters.FrankFitter(Rmax, N, geometry, nu=0, block_data=True, block_size=100000, alpha=1.05, p_0=None, weights_smooth=0.0001, tol=0.001, method='Normal', I_scale=100000.0, max_iter=2000, check_qbounds=True, store_iteration_diagnostics=False, assume_optically_thick=True, scale_height=None, verbose=True, convergence_failure='raise')

Fit a Gaussian process model using the Discrete Hankel Transform of Baddour & Chouinard (2015).

The GP model is based upon Oppermann et al. (2013), which use a maximum a posteriori estimate for the power spectrum as the GP prior for the real-space coefficients

Parameters:

Rmax (float, unit = arcsec) – Radius of support for the functions to transform, i.e., f(r) = 0 for R >= Rmax.
N (int) – Number of collaction points
geometry (SourceGeometry object) – Geometry used to deproject the visibilities before fitting
nu (int, default = 0) – Order of the discrete Hankel transform, given by J_nu(r)
block_data (bool, default = True) – Large temporary matrices are needed to set up the data. If block_data is True, we avoid this, limiting the memory requirement to block_size elements
block_size (int, default = 10**5) – Size of the matrices if blocking is used
alpha (float >= 1, default = 1.05) – Order parameter of the inverse gamma prior for the power spectrum coefficients
p_0 (float >= 0, default = None, unit=Jy^2) – Scale parameter of the inverse gamma prior for the power spectrum coefficients. If not provided p_0 = 1e-15 (method=”Normal”) or 1e-35 (method=”LogNormal”) will be used.
weights_smooth (float >= 0, default = 1e-4) – Spectral smoothness prior parameter. Zero is no smoothness prior
tol (float > 0, default = 1e-3) – Tolerence for convergence of the power spectrum iteration
method (string, default="Normal") – Model used for the brightness reconstrution. This must be one of “Normal” of “LogNormal”.
I_scale (float, default = 1e5, unit= Jy/Sr) – Brightness scale. Only used in the LogNormal model. Note the LogNormal model produces I(Rmax) = I_scale.
max_iter (int, default = 2000) – Maximum number of fit iterations
check_qbounds (bool, default = True) – Whether to check if the first (last) collocation point is smaller (larger) than the shortest (longest) deprojected baseline in the dataset
store_iteration_diagnostics (bool, default = False) – Whether to store the power spectrum parameters and brightness profile for each fit iteration
assume_optically_thick (bool, default = True) – Whether to correct the visibility amplitudes by a factor of 1 / cos(inclination); see frank.geometry.rescale_total_flux
scale_height (function R --> H, optional) – Specifies the vertical thickness of disc as a function of radius. Both R and H should be in arcsec. Assumes a Gaussian vertical structure. Only works with assume_optically_thick=False
verbose (bool) – Whether to print notification messages
convergence_failure (string, default = 'raise') – Decide what to do when the frank model does not converge within max_iter. Should be one of: ‘raise’ : raise an error ‘warn’ : print a warning message and continue ‘ignore’ : Ignore the error.

References

Baddour & Chouinard (2015): DOI: https://doi.org/10.1364/JOSAA.32.000611
Oppermann et al. (2013): DOI: https://doi.org/10.1103/PhysRevE.87.032136

property MAP_solution: Reconstruction for the maximum a posteriori power spectrum

property MAP_spectrum: Maximum a posteriori power spectrum

property MAP_spectrum_covariance: Covariance matrix of the maximum a posteriori power spectrum

property Qmax: Maximum frequency, unit = \(\lambda\)

property Rmax: Maximum radius, unit = arcsec

fit(u, v, V, weights=1)

Fit the visibilties

Parameters:

u (1D array, unit = \(\lambda\)) – uv-points of the visibilies
v (1D array, unit = \(\lambda\)) – uv-points of the visibilies
V (1D array, unit = Jy) – Visibility amplitudes at q
weights (1D array, optional, unit = J^-2) – Weights of the visibilities, weight = 1 / sigma^2, where sigma is the standard deviation

Returns:

sol – Least-squares Fourier-Bessel series fit

Return type:

FrankRadialFit

property geometry: Geometry object

property q: Frequency points, unit = \(\lambda\)

property r: Radius points, unit = arcsec

property size: Number of points in reconstruction

class frank.radial_fitters.FrankRadialFit(vis_map, info, geometry)

Base class for results of frank fits.

Parameters:

vis_map (VisibilityMapping object) – Mapping between image and visibility plane.
info (dict) – Dictionary containing useful quantities for reproducing a fit (such as the hyperparameters used)
geometry (SourceGeometry object, optional) – Geometry used to correct the visibilities for the source inclination. If not provided, the geometry determined during the fit will be used.

property Qmax: Maximum frequency, unit = \(\lambda\)

property Rmax: Maximum radius, unit = arcsec

property geometry: SourceGeometry object

predict(u, v, I=None, geometry=None)

Predict the visibilities in the sky-plane

Parameters:

u (array, unit = \(\lambda\)) – uv-points to predict the visibilities at
v (array, unit = \(\lambda\)) – uv-points to predict the visibilities at
I (array, optional, unit = Jy) – Intensity points to predict the vibilities of. If not specified, the mean will be used. The intensity should be specified at the collocation points, I[k] = \(I(r_k)\)
geometry (SourceGeometry object, optional) – Geometry used to correct the visibilities for the source inclination. If not provided, the geometry determined during the fit will be used

Returns:

V(u,v) – Predicted visibilties of a source with a radial flux distribution given by \(I\) and the position angle, inclination and phase centre determined by the geometry object

Return type:

array, unit = Jy

predict_deprojected(q=None, I=None, geometry=None, block_size=100000, assume_optically_thick=True)

Predict the visibilities in the deprojected-plane

Parameters:

q (array, default = self.q, unit = \(\lambda\)) – 1D uv-points to predict the visibilities at
I (array, optional, unit = Jy / sr) – Intensity points to predict the vibilities of. If not specified, the mean will be used. The intensity should be specified at the collocation points, I[k] = I(r_k)
geometry (SourceGeometry object, optional) – Geometry used to correct the visibilities for the source inclination. If not provided, the geometry determined during the fit will be used
block_size (int, default = 10**5) – Maximum matrix size used in the visibility calculation
assume_optically_thick (bool, default = True) – Whether to correct the visibility amplitudes for the source inclination

Returns:

V(q) – Predicted visibilties of a source with a radial flux distribution given by \(I\). The amplitude of the visibilities are reduced according to the inclination of the source, for consistency with uvplot

Return type:

array, unit = Jy

Notes

The visibility amplitudes are still reduced due to the projection, for consistency with uvplot

property q: Frequency points, unit = \(\lambda\)

property r: Radius points, unit = arcsec

property size: Number of points in reconstruction

class frank.radial_fitters.FrankGaussianFit(DHT, fit, info={}, geometry=None)

Result of a frank fit with a Gaussian brightness model.

Parameters:

DHT (DiscreteHankelTransform) – A DHT object with N bins that defines H(p). The DHT is used to compute \(S(p)\)
fit (GaussianModel object) – Result of fitting with MAP power spectrum.
info (dict, optional) – Dictionary containing useful quantities for reproducing a fit (such as the hyperparameters used)
geometry (SourceGeometry object, optional) – Geometry used to correct the visibilities for the source inclination. If not provided, the geometry determined during the fit will be used.

property MAP: Posterior maximum, unit = Jy / sr

property Qmax: Maximum frequency, unit = \(\lambda\)

property Rmax: Maximum radius, unit = arcsec

property covariance: Posterior covariance, unit = (Jy / sr)**2

property geometry: SourceGeometry object

log_likelihood(I=None)

Compute one of two types of likelihood.

If \(I\) is provided, this computes

Otherwise the marginalized likelihood is computed,

Parameters:: I (array, size = N, optional, unit = Jy / sr) – Intensity \(I(r)\) to compute the likelihood of
Returns:: log_P – Log likelihood, \(\log[P(I,V|p)]\) or \(\log[P(V|p)]\)
Return type:: float

Notes

The prior probability P(S) is not included.
The likelihoods take the form:

\[\log[P(I,V|p)] = \frac{1}{2} j^T I - \frac{1}{2} I^T D^{-1} I - \frac{1}{2} \log[\det(2 \pi S)] + H_0\]

and

\[\log[P(V|p)] = \frac{1}{2} j^T D j + \frac{1}{2} \log[\det(D)/\det(S)] + H_0\]

where

\[H_0 = -\frac{1}{2} V^T w V + \frac{1}{2} \sum \log(w /2 \pi)\]

is the noise likelihood.

property mean: Posterior mean, unit = Jy / sr

property power_spectrum: Power spectrum coefficients

predict(u, v, I=None, geometry=None)

Predict the visibilities in the sky-plane

Parameters:

u (array, unit = \(\lambda\)) – uv-points to predict the visibilities at
v (array, unit = \(\lambda\)) – uv-points to predict the visibilities at
I (array, optional, unit = Jy) – Intensity points to predict the vibilities of. If not specified, the mean will be used. The intensity should be specified at the collocation points, I[k] = \(I(r_k)\)
geometry (SourceGeometry object, optional) – Geometry used to correct the visibilities for the source inclination. If not provided, the geometry determined during the fit will be used

Returns:

V(u,v) – Predicted visibilties of a source with a radial flux distribution given by \(I\) and the position angle, inclination and phase centre determined by the geometry object

Return type:

array, unit = Jy

predict_deprojected(q=None, I=None, geometry=None, block_size=100000, assume_optically_thick=True)

Predict the visibilities in the deprojected-plane

Parameters:

q (array, default = self.q, unit = \(\lambda\)) – 1D uv-points to predict the visibilities at
I (array, optional, unit = Jy / sr) – Intensity points to predict the vibilities of. If not specified, the mean will be used. The intensity should be specified at the collocation points, I[k] = I(r_k)
geometry (SourceGeometry object, optional) – Geometry used to correct the visibilities for the source inclination. If not provided, the geometry determined during the fit will be used
block_size (int, default = 10**5) – Maximum matrix size used in the visibility calculation
assume_optically_thick (bool, default = True) – Whether to correct the visibility amplitudes for the source inclination

Returns:

V(q) – Predicted visibilties of a source with a radial flux distribution given by \(I\). The amplitude of the visibilities are reduced according to the inclination of the source, for consistency with uvplot

Return type:

array, unit = Jy

Notes

The visibility amplitudes are still reduced due to the projection, for consistency with uvplot

property q: Frequency points, unit = \(\lambda\)

property r: Radius points, unit = arcsec

property size: Number of points in reconstruction

solve_non_negative(): Compute the best fit solution with non-negative intensities

class frank.radial_fitters.FrankLogNormalFit(DHT, fit, info={}, geometry=None)

Result of a frank fit with a Gaussian brightness model.

Parameters:

DHT (DiscreteHankelTransform) – A DHT object with N bins that defines H(p). The DHT is used to compute \(S(p)\)
fit (LogNormalMAPModel object) – Result of fitting with MAP power spectrum.
info (dict, optional) – Dictionary containing useful quantities for reproducing a fit (such as the hyperparameters used)
geometry (SourceGeometry object, optional) – Geometry used to correct the visibilities for the source inclination. If not provided, the geometry determined during the fit will be used.

property MAP: Posterior maximum, unit = Jy / sr

property Qmax: Maximum frequency, unit = \(\lambda\)

property Rmax: Maximum radius, unit = arcsec

property covariance: Posterior covariance, unit = log[(Jy / sr)**2]

property geometry: SourceGeometry object

log_likelihood(I=None)

Compute the likelihood,

Parameters:: I (array, size = N, optional) – Intensity \(I(r)=exp(s0*s)\) to compute the likelihood of
Returns:: log_P – Log likelihood, \(\log[P(I,V|p)]\)
Return type:: float

Notes

The prior probability P(S) is not included.
The likelihood takes the form:

\[\log[P(I,V|p)] = j^T I - \frac{1}{2} I^T D^{-1} I - \frac{1}{2} \log[\det(2 \pi S)] + H_0\]

where

\[H_0 = -\frac{1}{2} V^T w V + \frac{1}{2} \sum \log(w /2 \pi)\]

is the noise likelihood.

property power_spectrum: Power spectrum coefficients

predict(u, v, I=None, geometry=None)

Predict the visibilities in the sky-plane

Parameters:

u (array, unit = \(\lambda\)) – uv-points to predict the visibilities at
v (array, unit = \(\lambda\)) – uv-points to predict the visibilities at
I (array, optional, unit = Jy) – Intensity points to predict the vibilities of. If not specified, the mean will be used. The intensity should be specified at the collocation points, I[k] = \(I(r_k)\)
geometry (SourceGeometry object, optional) – Geometry used to correct the visibilities for the source inclination. If not provided, the geometry determined during the fit will be used

Returns:

V(u,v) – Predicted visibilties of a source with a radial flux distribution given by \(I\) and the position angle, inclination and phase centre determined by the geometry object

Return type:

array, unit = Jy

predict_deprojected(q=None, I=None, geometry=None, block_size=100000, assume_optically_thick=True)

Predict the visibilities in the deprojected-plane

Parameters:

q (array, default = self.q, unit = \(\lambda\)) – 1D uv-points to predict the visibilities at
I (array, optional, unit = Jy / sr) – Intensity points to predict the vibilities of. If not specified, the mean will be used. The intensity should be specified at the collocation points, I[k] = I(r_k)
geometry (SourceGeometry object, optional) – Geometry used to correct the visibilities for the source inclination. If not provided, the geometry determined during the fit will be used
block_size (int, default = 10**5) – Maximum matrix size used in the visibility calculation
assume_optically_thick (bool, default = True) – Whether to correct the visibility amplitudes for the source inclination

Returns:

V(q) – Predicted visibilties of a source with a radial flux distribution given by \(I\). The amplitude of the visibilities are reduced according to the inclination of the source, for consistency with uvplot

Return type:

array, unit = Jy

Notes

The visibility amplitudes are still reduced due to the projection, for consistency with uvplot

property q: Frequency points, unit = \(\lambda\)

property r: Radius points, unit = arcsec

property size: Number of points in reconstruction

class frank.debris_fitters.FrankDebrisFitter(Rmax, N, geometry, scale_height, nu=0, block_data=True, block_size=100000, alpha=1.05, p_0=None, weights_smooth=0.0001, tol=0.001, method='Normal', I_scale=100000.0, max_iter=2000, check_qbounds=True, store_iteration_diagnostics=False, verbose=True, convergence_failure='raise')

Fit a Gaussian process model using the Discrete Hankel Transform of Baddour & Chouinard (2015).

The brightness model is \(I(R, z) = I(R) exp(-z^2/2H(R)^2)\), where \(H(R)\) is the (known) scale-height.

The GP model is based upon Oppermann et al. (2013), which use a maximum a posteriori estimate for the power spectrum as the GP prior for the real-space coefficients.

Parameters:

Rmax (float, unit = arcsec) – Radius of support for the functions to transform, i.e., f(r) = 0 for R >= Rmax.
N (int) – Number of collaction points
geometry (SourceGeometry object) – Geometry used to deproject the visibilities before fitting
scale_height (function R --> H) – Specifies the thickness of disc as a function of radius. Both units should be in arcsec.
nu (int, default = 0) – Order of the discrete Hankel transform, given by J_nu(r)
block_data (bool, default = True) – Large temporary matrices are needed to set up the data. If block_data is True, we avoid this, limiting the memory requirement to block_size elements
block_size (int, default = 10**5) – Size of the matrices if blocking is used
alpha (float >= 1, default = 1.05) – Order parameter of the inverse gamma prior for the power spectrum coefficients
p_0 (float >= 0, default = None, unit=Jy^2) – Scale parameter of the inverse gamma prior for the power spectrum coefficients. If not provided p_0 = 1e-15 (method=”Normal”) or 1e-35 (method=”LogNormal”) will be used.
weights_smooth (float >= 0, default = 1e-4) – Spectral smoothness prior parameter. Zero is no smoothness prior
tol (float > 0, default = 1e-3) – Tolerence for convergence of the power spectrum iteration
method (string, default="Normal") – Model used for the brightness reconstrution. This must be one of “Normal” of “LogNormal”.
I_scale (float, default = 1e5, unit= Jy/Sr) – Brightness scale. Only used in the LogNormal model. Notet the LogNormal model produces I(Rmax) = I_scale.
max_iter (int, default = 2000) – Maximum number of fit iterations
check_qbounds (bool, default = True) – Whether to check if the first (last) collocation point is smaller (larger) than the shortest (longest) deprojected baseline in the dataset
store_iteration_diagnostics (bool, default = False) – Whether to store the power spectrum parameters and brightness profile for each fit iteration
verbose – Whether to print notification messages
convergence_failure (string, default = 'raise') – Decide what to do when the frank model does not converge within max_iter. Should be one of: ‘raise’ : raise an error ‘warn’ : print a warning message and continue ‘ignore’ : Ignore the error.

References

Baddour & Chouinard (2015): DOI: https://doi.org/10.1364/JOSAA.32.000611
Oppermann et al. (2013): DOI: https://doi.org/10.1103/PhysRevE.87.032136

property MAP_solution: Reconstruction for the maximum a posteriori power spectrum

property MAP_spectrum: Maximum a posteriori power spectrum

property MAP_spectrum_covariance: Covariance matrix of the maximum a posteriori power spectrum

property Qmax: Maximum frequency, unit = \(\lambda\)

property Rmax: Maximum radius, unit = arcsec

fit(u, v, V, weights=1)

Fit the visibilties

Parameters:

u (1D array, unit = \(\lambda\)) – uv-points of the visibilies
v (1D array, unit = \(\lambda\)) – uv-points of the visibilies
V (1D array, unit = Jy) – Visibility amplitudes at q
weights (1D array, optional, unit = J^-2) – Weights of the visibilities, weight = 1 / sigma^2, where sigma is the standard deviation

Returns:

sol – Least-squares Fourier-Bessel series fit

Return type:

FrankRadialFit

property geometry: Geometry object

property q: Frequency points, unit = \(\lambda\)

property r: Radius points, unit = arcsec

property size: Number of points in reconstruction

Utility functions and classes

These are some useful functions and classes for various aspects of fitting and analysis.

Hankel transform

class frank.hankel.DiscreteHankelTransform(Rmax, N, nu=0)

Utilities for computing the discrete Hankel transform.

This class provides the necessary interface to compute a discrete version of the Hankel transform (DHT):

H[f](q) = int_0^R_{max} f(r) J_nu(2*pi*q*r) * 2*pi*r dr.

The DHT is based on [1].

Additionally this class provides coefficients of the DHT [1] transform matrix

Parameters:

Rmax (float) – Maximum radius beyond which f(r) is zero
N (integer) – Number of terms to use in the series
nu (integer, default = 0) – Order of the Bessel function, J_nu(r)

References

[1] Baddour & Chouinard (2015): DOI: https://doi.org/10.1364/JOSAA.32.000611 Note: the definition of the DHT used here differs by factors of 2*pi.

property Qmax: Maximum frequency

property Rmax: Maximum radius

coefficients(q=None, direction='forward')

Coefficients of the transform matrix, defined by: H[f](q) = np.dot(Y, f)

Parameters:

q (array or None) – Frequency points at which to evaluate the transform. If q = None, the points of the DHT are used. If direction=’backward’, these points should instead be the radius points
direction ({ 'forward', 'backward' }, optional) – Direction of the transform. If not supplied, the forward transform is used

Returns:

Y – The transformation matrix

Return type:

array, size = (len(q), N)

property order: Order of the Bessel function

property q: Frequency points

property r: Radius points

property size: Number of points used in the DHT

transform(f, q=None, direction='forward')

Compute the Hankel transform of an array

Parameters:

f (array, size = N) –

Function to Hankel transform, evaluated at the collocation points:
f[k] = f(r_k) or f[k] = f(q_k)
q (array or None) – The frequency points at which to evaluate the Hankel transform. If not specified, the conjugate points of the DHT will be used. For the backwards transform, q should be the radius points
direction ({ 'forward', 'backward' }, optional) – Direction of the transform. If not supplied, the forward transform is used

Returns:

H[f] – The Hankel transform of the array f

Return type:

array, size = N or len(q) if supplied

frank.utilities.generic_dht(x, f, Rmax=2.0, N=500, direction='forward', grid=None, inc=0.0)

Compute the visibilities or brightness of a model by directly applying the Discrete Hankel Transform.

The correction for inclination will also be applied, assuming an optically thick disc. For an optically thin disc, setting inc=0 (the default) will achieve the correct scaling.

Parameters:

x (array, unit = [arcsec] or [lambda]) – Radial or spatial frequency coordinates of f(x)
f (array, unit = [Jy / sr] or [Jy]) – Amplitude values of f(x)
Rmax (float, unit = [arcsec], default=2.0) – Maximum radius beyond which the real space function is zero
N (integer, default=500) – Number of terms to use in the Fourier-Bessel series
direction ({ 'forward', 'backward' }, default='forward') – Direction of the transform. ‘forward’ is real space -> Fourier space.
grid (array, unit = [arcsec] or [lambda], default=None) – The radial or spatial frequency points at which to sample the DHT. If None, the DHT collocation points will be used.
inc (float, unit = [deg], default = 0.0) – Source inclination. The total flux of the transform of f(x) will be scaled by cos(inc); this has no effect if inc=0.

Returns:

grid (array, size=N, unit = [arcsec] or [lambda]) – Spatial frequency or radial coordinates of the Hankel transform of f(x)
f_transform (array, size=N, unit = [Jy / sr] or [Jy]) – Hankel transform of f(x)

Notes

‘x’ and ‘f’ should be sufficiently sampled to ensure an interpolation will be accurate, otherwise ‘f_transform’ may be a poor estimate of their transform.

frank.utilities.get_collocation_points(Rmax=2.0, N=500, direction='forward')

Obtain the collocation points of a discrete Hankel transform for a given ‘Rmax’ and ‘N’ (see frank.hankel.DiscreteHankelTransform)

Parameters:

Rmax (float, unit = [arcsec], default=2.0) – Maximum radius beyond which the real space function is zero
N (integer, default=500) – Number of terms to use in the Fourier-Bessel series
direction ({ 'forward', 'backward' }, default='forward') – Direction of the transform. ‘forward’ is real space -> Fourier space, returning real space radial collocation points needed for the transform.

Returns:

coll_pts – The DHT collocation points in either real or Fourier space.

Return type:

array, unit = [lambda] or [arcsec]

Unit conversion

frank.utilities.arcsec_baseline(x)

Provide x as a radial scale [arcsec] to return the corresponding baseline [lambda], or vice-versa

Parameters:: x (float) – Radial scale [arcsec] or baseline [lambda]
Returns:: converted – Baseline [lambda] or radial scale [arcsec]
Return type:: float

frank.utilities.radius_convert(x, dist, conversion='arcsec_au')

Provide x as a radius/radii in [arcsec] to convert to [au] (or vice-versa), assuming a distance in [pc]

Parameters:

x (float or array) – Radius or radii ([arcsec] or [au])
dist (float) – Distance to source [pc]
conversion ({'arcsec_au', 'au_arcsec'}, default = 'arcsec_au') – The unit conversion to perform, e.g. ‘arcsec_au’ converts from [arcsec] to [au]

Returns:

converted – Radius or radii ([au] or [arcsec])

Return type:

float or array

frank.utilities.jy_convert(x, conversion, bmaj=None, bmin=None)

Provide x as a brightness in one of the units [Jy / beam], [Jy / arcsec^2], [Jy / sterad] to convert x to another of these units

Parameters:

x (array) – Brightness in one of: [Jy / beam], [Jy / arcsec^2], [Jy / sterad]
conversion ({ 'beam_sterad', 'beam_arcsec2', 'arcsec2_beam',) – ‘arcsec2_sterad’, ‘sterad_beam’, ‘sterad_arcsec2’} The unit conversion to perform, e.g., ‘beam_sterad’ converts [Jy / beam] to [Jy / sterad]
bmaj (float, optional) – Beam FWHM along the major axis [arcsec]
bmin (float, optional) – Beam FWHM along the minor axis [arcsec]

Returns:

converted – Brightness in converted units

Return type:

float

Data alteration

frank.utilities.normalize_uv(u, v, wle)

Normalize data u and v coordinates by the observing wavelength

Parameters:

u (array, unit = [m]) – u and v coordinates of observations
v (array, unit = [m]) – u and v coordinates of observations
wle (float or array, unit = [m]) – Observing wavelength of observations. If an array, it should be the pointwise wavelength for each (u,v) point

Returns:

u_normed, v_normed – u and v coordinates normalized by observing wavelength

Return type:

array, unit = \(\lambda\)

frank.utilities.cut_data_by_baseline(u, v, vis, weights, cut_range, geometry=None)

Truncate the data to be within a chosen baseline range.

The cut will be done in deprojected baseline space if the geometry is provided.

Parameters:

u (array, unit = \(\lambda\)) – u and v coordinates of observations
v (array, unit = \(\lambda\)) – u and v coordinates of observations
vis (array, unit = Jy) – Observed visibilities (complex: real + imag * 1j)
weights (array, unit = Jy^-2) – Weights assigned to observed visibilities, of the form \(1 / \sigma^2\)
cut_range (list of float, length = 2, unit = [lambda]) – Lower and upper baseline bounds outside of which visibilities are truncated
geometry (SourceGeometry object, optional) – Fitted geometry (see frank.geometry.SourceGeometry).

Returns:

u_cut, v_cut (array, unit = \(\lambda\)) – u and v coordinates in the chosen baseline range
vis_cut (array, unit = Jy) – Visibilities in the chosen baseline range
weights_cut (array, unit = Jy^-2) – Weights in the chosen baseline range

Visibility binning and weights estimation

class frank.utilities.UVDataBinner(uv, V, weights, bin_width)

Average uv-data into bins of equal size.

Compute the weighted mean of the visibilities in each bin

Parameters:

uv (array, unit = \(\lambda\)) – Baselines of the data to bin
V (array, unit = Jy) – Observed visibility. If complex, both the real and imaginary components will be binned. Else only the real part will be binned.
weights (array, unit = Jy^-2) – Weights on the visibility points
bin_width (float, unit = \(\lambda\)) – Width of the uv-bins

Notes

Uses numpy masked arrays to mask bins with no uv points.

frank.utilities.estimate_weights(u, v=None, V=None, nbins=300, log=True, use_median=False, verbose=True)

Estimate the weights using the variance of the binned visibilities.

The estimation is done assuming that the variation in each bin is dominated by the noise. This will be true if: 1) The source is axi-symmetric, 2) The uv-points have been deprojected, 3) The bins are not too wide, Otherwise the variance may be dominated by intrinsic variations in the visibilities.

Parameters:

u (array, unit = \(\lambda\)) – u and v coordinates of observations (deprojected). Data will be binned by baseline. If v is not None, np.hpot(u,v) will be used instead. Note that if V is None the argument v will be intepreted as V instead
v (array, unit = \(\lambda\)) – u and v coordinates of observations (deprojected). Data will be binned by baseline. If v is not None, np.hpot(u,v) will be used instead. Note that if V is None the argument v will be intepreted as V instead
V (array, unit = Jy, default = None) – Observed visibility. If complex, the weights will be computed from the average of the variance of the real and imaginary components, as in CASA’s statwt. Otherwise the variance of the real part is used.
nbins (int, default = 300) – Number of bins used.
log (bool, default = True) – If True, the uv bins will be constructed in log space, otherwise linear spaced bins will be used.
use_median (bool, default = False) – If True all of the weights will be set to the median of the variance estimated across the bins. Otherwise, the baseline dependent variance will be used.
verbose (bool, default = True) – If true, the logger will record calls to this function, along with whether the median estimate was used.

Returns:

weights – Estimate of the weight for each uv point.

Return type:

array

Notes

This function does not use the original weights in the estimation.
Bins with only one uv point do not have a variance estimate. Thus the mean of the variance in the two adjacent bins is used instead.

Examples

All of the following calls will work as expected:: estimate_weights(u, v, V) ` `estimate_weights(u, V) estimate_weights(u, V=V)

In each case the variance of V in the uv-bins is used to estimate the weights. The first call will use q = np.hypot(u, v) in the uv-bins. The second and third calls are equivalent to the first with u=0.

Imaging

frank.utilities.sweep_profile(r, I, project=False, phase_shift=False, geom=None, axis=0, xmax=None, ymax=None, dr=None)

Sweep a 1D radial brightness profile over \(2 \pi\) to yield a 2D brightness distribution. Optionally project this sweep by a supplied geometry.

Parameters:

r (array) – Radial coordinates at which the 1D brightness profile is defined
I (array) – Brightness values at r
project (bool, default = False) – Whether to project the swept profile by the supplied geom
phase_shift (bool, default = False) – Whether to phase shift the projected profile by the supplied geom. If False, the source will be centered in the image
geom (SourceGeometry object, default=None) – Fitted geometry (see frank.geometry.SourceGeometry). Here we use geom.inc [deg], geom.PA [deg], geom.dRA [arcsec], geom.dDec [arcsec] if project=True
axis (int, default = 0) – Axis over which to interpolate the 1D profile
xmax (float, optional, default = None) – Value setting the x- and y-bounds of the image (same units as r). The positive and negative bounds are both set to this value (modulo sign). If not provided, these will be set to r.max()
ymax (float, optional, default = None) – Value setting the x- and y-bounds of the image (same units as r). The positive and negative bounds are both set to this value (modulo sign). If not provided, these will be set to r.max()
dr (float, optional, default = None) – Pixel size (same units as r). If not provided, it will be set at the same spatial scale as r

Returns:

I2D (array, shape = (len(r), len(r))) – 2D brightness distribution (projected if project=True)
xmax (float) – Maximum x-value of the 2D grid
ymax (float) – Maximum y-value of the 2D grid

Notes

Sign convention: a negative geom.dRA shifts the source to the right in the image

frank.utilities.make_image(fit, Npix, xmax=None, ymax=None, project=True)

Make an image of a model fit.

Parameters:

fit (FrankFitter result object) – Fitted profile to make an image of
Npix (int or list) – Number of pixels in the x-direction, or [x-,y-] direction
xmax (float or None, unit=arcsec) – Size of the image is [-xmax, xmax]. By default this is twice fit.Rmax to avoid aliasing.
ymax (float or None, unit=arcsec) – Size of the image is [-ymax,ymax]. Defaults to xmax if ymax=None
project (bool, default=True) – Whether to produce a projected image.

Returns:

x (array, 1D; unit=arcsec) – Locations of the x-points in the image.
y (array, 1D; unit=arcsec) – Locations of the y-points in the image.
I (array, 2D; unit=Jy/Sr) – Image of the surface brightness.

frank.utilities.convolve_profile(r, I, disc_i, disc_pa, clean_beam, n_per_sigma=5, axis=0)

Convolve a 1D radial brightness profile with a 2D Gaussian beam, degrading the profile’s resolution

Parameters:

r (array) – Radial coordinates at which the 1D brightness profile is defined
I (array) – Brightness values at r
disc_i (float, unit = deg) – Disc inclination
disc_pa (float, unit = deg) – Disc position angle
clean_beam (dict) – Dictionary with beam bmaj (FWHM of beam along its major axis) [arcsec], bmin (FWHM of beam along its minor axis) [arcsec], pa (beam position angle) [deg]
n_per_sigma (int, default = 50) – Number of points per standard deviation of the Gaussian kernel (used for gridding)
axis (int, default = 0) – Axis over which to interpolate the 1D profile

Returns:

I_smooth – Convolved brightness profile I at coordinates r

Return type:

array, shape = (len(r), len(r))

Mock data routines

frank.utilities.add_vis_noise(vis, weights, seed=None)

Add Gaussian noise to visibilities

Parameters:

vis (array, unit = [Jy]) – Visibilities to add noise to. Can be complex (real + imag * 1j) or purely real.
weights (array, unit = [Jy^-2]) – Weights on the visibilities, of the form \(1 / \sigma^2\). Injected noise will be scaled proportional to sigma.
seed (int, default = None) – Number to initialize a pseudorandom number generator for the noise draws

Returns:

vis_noisy – Visibilities with added noise

Return type:

array, shape = vis.shape

frank.utilities.make_mock_data(r, I, Rmax, u, v, projection=None, geometry=None, N=500, add_noise=False, weights=None, seed=None)

Generate mock visibilities from a provided brightness profile and (u,v) distribution.

Parameters:

r (array, unit = [arcsec]) – Radial coordinates of I(r)
I (array, unit = [Jy / sr]) – Brightness values at r
Rmax (float, unit = [arcsec], default=2.0) – Maximum radius beyond which I(r) is zero. This should be larger than the disk size
u (array, unit = \(\lambda\)) – u and v coordinates of observations
v (array, unit = \(\lambda\)) – u and v coordinates of observations
projection (str, default = None) – One of [None, ‘deproject’, ‘reproject’] If None, the visibilities will be neither deprojected nor reprojected. If ‘deproject’ or ‘reproject’, the visibilites will be accordingly deprojected or reprojected by the supplied geometry and their total flux scaled by the inclination.
geometry (SourceGeometry object, default=None) – Source geometry (see frank.geometry.SourceGeometry). Must be supplied if projection is ‘deproject’ or ‘reproject’.
N (integer, default=500) – Number of terms to use in the Fourier-Bessel series
add_noise (bool, default = False) – Whether to add noise to the mock visibilities
weights (array, unit = Jy^-2) – Visibility weights, of the form \(1 / \sigma^2\). If provided, injected noise will be scaled proportional to sigma.
seed (int, default = None) – Number to initialize a pseudorandom number generator for the noise draws

Returns:

baselines (array, unit = \(\lambda\)) – Baseline coordinates of the mock visibilities. These will be equal to np.hypot(u, v) if ‘geometry’ is None (or if its keys are all equal to 0)
vis (array, unit = Jy) – Mock visibility amplitudes, including noise if ‘add_noise’ is True

Notes

‘r’ and ‘I’ should be sufficiently sampled to ensure an interpolation will be accurate, otherwise ‘vis’ may be a poor estimate of their transform.

Statistical analysis

frank.utilities.draw_bootstrap_sample(u, v, vis, weights)

Obtain the sample for a bootstrap, drawing, with replacement, N samples from a length N dataset

Parameters:

u (array, unit = \(\lambda\)) – u and v coordinates of observations
v (array, unit = \(\lambda\)) – u and v coordinates of observations
vis (array, unit = Jy) – Observed visibilities (complex: real + imag * 1j)
weights (array, unit = Jy^-2) – Weights on the visibilities, of the form \(1 / \sigma^2\)

Returns:

u_boot, v_boot (array, unit = \(\lambda\)) – Bootstrap sampled u and v coordinates
vis_boot (array, unit = Jy) – Bootstrap sampled visibilities
weights_boot (array, unit = Jy^-2) – Boostrap sampled weights on the visibilities