This repository contains code to do post-correlation beamforming and imaging. It is designed for compact low-frequency arrays like the EDA2.
astropy
h5py
healpy
matplotlib
numpy
pyephem
pygdsm
A beamformed voltage beam is given by the sum of weights
The power beam is the voltage beam after squaring and averaging in time. when averaged with time (denote with
and so
Where
In terms of computations, the most efficient way to form a power beam, with order
is conceptually useful (and useful for visualization), but comes at even more computational expense and memory requirements.
However, post-correlation beamforming is incredibly flexible: a user can form as many beams as they feel like, and can re-point the beam in any desired direction. Also, as interferometers natively output visbility matrices, post-correlation beamforming can be a useful approach if a real-time voltage beamformer does not exist.
The post-correlation approach allows a compact summation notation to forming a grid of beams on the sky (coordinate subscripts
Summation here is implied over all indices that appear in an upper and lower index. So, we sum across indexes
The visibility matrix can itself be written in summation notation as
where
In Python (Numpy and Cupy), there is a einsum
command, which has slightly different syntax. The visibility matrix is compute via:
# v: voltage array with shape (N_ant, N_timesteps), complex dtype
# V: visibility matrix with shape (N_ant, N_ant), complex dtype
V = np.einsum('pt,qt->pq', v, v.conj(), optimize=True)
For a given frequency
# weights: weight array with shape (N_pix, N_pix, N_ant), complex dtype
# V: visibility matrix with shape (N_ant, N_ant), complex dtype
B = np.einsum('ijp,pq,ijq->ij', weights, V, weights.conj(), optimize=True)
Note that the optimize=True
argument is key to good performance, by optimizing the contraction order can dispatch to BLAS packages (see opt_einsum).
To phase up an array, we need to generate a set of geometric time delays,
The geometric delay is measured from a reference antenna, which we will set to be
-
Use
$(l, m)$ direction cosines of the target from the phase center. Useful for making 2D images around the phase center. Remember$(l, m)$ are cartesian coordinates, not spherical. -
Use the hourangle and declination
$(H, d)$ of the target from the phase center. Useful for pointing at a Celestial source with known RA/DEC, or za/az.
In both instances, the geometric delay is found from the W-component:
-
Local XYZ: Relative to the surface of the Earth, at the observatory site. Normally X=East, Y=North, Z=Up. That is: The Z-direction points toward zenith.
-
ECEF XYZ: Abstract coordinates relative to the center of the Earth. The Z direction points toward the the North celestial pole**. The X-axis is in the plane of the Equator and runs through
$0^\circ$ and$180^\circ$ degrees. The Y-axis runs through$90^\circ$ W to$90^\circ$ E.
** It's actually the 'international reference pole (IRP)', which I'm not 100% sure is the same as the NCP.
This is where you are pointing.
-
Local za/az: zenith angle (0 to
$90^\circ$ ) and azimuthal angle (full 360 degrees). Also called horizontal coordinates. -
Celestial RA/DEC: Right ascension and declination, as measured from the North Celestial pole.
When making an image with an interferometer, we have to project the 3D celestial sphere into points on a 2D plane. To do so, we choose a phase center pointing
In UVW coordinates, The W-axis is aligned along the pointing direction. It is extremely important to realize that later on the W-component of a baseline is equivalent to the geometric delay due to light travel time. The U and V axes are then in a plane.
-
Local UVW: The W-direction is chosen to be zenith, and then U and V axes can be lined up with local XY coordinates.
-
Celestial UVW: The W-direction is still chosen to be zenith, but we align the XY-plane with the North celestial pole.
When projecting from 3D spherical to 2D, the two plane axes are called
$l = cos(X)$ $m = cos(Y)$
Which both range from (-1, 1). Note these are cartesian, not spherical coordinates! The phase center is at the origin (0, 0). The pointing vector that corresponds to this is: