Libs

Libs contains functions that support the processing components but are not expected to interface to workflows.

Image

Iterators

Functions that define and manipulate images. Images are just data and a World Coordinate System.

image_channel_iter(im: data_models.memory_data_models.Image, subimages=1) → collections.abc.Iterable

Create a image_channel_iter generator, returning images

The WCS is adjusted appropriately for each raster element. Hence this is a coordinate-aware way to iterate through an image.

Provided we don’t break reference semantics, memory should be conserved

To update the image in place:

for r in raster(im, facets=2)::

r.data[…] = numpy.sqrt(r.data[…])

Parameters:

• im – Image
• subimages – Number of subimages

image_null_iter(im: data_models.memory_data_models.Image, facets=1, overlap=0) → collections.abc.Iterable

One time iterator

Parameters:

• im –
• facets – Number of image partitions on each axis (2)
• overlap – overlap in pixels

Returns:

image_raster_iter(im: data_models.memory_data_models.Image, facets=1, overlap=0, taper='flat', make_flat=False) → collections.abc.Iterable

Create an image_raster_iter generator, returning images, optionally with overlaps

The WCS is adjusted appropriately for each raster element. Hence this is a coordinate-aware way to iterate through an image.

Provided we don’t break reference semantics, memory should be conserved. However make_flat creates a new set of images and thus reference semantics dont hold.

To update the image in place:

for r in raster(im, facets=2)::

r.data[…] = numpy.sqrt(r.data[…])

If the overlap is greater than zero, we choose to keep all images the same size so the other ring of facets are ignored. So if facets=4 and overlap > 0 then the iterator returns (facets-2)**2 = 4 images.

A taper is applied in the overlap regions. None implies a constant value, linear is a ramp, and quadratic is parabolic at the ends.

Parameters:

• im – Image
• facets – Number of image partitions on each axis (2)
• overlap – overlap in pixels
• taper – method of tapering at the edges: ‘flat’ or ‘linear’ or ‘quadratic’ or ‘tukey’
• make_flat – Make the flat images

Cleaners

Image Deconvolution functions

argmax(a)

Return unravelled index of the maximum

param: a: array to be searched

calculate_scale_inverse_moment_moment_hessian(scale_scale_moment_moment_psf)

Calculate inverse_scale dependent moment moment hessian

Part of the initialisation for Algorithm 1. Lines 7 - 9

Parameters:scale_scale_moment_moment_psf – scale_moment_psf [nscales, nscales, nmoments, nmoments] Returns:scale-dependent moment-moment inverse hessian

calculate_scale_moment_principal_solution(smresidual, ihsmmpsf)

Calculate the principal solution in moment space for each scale

Lines 20 - 26

Parameters:

• smresidual – scale-dependent moment residual [nscales, nmoments, nx, ny]
• ihsmmpsf – Inverse of scale dependent moment moment Hessian

Returns:

Decoupled residual images [nscales, nmoments, nx, ny]

calculate_scale_moment_residual(residual, scalestack)

Calculate scale-dependent moment residuals

Part of the initialisation for Algorithm 1: lines 12 - 17

Parameters:

• scalestack –
• residual – residual [nmoments, nx, ny]

Returns:

scale-dependent moment residual [nscales, nmoments, nx, ny]

calculate_scale_scale_moment_moment_psf(psf, scalestack)

Calculate scale-dependent moment psfs

Part of the initialisation for Algorithm 1

Parameters:

• scalestack –
• psf – psf

Returns:

scale-dependent moment psf [nscales, nscales, nmoments, nmoments, nx, ny]

convolve_convolve_scalestack(scalestack, img)

Convolve img by the specified scalestack, returning the resulting stack

Parameters:

• scalestack – stack containing the scales
• img – Image to be convolved

Returns:

Twice convolved image [nscales, nscales, nx, ny]

convolve_scalestack(scalestack, img)

Convolve img by the specified scalestack, returning the resulting stack

Parameters:

• scalestack – stack containing the scales
• img – Image to be convolved

Returns:

stack

create_scalestack(scaleshape, scales, norm=True)

Create a cube consisting of the scales

Parameters:

• scaleshape – desired shape of stack
• scales – scales (in pixels)
• norm – Normalise each plane to unity?

Returns:

stack

find_global_optimum(hsmmpsf, ihsmmpsf, smresidual, windowstack, findpeak)

Find the optimum peak using one of a number of algorithms

find_max_abs_stack(stack, windowstack, couplingmatrix)

Find the location and value of the absolute maximum in this stack :param stack: stack to be searched :param windowstack: Window for the search :param couplingmatrix: Coupling matrix between difference scales :return: x, y, scale

find_optimum_scale_zero_moment(smpsol, windowstack)

Find the optimum scale for moment zero

Line 27 of Algorithm 1

Parameters:

• windowstack –
• smpsol – Decoupled residual images for each scale and moment

Returns:

x, y, optimum scale for peak

hogbom(dirty, psf, window, gain, thresh, niter, fracthresh, prefix='')

Clean the point spread function from a dirty image

See Hogbom CLEAN (1974A&AS…15..417H)

This version operates on numpy arrays.

Parameters:

• fracthresh –
• prefix –
• dirty – The dirty Image, i.e., the Image to be deconvolved
• psf – The point spread-function
• window – Regions where clean components are allowed. If True, entire dirty Image is allowed
• gain – The “loop gain”, i.e., the fraction of the brightest pixel that is removed in each iteration
• thresh – Cleaning stops when the maximum of the absolute deviation of the residual is less than this value
• niter – Maximum number of components to make if the threshold thresh is not hit

Returns:

clean component Image, residual Image

msclean(dirty, psf, window, gain, thresh, niter, scales, fracthresh, prefix='')

Perform multiscale clean

Multiscale CLEAN (IEEE Journal of Selected Topics in Sig Proc, 2008 vol. 2 pp. 793-801)

This version operates on numpy arrays.

Parameters:

• prefix –
• fracthresh –
• dirty – The dirty image, i.e., the image to be deconvolved
• psf – The point spread-function
• window – Regions where clean components are allowed. If True, all of the dirty image is allowed
• gain – The “loop gain”, i.e., the fraction of the brightest pixel that is removed in each iteration
• thresh – Cleaning stops when the maximum of the absolute deviation of the residual is less than this value
• niter – Maximum number of components to make if the threshold “thresh” is not hit
• scales – Scales (in pixels width) to be used

Returns:

clean component image, residual image

msmfsclean(dirty, psf, window, gain, thresh, niter, scales, fracthresh, findpeak='ARL', prefix='')

Perform image plane multiscale multi frequency clean

This algorithm is documented as Algorithm 1 in: U. Rau and T. J. Cornwell, “A multi-scale multi-frequency deconvolution algorithm for synthesis imaging in radio interferometry,” A&A 532, A71 (2011). Note that this is only the image plane parts.

Specific code is linked to specific lines in that algorithm description.

This version operates on numpy arrays that have been converted to moments on the last axis.

Parameters:

• fracthresh –
• dirty – The dirty image, i.e., the image to be deconvolved
• psf – The point spread-function
• window – Regions where clean components are allowed. If True, all of the dirty image is allowed
• gain – The “loop gain”, i.e., the fraction of the brightest pixel that is removed in each iteration
• thresh – Cleaning stops when the maximum of the absolute deviation of the residual is less than this value
• niter – Maximum number of components to make if the threshold “thresh” is not hit
• scales – Scales (in pixels width) to be used
• fracthresh – Fractional stopping threshold
• findpeak – Method of finding peak in mfsclean: ‘Algorithm1’|’CASA’|’ARL’, Default is ARL.
• prefix – Prefix to log messages to provide context

Returns:

clean component image, residual image

overlapIndices(res, psf, peakx, peaky)

Find the indices where two arrays overlap

Parameters:

• res –
• psf –
• peakx – peak in x
• peaky – peak in y

Returns:

(limits in a1, limits in a2)

spheroidal_function(vnu)

Evaluates the PROLATE SPHEROIDAL WAVEFUNCTION

m=6, alpha = 1 from Schwab, Indirect Imaging (1984). This is one factor in the basis function.

update_moment_model(m_model, scalestack, lhs, rhs, gain, mscale, mval)

Update model with an appropriately scaled and centered blob for each moment

update_scale_moment_residual(smresidual, ssmmpsf, lhs, rhs, gain, mscale, mval)

Update residual by subtracting the effect of model update for each moment

Calibration

Solvers

Functions to solve for antenna/station gain

This uses an iterative substitution algorithm due to Larry D’Addario c 1980’ish. Used in the original VLA Dec-10 Antsol.

For example:

gtsol = solve_gaintable(vis, originalvis, phase_only=True, niter=niter, crosspol=False, tol=1e-6)
vis = apply_gaintable(vis, gtsol, inverse=True)

solution_residual_matrix(gain, x, xwt)

Calculate residual across all baselines of gain for point source equivalent visibilities

Parameters:

• gain – gain [nant, …]
• x – Point source equivalent visibility [nant, …]
• xwt – Point source equivalent weight [nant, …]

Returns:

residual[…]

solution_residual_scalar(gain, x, xwt)

Calculate residual across all baselines of gain for point source equivalent visibilities

Parameters:

• gain – gain [nant, …]
• x – Point source equivalent visibility [nant, …]
• xwt – Point source equivalent weight [nant, …]

Returns:

residual[…]

solution_residual_vector(gain, x, xwt)

Calculate residual across all baselines of gain for point source equivalent visibilities

Vector case i.e. off-diagonals of gains are zero

Parameters:

• gain – gain [nant, …]
• x – Point source equivalent visibility [nant, …]
• xwt – Point source equivalent weight [nant, …]

Returns:

residual[…]

solve_antenna_gains_itsubs_matrix(gain, gwt, x, xwt, niter=30, tol=1e-08, phase_only=True, refant=0)

Solve for the antenna gains using full matrix expressions

x(antenna2, antenna1) = gain(antenna1) conj(gain(antenna2))

See Appendix D, section D.1 in:

J. P. Hamaker, “Understanding radio polarimetry - IV. The full-coherency analogue of scalar self-calibration: Self-alignment, dynamic range and polarimetric fidelity,” Astronomy and Astrophysics Supplement Series, vol. 143, no. 3, pp. 515–534, May 2000.

Parameters:

• gain – gains
• gwt – gain weight
• x – Equivalent point source visibility[nants, nants, …]
• xwt – Equivalent point source weight [nants, nants, …]
• niter – Number of iterations
• tol – tolerance on solution change
• phase_only – Do solution for only the phase? (default True)
• refant – Reference antenna for phase (default=0.0)

Returns:

gain [nants, …], weight [nants, …]

solve_antenna_gains_itsubs_scalar(gain, gwt, x, xwt, niter=30, tol=1e-08, phase_only=True, refant=0)

Solve for the antenna gains

x(antenna2, antenna1) = gain(antenna1) conj(gain(antenna2))

This uses an iterative substitution algorithm due to Larry D’Addario c 1980’ish (see ThompsonDaddario1982 Appendix 1). Used in the original VLA Dec-10 Antsol.

Parameters:

• gain – gains
• gwt – gain weight
• x – Equivalent point source visibility[nants, nants, …]
• xwt – Equivalent point source weight [nants, nants, …]
• niter – Number of iterations
• tol – tolerance on solution change
• phase_only – Do solution for only the phase? (default True)
• refant – Reference antenna for phase (default=0.0)

Returns:

gain [nants, …], weight [nants, …]

solve_antenna_gains_itsubs_vector(gain, gwt, x, xwt, niter=30, tol=1e-08, phase_only=True, refant=0)

Solve for the antenna gains using full matrix expressions

x(antenna2, antenna1) = gain(antenna1) conj(gain(antenna2))

See Appendix D, section D.1 in:

J. P. Hamaker, “Understanding radio polarimetry - IV. The full-coherency analogue of scalar self-calibration: Self-alignment, dynamic range and polarimetric fidelity,” Astronomy and Astrophysics Supplement Series, vol. 143, no. 3, pp. 515–534, May 2000.

Parameters:

• gain – gains
• gwt – gain weight
• x – Equivalent point source visibility[nants, nants, …]
• xwt – Equivalent point source weight [nants, nants, …]
• niter – Number of iterations
• tol – tolerance on solution change
• phase_only – Do solution for only the phase? (default True)
• refant – Reference antenna for phase (default=0.0)

Returns:

gain [nants, …], weight [nants, …]

solve_from_X(gt: data_models.memory_data_models.GainTable, x: numpy.ndarray, xwt: numpy.ndarray, chunk, crosspol, niter, phase_only, tol, npol) → data_models.memory_data_models.GainTable

Solve for gains from the point source equivalents

Parameters:

• gt –
• x – point source visibility
• xwt – point source weight
• chunk – which chunk of the gaintable?
• crosspol –
• niter –
• phase_only –
• tol –
• npol –

Returns:

Fourier transform

FFT support

FFT support functions

extract_mid(a, npixel)

Extract a section from middle of a map

Suitable for zero frequencies at npixel/2. This is the reverse operation to pad.

Note

Only the two innermost axes are transformed

Parameters:

• npixel – desired size of the section to extract
• a – grid from which to extract

extract_oversampled(a, xf, yf, kernel_oversampling, kernelwidth)

Extract the (xf-th,yf-th) w-kernel from the oversampled parent

Offsets are suitable for correcting of fractional coordinates, e.g. an offset of (xf,yf) results in the kernel for an (-xf,-yf) sub-grid offset.

We do not want to make assumptions about the source grid’s symmetry here, which means that the grid’s side length must be at least kernel_oversampling*(npixel+2) to contain enough information in all circumstances

Parameters:

• xf –
• yf –
• a – grid from which to extract
• kernel_oversampling – oversampling factor
• kernelwidth – size of section

fft(a)

Fourier transformation from image to grid space

Note

If there are four axes then the last outer axes are not transformed

Parameters:a – image in lm coordinate space Returns:uv grid

ifft(a)

Fourier transformation from grid to image space

Note

If there are four axes then the last outer axes are not transformed

Parameters:a – uv grid to transform Returns:an image in lm coordinate space

pad_mid(ff, npixel)

Pad a far field image with zeroes to make it the given size.

Effectively as if we were multiplying with a box function of the original field’s size, which is equivalent to a convolution with a sinc pattern in the uv-grid.

Note

Only the two innermost axes are transformed

This function does not handle odd-sized dimensions properly

Parameters:

• ff – The input far field. Should be smaller than npixelxnpixel.
• npixel – The desired far field size

Convolutional Gridding

Imaging is based on used of the FFT to perform Fourier transforms efficiently. Since the observed visibility data_models do not arrive naturally on grid points, the sampled points are resampled on the FFT grid using a convolution function to smear out the sample points. The resulting grid points are then FFT’ed. The result can be corrected for the gridding convolution function by division in the image plane of the transform.

This approach may be extended to include image plane effect such as the w term and the antenna/station primary beam.

This module contains functions for performing the gridding process and the inverse degridding process.

anti_aliasing_calculate(shape, oversampling=1, support=3)

Compute the prolate spheroidal anti-aliasing function

The kernel is to be used in gridding visibility data onto a grid on for degridding from a grid. The gridding correction function (gcf) is used to correct the image for decorrelation due to gridding.

Return the 2D grid correction function (gcf), and the convolving kernel (kernel

See VLA Scientific Memoranda 129, 131, 132 :param shape: (height, width) pair :param oversampling: Number of sub-samples per grid pixel :param support: Support of kernel (in pixels) width is 2*support+2

convert_to_tuple3(x_ary)

Numba cannot do this conversion itself. Hardcode 3 for speed

convolutional_degrid(kernel_list, vshape, uvgrid, vuvwmap, vfrequencymap)

Convolutional degridding with frequency and polarisation independent

Takes into account fractional uv coordinate values where the GCF is oversampled

Parameters:

• kernel_list – list of oversampled convolution kernel
• vshape – Shape of visibility
• uvgrid – The uv plane to de-grid from
• vuvwmap – function to map uvw to grid fractions
• vfrequencymap – function to map frequency to image channels

Returns:

Array of visibilities.

convolutional_grid(kernel_list, uvgrid, vis, visweights, vuvwmap, vfrequencymap)

Grid after convolving with frequency and polarisation independent gcf

Takes into account fractional uv coordinate values where the GCF is oversampled

Parameters:

• kernel_list – List of oversampled convolution kernels
• uvgrid – Grid to add to [nchan, npol, npixel, npixel]
• vis – Visibility values
• visweights – Visibility weights
• vuvwmap – map uvw to grid fractions
• vfrequencymap – map frequency to image channels

Returns:

uv grid[nchan, npol, ny, nx], sumwt[nchan, npol]

coordinateBounds(npixel)

Returns lowest and highest coordinates of an image/grid given:

1. Step size is \(1/npixel\):

   \[\frac{high-low}{npixel-1} = \frac{1}{npixel}\]

2. The coordinate \(\lfloor npixel/2\rfloor\) falls exactly on zero:

   \[low + \left\lfloor\frac{npixel}{2}\right\rfloor * (high-low) = 0\]

This is the coordinate system for shifted FFTs.

coordinates(npixel: int)

1D array which spans [-.5,.5[ with 0 at position npixel/2

coordinates2(npixel: int)

Two dimensional grids of coordinates spanning -1 to 1 in each dimension

1. a step size of 2/npixel and
2. (0,0) at pixel (floor(n/2),floor(n/2))

coordinates2Offset(npixel: int, cx: int, cy: int)

Two dimensional grids of coordinates centred on an arbitrary point.

This is used for A and w beams.

1. a step size of 2/npixel and
2. (0,0) at pixel (cx, cy,floor(n/2))

frac_coord(npixel, kernel_oversampling, p)

Compute whole and fractional parts of coordinates, rounded to kernel_oversampling-th fraction of pixel size

The fractional values are rounded to nearest 1/kernel_oversampling pixel value. At fractional values greater than (kernel_oversampling-0.5)/kernel_oversampling coordinates are rounded to next integer index.

Parameters:

• npixel – Number of pixels in total
• kernel_oversampling – Fractional values to round to
• p – Coordinate in range [-.5,.5[

grdsf(nu)

Calculate PSWF using an old SDE routine re-written in Python

Find Spheroidal function with M = 6, alpha = 1 using the rational approximations discussed by Fred Schwab in ‘Indirect Imaging’. This routine was checked against Fred’s SPHFN routine, and agreed to about the 7th significant digit. The gridding function is (1-NU**2)*GRDSF(NU) where NU is the distance to the edge. The grid correction function is just 1/GRDSF(NU) where NU is now the distance to the edge of the image.

gridder(uvgrid, vis, xs, ys, kernel=array([[1.]]), kernel_ixs=None)

Grids visibilities at given positions. Convolution kernels are selected per visibility using kernel_ixs.

Parameters:

• uvgrid – Grid to update (two-dimensional complex array)
• vis – Visibility values (one-dimensional complex array)
• xs – Visibility position (one-dimensional int array)
• ys – Visibility values (one-dimensional int array)
• kernel – Convolution kernel (minimum two-dimensional complex array). If the kernel has more than two dimensions, additional indices must be passed in kernel_ixs. Default: Fixed one-pixel kernel with value 1.
• kernel_ixs – Map of visibilities to kernel indices (maximum two-dimensional int array). Can be omitted if kernel requires no indices, and can be one-dimensional if only one index is needed to identify kernels

gridder_numba(uvgrid, vis, xs, ys, kernel=array([[1.]]), kernel_ixs=None)

Grids visibilities at given positions. Convolution kernels are selected per visibility using kernel_ixs.

Parameters:

• uvgrid – Grid to update (two-dimensional complex array)
• vis – Visibility values (one-dimensional complex array)
• xs – Visibility position (one-dimensional int array)
• ys – Visibility values (one-dimensional int array)
• kernel – Convolution kernel (minimum two-dimensional complex array). If the kernel has more than two dimensions, additional indices must be passed in kernel_ixs. Default: Fixed one-pixel kernel with value 1.
• kernel_ixs – Map of visibilities to kernel indices (maximum two-dimensional int array). Can be omitted if kernel requires no indices, and can be one-dimensional if only one index is needed to identify kernels

visibility_recentre(uvw, dl, dm)

Compensate for kernel re-centering - see w_kernel_function.

Parameters:

• uvw – Visibility coordinates
• dl – Horizontal shift to compensate for
• dm – Vertical shift to compensate for

Returns:

Visibility coordinates re-centrered on the peak of their w-kernel

w_beam(npixel, field_of_view, w, cx=None, cy=None, remove_shift=False)

W beam, the fresnel diffraction pattern arising from non-coplanar baselines

Parameters:

• npixel – Size of the grid in pixels
• field_of_view – Field of view
• w – Baseline distance to the projection plane
• cx – location of delay centre def :npixel//2
• cy – location of delay centre def :npixel//2
• remove_shift – Remove overall phase shift at the centre of the image

Returns:

npixel x npixel array with the far field

weight_gridding(shape, visweights, vuvwmap, vfrequencymap, vpolarisationmap=None, weighting='uniform')

Reweight data using one of a number of algorithms

Parameters:

• shape –
• visweights – Visibility weights
• vuvwmap – map uvw to grid fractions
• vfrequencymap – map frequency to image channels
• vpolarisationmap – map polarisation to image polarisation
• weighting – ‘’ | ‘uniform’

Returns:

visweights, density, densitygrid

Imaging

Parameters

Functions that aid definition of fourier transform processing.

get_frequency_map(vis, im: data_models.memory_data_models.Image = None)

Map channels from visibilities to image

get_kernel_list(vis: data_models.memory_data_models.Visibility, im: data_models.memory_data_models.Image, **kwargs)

Get the list of kernels, one per visibility

get_polarisation_map(vis: data_models.memory_data_models.Visibility, im: data_models.memory_data_models.Image = None)

Get the mapping of visibility polarisations to image polarisations

get_rowmap(col, ucol=None)

Map to unique cols

Parameters:

• col – Data column
• ucol – Unique values in col

get_uvw_map(vis: data_models.memory_data_models.Visibility, im: data_models.memory_data_models.Image, padding=2)

Get the generators that map channels uvw to pixels

Parameters:

• vis –
• im –
• padding –

Returns:

uvw mode, shape, padding, uvw mapping

standard_kernel_list(vis: data_models.memory_data_models.Visibility, shape, oversampling=8, support=3)

Return a generator to calculate the standard visibility kernel

Parameters:

• vis – visibility
• shape – tuple with 2D shape of grid
• oversampling – Oversampling factor
• support – Support of kernel

Returns:

Function to look up gridding kernel

w_kernel_list(vis: data_models.memory_data_models.Visibility, im: data_models.memory_data_models.Image, oversampling=1, wstep=50.0, kernelwidth=16, **kwargs)

Calculate w convolution kernels

Uses create_w_term_like to calculate the w screen. This is exactly as wstacking does.

Returns (indices to the w kernel for each row, kernels)

Each kernel has axes [centre_v, centre_u, offset_v, offset_u]. We currently use the same convolution function for all channels and polarisations. Changing that behaviour would require modest changes here and to the gridding/degridding routines.

Parameters:

• im –
• kernelwidth –
• vis – visibility
• oversampling – Oversampling factor
• wstep – Step in w between cached functions

Returns:

(indices to the w kernel for each row, kernels)

Util

Array Functions

Useful array functions.

average_chunks(arr, wts, chunksize)

Average the array arr with weights by chunks

Array len does not have to be multiple of chunksize

This version is optimised for plain numpy. It is roughly ten times faster that average_chunks_jit when used without numba jit. It cannot (yet) be used with numba because the add.reduceat is not support in numba 0.31

Parameters:

• arr – 1D array of values
• wts – 1D array of weights
• chunksize – averaging size

Returns:

1D array of averaged data_models, 1d array of weights

average_chunks2(arr, wts, chunksize)

Average the two dimensional array arr with weights by chunks

Array len does not have to be multiple of chunksize.

Parameters:

• arr – 2D array of values
• wts – 2D array of weights
• chunksize – 2-tuple of averaging region e.g. (2,3)

Returns:

2D array of averaged data_models, 2d array of weights

average_chunks_jit(arr, wts, chunksize)

Average the array arr with weights by chunks

Array len does not have to be multiple of chunksize

This is a version written for numba. When used with numba.jit, it’s about 25 - 30% faster than the numpy version without jit.

Parameters:

• arr – 1D array of values
• wts – 1D array of weights
• chunksize – averaging size

Returns:

1D array of averaged data_models, 1d array of weights

insert_array(im, x, y, flux, bandwidth=1.0, support=7, insert_function=<function insert_function_L>)

Insert point into image using specified function

Parameters:

• im – Image
• x – x in float pixels
• y – y in float pixels
• flux – Flux[nchan, npol]
• bandwidth – Support of data in uv plane
• support – Support of function in image space
• insert_function – insert_function_L or insert_function_Sinc or insert_function_pswf

Returns:

tukey_filter(x, r)

Calculate the Tukey (tapered cosine) filter

See e.g. https://uk.mathworks.com/help/signal/ref/tukeywin.html

Parameters:

• x – x coordinate (float)
• r – transition point of filter (float)

Returns:

Value of filter for x

Coordinate Support

Coordinate support

We follow the casa definition of coordinate systems http://casa.nrao.edu/Memos/CoordConvention.pdf :

UVW is a right-handed coordinate system, with W pointing towards the source, and a baseline convention of \(ant2 - ant1\) where \(index(ant1) < index(ant2)\). Consider an XYZ Celestial coordinate system centered at the location of the interferometer, with \(X\) towards the East, \(Z\) towards the NCP and \(Y\) to complete a right-handed system. The UVW coordinate system is then defined by the hour-angle and declination of the phase-reference direction such that

1. when the direction of observation is the NCP (ha=0,dec=90), the UVW coordinates are aligned with XYZ,
2. V, W and the NCP are always on a Great circle,
3. when W is on the local meridian, U points East
4. when the direction of observation is at zero declination, an hour-angle of -6 hours makes W point due East.

The \((l,m,n)\) coordinates are parallel to \((u,v,w)\) such that \(l\) increases with Right-Ascension (or increasing longitude coordinate), \(m\) increases with Declination, and \(n\) is towards the source. With this convention, images will have Right Ascension increasing from Right to Left, and Declination increasing from Bottom to Top.

baselines(ants_uvw)

Compute baselines in uvw co-ordinate system from uvw co-ordinate system station positions

Parameters:ants_uvw – (u,v,w) co-ordinates of antennas in array

lmn_to_skycoord(lmn, phasecentre: astropy.coordinates.sky_coordinate.SkyCoord)

Convert l,m,n coordinate system + phascentre to astropy sky coordinate relative to a phase centre.

The l,m,n is a RHS coordinate system with * its origin on the sky sphere * m,n and the celestial north on the same plane * l,m a tangential plane of the sky sphere

Note that this means that l increases east-wards

simulate_point(dist_uvw, l, m)

Simulate visibilities for unit amplitude point source at direction cosines (l,m) relative to the phase centre.

This includes phase tracking to the centre of the field (hence the minus 1 in the exponent.)

Note that point source is delta function, therefore the FT relationship becomes an exponential, evaluated at (uvw.lmn)

Parameters:

• dist_uvw – \((u,v,w)\) distribution of projected baselines (in wavelengths)
• l – horizontal direction cosine relative to phase tracking centre
• m – orthogonal directon cosine relative to phase tracking centre

skycoord_to_lmn(pos: astropy.coordinates.sky_coordinate.SkyCoord, phasecentre: astropy.coordinates.sky_coordinate.SkyCoord)

Convert astropy sky coordinates into the l,m,n coordinate system relative to a phase centre.

The l,m,n is a RHS coordinate system with * its origin on the sky sphere * m,n and the celestial north on the same plane * l,m a tangential plane of the sky sphere

Note that this means that l increases east-wards

uvw_to_xyz(uvw, ha, dec)

Rotate \((x,y,z)\) positions relative to a sky position at \((ha, dec)\) to earth coordinates. Can be used for both antenna positions as well as for baselines.

Hour angle and declination can be given as single values or arrays of the same length. Angles can be given as radians or astropy quantities with a valid conversion.

Parameters:

• uvw – \((u,v,w)\) co-ordinates of antennas in array
• ha – hour angle of phase tracking centre (\(ha = ra - lst\))
• dec – declination of phase tracking centre

uvw_transform(uvw, transform_matrix)

Transforms UVW baseline coordinates such that the image is transformed with the given matrix. Will require kernels to be suitably transformed to work correctly.

Reference: Sault, R. J., L. Staveley-Smith, and W. N. Brouw. “An approach to interferometric mosaicing.” Astronomy and Astrophysics Supplement Series 120 (1996): 375-384.

Parameters:

• uvw – \((u,v,w)\) distribution of projected baselines (in wavelengths)
• transform_matrix – 2x2 matrix for image transformation

Returns:

New baseline coordinates

visibility_shift(uvw, vis, dl, dm)

Shift visibilities by the given image-space distance. This is based on simple FFT laws. It will require kernels to be suitably shifted as well to work correctly.

Parameters:

• uvw –
• vis – \((u,v,w)\) distribution of projected baselines (in wavelengths)
• vis – Input visibilities
• dl – Horizontal shift distance as directional cosine
• dm – Vertical shift distance as directional cosine

Returns:

New visibilities

xyz_at_latitude(local_xyz, lat)

Rotate local XYZ coordinates into celestial XYZ coordinates. These coordinate systems are very similar, with X pointing towards the geographical east in both cases. However, before the rotation Z points towards the zenith, whereas afterwards it will point towards celestial north (parallel to the earth axis).

Parameters:

• lat – target latitude (radians or astropy quantity)
• local_xyz – Array of local XYZ coordinates

Returns:

Celestial XYZ coordinates

xyz_to_baselines(ants_xyz, ha_range, dec)

Calculate baselines in \((u,v,w)\) co-ordinate system for a range of hour angles (i.e. non-snapshot observation) to create a uvw sampling distribution

Parameters:

• ants_xyz – \((x,y,z)\) co-ordinates of antennas in array
• ha_range – list of hour angle values for astronomical source as function of time
• dec – declination of astronomical source [constant, not \(f(t)\)]

xyz_to_uvw(xyz, ha, dec)

Rotate \((x,y,z)\) positions in earth coordinates to \((u,v,w)\) coordinates relative to astronomical source position \((ha, dec)\). Can be used for both antenna positions as well as for baselines.

Hour angle and declination can be given as single values or arrays of the same length. Angles can be given as radians or astropy quantities with a valid conversion.

Parameters:

• xyz – \((x,y,z)\) co-ordinates of antennas in array
• ha – hour angle of phase tracking centre (\(ha = ra - lst\))
• dec – declination of phase tracking centre.