<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Outline
Glossary
5. Imaging
Previous
Step1:  Import section specific modules
Step2:  5.5 The Break Down of the Small Angle Approximation and the W-Term
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
Step3:  Figure
Step4:  Figure
Step5:  Figure
Step6:  Figure
Step7:  Figure
Step8:  Figure
Step9:  Figure
Step10:  Figure
Step11:  Figure
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from IPython.display import HTML 
HTML('../style/course.css') #apply general CSS
Explanation: Outline
Glossary
5. Imaging
Previous: 5.4 Imaging weights 
Next: 5.5 References and further reading
Import standard modules:
End of explanation
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)
HTML('../style/code_toggle.html')
Explanation: Import section specific modules:
End of explanation
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
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()
Explanation: 5.5 The Break Down of the Small Angle Approximation and the W-Term
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.
<img src="figures/orthogonal_projection_difference.png" alt="Smiley face" width="512">
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.
5.5.1 Coplanar Sampling
Consider the following two hypothetical arrays: a perfectly coplanar array that only has baselines along the east-west direction, and a second perfectly coplanar two-dimensional array with some baselines in a non-east-west direction. An array is coplanar if there exists (approximately) a 2-D plane all visibility measurement lie on in the 3-D visibility space.
End of explanation
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()
Explanation: 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:
1. Improved u,v coverage at lower declinations, as plotted below.
2. Recall that the interferometer response is maximum when the phase-reference centre is orthogonal to the baseline. At lower observation angles it is desirable to have baseline components that are not aligned from east-to-west.
End of explanation
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()
Explanation: 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 co-planar, 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. Alos, if an observation is sufficiently short, called a snapshot observation, then the rotation of the Earth is short enough to approximate the measurements as coplanar.
End of explanation
Image(filename="figures/tilted_interferometer.png")
Explanation: 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}
5.5.2 Non-coplanar Sampling
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.
End of explanation
Image(filename="figures/vla_uncorrected.png")
Explanation: 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.
End of explanation
Image(filename="figures/vla_wproj.png")
Explanation: 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.
End of explanation
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()
Explanation: Figure: W-projection image of the 8 hour observation of the supernova reminant G55.7+3.4 on the JVLA in D-configuration.
End of explanation
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()
Explanation: 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.
End of explanation
Image(filename="figures/coplanar-faceting.png")
Explanation: 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.
5.3.3 Correcting Non-coplanar Baselines Effects
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 <cite data-cite='taylor1999synthesis'>Synthesis Imaging in Radio Astronomy II</cite> &#10548; 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:
End of explanation
Image(filename="figures/non-coplanar-faceting.png")
Explanation: 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.
End of explanation
<END_TASK>