Import standard modules:
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from IPython.display import HTML
HTML('../style/course.css') #apply general CSS
Import section specific modules:
In [ ]:
from IPython.display import Image
from mpl_toolkits.mplot3d import Axes3D
import track_simulator
from astropy.io import fits
import aplpy
#Disable astropy/aplpy logging
import logging
logger0 = logging.getLogger('astropy')
logger0.setLevel(logging.CRITICAL)
logger1 = logging.getLogger('aplpy')
logger1.setLevel(logging.CRITICAL)
In [ ]:
HTML('../style/code_toggle.html')
Up to this point we used a resampling step and the Fast Fourier Transform to move between the image and visibility domains. Recall that we used the following simplified Fourier relationship to justify this synthesis process:
\begin{equation} \begin{split} V(u,v) &= \int_\text{sky}{I(l,m)e^{-2\pi i/\lambda(\vec{b}\dot(\vec{s}-\vec{s}_0))}}dS\\ &= \int\int{I(l,m)e^{-2\pi i/\lambda(ul+vm+w(\sqrt{1-l^2-m^2}-1))}}\frac{dldm}{\sqrt{1-l^2-m^2}}\\ &\approx\int\int{I(l,m)e^{-2\pi i/\lambda(ul+vm)}}dldm\\ \end{split} \end{equation}The last approximation to the model is just a Fourier transform by definition and is the one used when we were imaging up to this point. However, the more accurate version that relates the measurement to the brightness distribution along the celestial sphere is not the classical Fourier transform. The approximation is only valid when $n - 1 = \sqrt{1-l^2-m^2} - 1 \ll 1$ (ie. images of small regions of the sky) and/or $w \approx 0$ (the array is coplanar). Here $(n-1)$ is the projection height difference between the planar approximation tangent to the celestial sphere and a source's true position on the sphere, see the illustration below.
Figure: The direction cosines (here $l$ is plotted against $n$) lie along the unit celestial sphere. $n$ is given by $n=\sqrt{1-l^2-m^2}$. If the projection pole (tangent point of the image) is at the same point as the phase reference centre, $n_0 = 1$. The total error between the orthogonal (SIN) projection of the source onto the tangent image plane and the source position on the celestial sphere is given as $\epsilon=(n-n_0)=(\sqrt{1-l^2-m^2} - 1)$.
Under the assumptions of a narrow field of view and coplanar measurements it is valid to use the FFT to construct a planar approximation to the sky. This section discusses the problem of wide-field imaging using non-coplanar baselines that arises when these assumptions are broken.
Consider the following two hypothetical arrays: a perfectly flat array that only has baselines along the east-west direction, and a second perfectly flat two-dimensional array with some baselines in a non-east-west direction.
In [ ]:
NO_ANTENNA = 4
NO_BASELINES = NO_ANTENNA * (NO_ANTENNA - 1) / 2 + NO_ANTENNA
CENTRE_CHANNEL = 1e9 / 299792458 #Wavelength of 1 GHz
#Create a perfectly planar array with both a perfectly East-West baseline and 2 2D baselines
ENU_2d = np.array([[5,0,0],
[-5,0,0],
[10,0,0],
[0,23,0]]);
ENU_ew = np.array([[5,0,0],
[-5,0,0],
[10,0,0],
[0,0,0]]);
ARRAY_LATITUDE = 30 #Equator->North
ARRAY_LONGITUDE = 0 #Greenwitch->East, prime -> local meridian
In [ ]:
fig = plt.figure(figsize=(10, 5))
ax=fig.add_subplot(121)
ax.set_title("2D Array")
ax.plot(ENU_2d[:,0],ENU_2d[:,1],"ko")
ax.set_xlabel("East")
ax.set_ylabel("North")
ax.set_xlim(-30,30)
ax.set_ylim(-30,30)
ax=fig.add_subplot(122)
ax.set_title("East-west array")
ax.plot(ENU_ew[:,0],ENU_ew[:,1],"ko")
ax.set_xlabel("East")
ax.set_ylabel("North")
ax.set_xlim(-30,30)
ax.set_ylim(-30,30)
plt.show()
Figure: ENU coordinates for two hypothetical flat arrays: a 2D array and an east-west array
The two-dimensional interferometer has two major advantages over its one-dimensional east-west counterpart:
In [ ]:
DECLINATION = 0
T_OBS = 12
T_INT = 1/60.0
uw_2hr_2d = track_simulator.sim_uv(0.0,DECLINATION,T_OBS,T_INT,ENU_2d,ARRAY_LATITUDE,False)/CENTRE_CHANNEL
uv_2hr_ew = track_simulator.sim_uv(0.0,DECLINATION,T_OBS,T_INT,ENU_ew,ARRAY_LATITUDE,False)/CENTRE_CHANNEL
fig = plt.figure(figsize=(10, 5))
ax=fig.add_subplot(121)
ax.set_title("2D Array")
ax.plot(uw_2hr_2d[:,0],uw_2hr_2d[:,1],'k.')
ax.set_xlabel("u")
ax.set_ylabel("v")
ax.set_xlim(-10,10)
ax.set_ylim(-10,10)
ax=fig.add_subplot(122)
ax.set_title("East-west Array")
ax.plot(uv_2hr_ew[:,0],uv_2hr_ew[:,1],'k.')
ax.set_xlabel("u")
ax.set_ylabel("v")
ax.set_xlim(-10,10)
ax.set_ylim(-10,10)
plt.show()
Figure: u,v coverage at declinatin $\delta=0$ for both a 2D and east-west array
The one drawback to using these two-dimensional layouts is that the measurements taken over the duration of the observation do not remain coplanar, even though the array layout is perfectly flat. The uvw tracks and their projections are plotted in 3-space below to illustrate this. This is opposed to the tracks created by the east-west interferometer which all remain in the same plane parallel to the Earth's equator. Alas, if an observation is sufficiently short, called a snapshot observation, then the rotation of the Earth is short enough to approximate the measurements of a two-dimensional interferometer as coplanar.
In [ ]:
DECLINATION = 45
T_INT = 1/60.0
T_OBS = 12
uvw_2d = track_simulator.sim_uv(0.0,DECLINATION,T_OBS,T_INT,ENU_2d,ARRAY_LATITUDE,False)/CENTRE_CHANNEL
uvw_ew = track_simulator.sim_uv(0.0,DECLINATION,T_OBS,T_INT,ENU_ew,ARRAY_LATITUDE,False)/CENTRE_CHANNEL
fig=plt.figure(figsize=(10, 5))
ax=fig.add_subplot(121,projection='3d')
ax.set_title("2D Array")
ax.view_init(elev=10, azim=160)
ax.plot(uvw_2d[:,0],uvw_2d[:,1],uvw_2d[:,2],'k.')
ax.set_xlabel("u")
ax.set_ylabel("v")
ax.set_zlabel("w")
ax=fig.add_subplot(122,projection='3d')
ax.set_title("East-west array")
ax.view_init(elev=10, azim=160)
ax.plot(uvw_ew[:,0],uvw_ew[:,1],uvw_ew[:,2],'k.')
ax.set_xlabel("u")
ax.set_ylabel("v")
ax.set_zlabel("w")
plt.show()
fig = plt.figure(figsize=(10, 5))
ax=fig.add_subplot(121)
ax.set_title("2D Array")
ax.plot(uvw_2d[:,0],uvw_2d[:,1],'k.')
ax.set_xlabel("u")
ax.set_ylabel("v")
ax=fig.add_subplot(122)
ax.set_title("East-west array")
ax.plot(uvw_ew[:,0],uvw_ew[:,1],'k.')
ax.set_xlabel("u")
ax.set_ylabel("v")
plt.show()
Figure: $u,v,w$ tracks and their projections onto ($u,v,w=0$) for a 2D and east-west interferometer
When the measurement domain is sampled along a single plane, as is true for the east-west interferometer, then all $w$ can be written as the same linear combination of u and v: $w = \alpha{u}+\beta{v}$. Although this introduces a slight distortion of the u,v coordinates in the Fourier relationship between the sky and the measurements, the distorted relationship remains a valid two-dimensional Fourier transform. It can be stated as:
\begin{equation} \begin{split} V(u,v,w) &= \int\int{I(l,m)e^{-2\pi i/\lambda(ul' + vm')}\frac{dldm}{\sqrt{1-l^2-m^2}}}\\ l' &= l + \alpha(\sqrt{1-l^2-m^2} - 1)\\ m' &= m + \beta(\sqrt{1-l^2-m^2} - 1)\\ \end{split} \end{equation}The same can not be said for two-dimensional arrays. There is no fixed relationship between $w$ and $u,v$. Instead the relationship depends both on the time-variant zenithal and parallactic angles, and the $u,v$ coverage only remains co-planar for instantaneous observations, provided the array layout is approximately flat.
Neglecting the $w(n-1)$ term by synthesizing wide-field images with two-dimensional arrays, using a planar approximation, introduces a direction-dependent error in the measurement. This phase error depends on the height-difference between antennae, as is illustrated by the tilted interferometer below.
In [ ]:
Image(filename="figures/tilted_interferometer.png")
Figure: As the two figures show the projection of the source vector onto the two baselines are different for the coplanar and tilted interferometers. The phase for signals taken by co-planar interferometer baselines along some line of sight, $\vec{s}$, is given as $\phi = \frac{-2\pi i}{\lambda}(ul + vm)$, as opposed to tilted baselines that measure this same phase as $\phi_\text{tilt} = \frac{-2\pi i}{\lambda}{[ul + vm + w(n-1)]}$. The signal propagation delay is worse on the longest baselines and along the direction of sources far away from the phase centre of the interferometer.
It is important to realize that this phase term is purely geometric in origin; inserting a delay to correct $\Delta{w}$ for non-coplanar measurements only serves to correct the error in a single line of sight. In other words only the phase centre of the interferometer is changed by such a correction.
The small angle approximation $\sqrt{1+x} \approx 1+ \frac{x}{2}$ gives some intuition on how this phase error effects the brightness of sources away from the phase centre. It can be shown that: \begin{equation} V(u,v,w) \approx {\int\int{I(l,m)(e^{2\pi i/\lambda wl^2/2}e^{2\pi i /\lambda wm^2/2})e^{-2\pi i /\lambda(ul+vm)}\frac{dldm}{n}}} \end{equation}
Since $w$ can be rewritten as a complex relationship of $u,v$ and time-variant elevation and azimuth angles we expect to see a time- and baseline-variant shift in source position. This relative position shift also grows roughly quadratically with the source offsets in l and m. The images below show how sources are smeared over large areas during long observations first on the data captured with the JVLA and then on a simulated MeerKAT observation.
In [ ]:
Image(filename="figures/vla_uncorrected.png")
Figure: Uncorrected image of the 8 hour observation of the supernova reminant G55.7+3.4 on the JVLA in D-configuration. Notice the eliptical smearing around the point source.
In [ ]:
Image(filename="figures/vla_wproj.png")
Figure: W-projection image of the 8 hour observation of the supernova reminant G55.7+3.4 on the JVLA in D-configuration.
In [ ]:
gc1 = aplpy.FITSFigure('../data/fits/wterm/MeerKAT_6h60s_dec-30_10MHz_10chans_uniform_n3000_w0-image.fits')
cpx = gc1.pixel2world(256, 256)
gc1.recenter(cpx[0], cpx[1], radius=0.2)
gc1.show_colorscale(vmin=-0.2, vmax=1., cmap='viridis')
gc1.hide_axis_labels()
gc1.hide_tick_labels()
plt.title('MeerKAT Observation (Not Corrected)')
gc1.add_colorbar()
Figure: A quadrant of an image, not w-projection corrected, from a MeerKAT simulated observation, the phase centre is at the top right corner. Sources further from the phase centre have a larger amount of smearing compared to closer in.
In [ ]:
gc1 = aplpy.FITSFigure('../data/fits/wterm/MeerKAT_6h60s_dec-30_10MHz_10chans_uniform_n3000-image.fits')
cpx = gc1.pixel2world(256, 256)
gc1.recenter(cpx[0], cpx[1], radius=0.2)
gc1.show_colorscale(vmin=-0.2, vmax=1., cmap='viridis')
gc1.hide_axis_labels()
gc1.hide_tick_labels()
plt.title('MeerKAT Observation (W-Corrected)')
gc1.add_colorbar()
Figure: A quadrant of an image, w-projection corrected, from a MeerKAT simulated observation, the phase centre is at the top right corner. Sources far from the phase centre remain point-like when the correction is accounted for as compared to the un-corrected image above.
There are various ways the delay error introduced when discarding the $w(n-1)$ term during resampling and 2D Fast Fourier Transform can be corrected for, these include:
Full 3D transform: Similar to the 2D Direct Fourier Transform the Fourier transform can be computed for every element in a cube of $l,m,n$ values. The sky lies along the unit sphere within this cube. See Perley's discussion in Synthesis Imaging in Radio Astronomy II ⤴ for a full derivation of this usually computationally and memory prohibitive technique.
Snapshot imaging: As alluded to earlier the visibility measurements taken during very short observations are co-planar, assuming the physical array lies on a flat plane. During each observation the $l,m$ coordnates are slightly distorted and the images have to be interpolated to the same coordinates before the images can be averaged into a single map of the sky.
Facet imaging: In facet imaging the goal is to drive the $(n-1)$ factor down to 0; satisfying the narrow-field assumption that makes the 2D Fourier inversion valid. There are a few ways in which the sky can be split into smaller images, but the classical faceting approach is to tile the celestial sphere with many small tangent images, approximating the sky by a polyhedron.
The algorithm behind tangent (polyhedron) facet imaging is simple to implement. First the sky is recentred at the image centres $l_i,m_i$ of each of the narrow-field facets, by phase rotating the measured visibilities. Each of the facet-images is then rotated to be tangent to the sky sphere. As the Fourier transform preserves rotations, the facets can be tilted by rotating the u,v coordinates of the measurements to the tracks that would have been produced if the interferometer was pointing at $\alpha_i,\delta_i$, instead of the original phase tracking centre. Let $(l_\Delta,m_\Delta,n_\Delta) = (l_i-l_0,m_i-m_0,n_i-n_0)$, then:
\begin{equation} \begin{split} V(u,v,w)&\approx\int{\int{B(l-l_i,m-m_i,n-n_i)e^{-2{\pi}i[u(l-l_i)+v(m-m_i)+w(n-n_i)]}\frac{dldm}{n}}}\\ &\approx\int{\int{B(l-l_i,m-m_i,n-n_i)e^{-2{\pi}i[u(l-l_0-l_\Delta)+v(m-m_0-m_\Delta)+w(n-n_0-n_\Delta)]}\frac{dldm}{n}}}\\ &\approx\left[\int{\int{B(l-l_0,m-m_0,n-n_0)e^{-2{\pi}i[u(l-l_0)+v(m-m_0)+w(n-n_0)]}\frac{dldm}{n}}}\right]e^{2{\pi}i[ul_\Delta,vm_\Delta,wn_\Delta]}\\ \end{split} \end{equation}Note that if only the phase rotation is performed without rotating the facet geometry the effective field of view of individual facets that are far away from the phase centre will decrease. In order to keep the projection error at the edge of all the facets comparable this means that the facets closer to the edge of the field must be significantly smaller, increasing the computational demands of such an approach. Instead if the facets are rotated to form a polyhedron around the celestial sphere the facets can all be the same size. A simple visual proof of this is given by the following two cartoons:
In [ ]:
Image(filename="figures/coplanar-faceting.png")
Figure: Only phase steering the visibilities to new phase centres without tilting the u,v,w coordinates to correspond to the new phase tracking centre significantly reduces the achievable field of view. Here instead each facet is parallel to the original tangent plane. As the new centre is taken further away from the original phase tracking centre the effective facet size must be shrunk down to achieve a comparable projection error at the edge of the synthesized facets.
In [ ]:
Image(filename="figures/non-coplanar-faceting.png")
Figure: When the geometry of the facets are rotated the facets can all be the same size
In order to make each facet tangent to the celestial sphere at $(l_i,m_i)$ it is necessary to employ the following rotation matrices to compute new u',v',w' coordinates after the visibilities have been phase shifted using the old u,v,w coordinates. \begin{equation} \left[\begin{array}{c} u'\\ v'\\ w'\\ \end{array} \right] = R(\alpha_i,\delta_i)R^{T}(\alpha_0,\delta_0)\left[\begin{array}{c} u\\ v \\ w \\ \end{array} \right] \end{equation}
\begin{equation} R(\alpha,\delta) = \left[\begin{array}{c c c} -\sin{\alpha} & \cos{\alpha} & 0 \\ -\sin{\delta}\cos{\alpha} & -\sin{\delta}\sin{\alpha} & \cos{\delta}\\ \cos{\delta}\cos{\alpha} & \cos{\delta}\sin{\alpha} & \sin{\delta}\\ \end{array}\right] \end{equation}The number of facets roughly needed to satisfy the sampling criterion is given as:
\begin{equation} N_\text{facets} = \frac{2L\lambda}{{\xi}D^2}, \; \xi\ll{1} \end{equation}
$L$ is the magnitude of the longest baseline and $D$ is the diametre of the antenna apertures. $\xi$ is a quality factor that specifies the furthest seperation allowed for the facet images from the celestial sphere.
W-projection: W-projection eliminates $w$ from the phase term by relating all non-coplanar measurements to measurements taken at $w=0$. Employing the convolution theorem the following relationship between $V(u,v,w)$ and $V(u,v,0)$ can be obtained:
\begin{equation} \begin{split} V &= \int\int{I(l,m)e^{-2\pi i[ul+vm]}e^{-2\pi i[w(n-1)]}\frac{dldm}{n}}\\ V &= \int\int{I(l,m)e^{-2\pi i[ul+vm]}\mathcal{w}_w(l,m)\frac{dldm}{n}}\\ V(u,v,w) &= V(u,v,w=0)\circ\mathcal{W}_w(u,v)\\ \end{split} \end{equation}This says that any measurement can be related to the the same $u,v$ plane during the resampling step by picking the the relevant $w$-dependent filter from a stack of precomputed filters. Normally these w-filters are combined with an anti-aliasing filter, as used in the narrow-field imaging approach.
The support size of these $\mathcal{W}_w$ filters are dependent on the size of the image and given by the following relation:
\begin{equation} W_\text{sup} = \frac{4\pi w_\text{max}}{\lambda} \frac{D_\text{im}^2}{\sqrt{2-D_\text{im}^2}} \end{equation}Here $D_\text{im}$ is the diametre of the image and shows that the computational complexity of the method depends on the size of the image.
W-stacking: The alternative image-domain approach (called "w-stacking") multiplies the w-dependent phase directly into the image:
The number of discretized w-filters (w-projection) or image layers (w-stacking) are given by:
\begin{equation} N_\text{planes} = \frac{2\pi w_\text{max}(n_\text{edge}-1)}{\lambda_\text{min}\xi}, \; \xi\ll{1} \end{equation}The more discretized layers of w-values (as controlled by the quality control factor $\xi$) used during inversion, the more accurate the phase correction becomes.
Take home point: When imaging over field of views several degrees in size, especially if there are long baselines in the array you are observing with, you have to enable wide-field corrections in the imaging tool you're using. Check your imager's documentation on how to enable its wide-field mode.