<jupyter_start><jupyter_text>From aperture masking to kernel-phaseAuthor: Dr. <NAME>In the beginning of this notebook, you will find a few examples which we will go through during the lecture. Further down, you will find some problems that you should try to solve on your own.The goal of this notebook is to recall some basics about Fourier plane imaging and then to develop an understanding of the aperture masking and kernel-phase interferometry techniques. The point-spread function (PSF)The object intensity distribution of a single star is a Dirac delta function:<jupyter_code>import astropy.io.fits as pyfits
import matplotlib.pyplot as plt
import numpy as np
OID = np.zeros((100, 100)) # an empty 100 x 100 image
OID[50, 50] = 1. # we put a point source with flux normalized to 1 in the center

# Plotting...
plt.imshow(OID, origin='lower')
plt.colorbar()
plt.show()<jupyter_output><empty_output><jupyter_text>Its Fourier transform is the Fourier transform of a Dirac delta function which is 1, i.e., has an absolute value of 1 and a phase value of 0 (remember, Fourier transforms are complex):<jupyter_code>OID_FT = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(OID))) # compute FT assuming that central pixel is origin

# Plotting...
f, ax = plt.subplots(1, 2, figsize=(2*6.4, 1*4.8))
p0 = ax[0].imshow(np.abs(OID_FT), origin='lower', vmin=0, vmax=1)
plt.colorbar(p0, ax=ax[0])
ax[0].set_title('Fourier amplitude')
p1 = ax[1].imshow(np.angle(OID_FT), origin='lower', vmin=-np.pi, vmax=np.pi)
plt.colorbar(p1, ax=ax[1])
ax[1].set_title('Fourier phase')
plt.show()<jupyter_output><empty_output><jupyter_text>NOTE: fftshift considers sz/2 as the image center, although technically the image center is (sz-1)/2.The object intensity distribution of a binary star are two Dirac delta functions:<jupyter_code>OID = np.zeros((100, 100)) # an empty 100 x 100 image
OID[50, 50] = 1. # we put a point source with flux normalized to 1 in the center
OID[45, 50] = 0.1 # we put a second point source with smaller flux somewhere near the primary

# Plotting...
plt.imshow(OID, origin='lower')
plt.colorbar()
plt.show()<jupyter_output><empty_output><jupyter_text>In the high contrast regime, its Fourier transform is a sine wave whose amplitude depends on the contrast and whose wavelength and direction depend on the separation:<jupyter_code>OID_FT = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(OID))) # compute FT assuming that central pixel is origin

# Plotting...
f, ax = plt.subplots(1, 2, figsize=(2*6.4, 1*4.8))
p0 = ax[0].imshow(np.abs(OID_FT), origin='lower')
plt.colorbar(p0, ax=ax[0])
ax[0].set_title('Fourier amplitude')
p1 = ax[1].imshow(np.angle(OID_FT), origin='lower')
plt.colorbar(p1, ax=ax[1])
ax[1].set_title('Fourier phase')
plt.show()<jupyter_output><empty_output><jupyter_text>The wavelength is given by 1/n, where n is the fractional separation with respect to the image size. In this case, the image size is 100 and the separation is 5, so that the fractional separation is 5/100 = 0.05 and the wavelength is 1/0.05 = 20. The amplitude of the wave is roughly the contrast, here 0.1/1 = 0.1 (in both amplitude and phase). The direction of the wave is given by the direction of the separation vector.The following data contains the image of a PSF reference star (i.e., a single star) observed with one of the world's largest telescopes. How does its OTF look like?<jupyter_code># Opening data from a FITS file...
hdul = pyfits.open('n0010.fits')
data = hdul[0].data
hdul.close()

# Plotting...
f, ax = plt.subplots(1, 2, figsize=(2*6.4, 1*4.8))
p0 = ax[0].imshow(data, origin='lower')
plt.colorbar(p0, ax=ax[0])
ax[0].set_title('Image')
p1 = ax[1].imshow(np.log10(np.abs(data)), origin='lower')
plt.colorbar(p1, ax=ax[1])
ax[1].set_title('Image in log stretch')
plt.show()<jupyter_output><empty_output><jupyter_text>Since the data contains only the image of a single star, its Fourier transform is directly equivalent to the optical transfer function of the telescope (remember the auto-correlation theorem):<jupyter_code>OTF = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(data))) # compute FT assuming that central pixel is origin

# Plotting...
f, ax = plt.subplots(1, 2, figsize=(2*6.4, 1*4.8))
p0 = ax[0].imshow(np.log10(np.abs(OTF)), origin='lower')
plt.colorbar(p0, ax=ax[0])
ax[0].set_title('Fourier amplitude in log stretch')
p1 = ax[1].imshow(np.angle(OTF), origin='lower')
plt.colorbar(p1, ax=ax[1])
ax[1].set_title('Fourier phase')
plt.show()<jupyter_output><empty_output><jupyter_text>The Fourier power (in logarithmic color scale) and the Fourier phase clearly show a hexagonal shape. The data comes from Keck! Aperture masking interferometry (AMI)As shown during the lecture, if we write the equations for the Fourier phase in matrix form, we can also search for linear combinations that cancel out the phase error terms by computing the left null space (or kernel) of the matrix A. Numerically, this can be done using a singular value decomposition:<jupyter_code>A = np.array([[-1, 1, 0], [0, -1, 1], [1, 0, -1]]) # this is the A matrix from the closure-phase equation
U, s, Vh = np.linalg.svd(A.T, full_matrices=True) # using numpy, we compute its singular value decomposition (SVD)
print('The unitary matrix U:')
print(U)
print('The singular values of A:')
print(s)
print('The unitary matrix Vh:')
print(Vh)<jupyter_output><empty_output><jupyter_text>The left null space (or kernel) of the matrix A is given by the columns of Vh that correspond to the zero singular values:<jupyter_code>ww = np.where(np.abs(s) < 1e-10)[0] # find the zero singular values of A
print(ww)
K = Vh[ww, :] # this is the kernel of A
print(K)<jupyter_output><empty_output><jupyter_text>The matrix K now contains basis vectors for the Fourier phase that cancel out the phase error terms. Here, there is only one such basis vector, and it is (0.577, 0.577, 0.577). This means that t*(0.577*phi1+0.577phi2+0.577phi3) is something that is independent of the phase error terms. We can run the check:<jupyter_code>np.dot(K, A)<jupyter_output><empty_output><jupyter_text>which returns the zero vector (+ some numerical error). Indeed, we found the left null space of A!Let's have a look at how interferometry increases the spatial resolution from $\sim\lambda/D$ to $\sim\lambda/(2D)$:<jupyter_code>ramp = np.arange(512)-512/2.
# print(ramp)
xx, yy = np.meshgrid(ramp, ramp)
dist = np.sqrt(xx**2+yy**2)

# Compute single dish pupil...
pupil_single = dist < 32
psf_single = np.abs(np.fft.fftshift(np.fft.fft2(np.fft.fftshift(pupil_single))))**2

# Plotting...
f, ax = plt.subplots(1, 2, figsize=(2*6.4, 1*4.8))
p0 = ax[0].imshow(pupil_single[206:306, 206:306], origin='lower')
plt.colorbar(p0, ax=ax[0])
ax[0].set_title('Pupil')
p1 = ax[1].imshow(np.log10(psf_single[236:276, 236:276]), origin='lower')
plt.colorbar(p1, ax=ax[1])
ax[1].set_title('Point-spread function')
plt.show()
ramp2 = np.arange(512)-512/2.-28.
xx, yy = np.meshgrid(ramp2, ramp)
dist1 = np.sqrt(xx**2+yy**2)
ramp2 = np.arange(512)-512/2.+28.
xx, yy = np.meshgrid(ramp2, ramp)
dist2 = np.sqrt(xx**2+yy**2)

# Compute two aperture interferometer...
pupil_double = (dist1 < 4) | (dist2 < 4)
psf_double = np.abs(np.fft.fftshift(np.fft.fft2(np.fft.fftshift(pupil_double))))**2

# Plotting...
f, ax = plt.subplots(1, 2, figsize=(2*6.4, 1*4.8))
p0 = ax[0].imshow(pupil_double[206:306, 206:306], origin='lower')
plt.colorbar(p0, ax=ax[0])
ax[0].set_title('Pupil')
p1 = ax[1].imshow(np.log10(psf_double[236:276, 236:276]), origin='lower')
plt.colorbar(p1, ax=ax[1])
ax[1].set_title('Point-spread function')
plt.show()<jupyter_output><empty_output><jupyter_text>The peak to zero distance for the Airy pattern is 1.2 lambda/D, whereas it is lambda/2B for an interferometer!<jupyter_code>pupil_temp = np.zeros(pupil_single.shape)
pupil_temp[:, :256] = pupil_single[:, :256]/np.max(pupil_single[:, :256])
pupil_temp[:, 256:] = pupil_double[:, 256:]/np.max(pupil_double[:, 256:])
psf_temp = np.zeros(psf_single.shape)
psf_temp[:, :256] = psf_single[:, :256]
psf_temp[:, 256:] = psf_double[:, 256:]

# Plotting...
f, ax = plt.subplots(1, 2, figsize=(2*6.4, 1*4.8))
p0 = ax[0].imshow(pupil_temp[206:306, 206:306], origin='lower')
plt.colorbar(p0, ax=ax[0])
ax[0].set_title('Pupil')
p1 = ax[1].imshow(np.log10(psf_temp[236:276, 236:276]), origin='lower')
plt.colorbar(p1, ax=ax[1])
ax[1].set_title('Point-spread function')
plt.show()<jupyter_output><empty_output><jupyter_text>Kernel-phase interferometry (KPI)Next, we will look at some kernel-phase data and investigate the properties of kernel-phase. We have kernel-phase data from two different objects, one is a PSF reference star and one is a binary star. Which one is which?As we have learnt during the lecture, the kernel-phase of a point source should theoretically be zero. Instead, the kernel-phase of a binary star should be non-zero. The easiest way to find out which one is which is to make a histogram of the kernel-phase data and to check which one is closer to zero:<jupyter_code># Opening data from a FITS file...
hdul1 = pyfits.open('n0258.fits')
hdul2 = pyfits.open('n0438.fits')

# The kernel-phase data is stored in the KP-DATA extension...
kp1 = hdul1['KP-DATA'].data[0, 0]
kp2 = hdul2['KP-DATA'].data[0, 0]

# Plotting...
f = plt.figure()
ax = plt.gca()
ax.hist(kp1, bins=20, range=(-3, 3), alpha=0.5, label='first')
ax.hist(kp2, bins=20, range=(-3, 3), alpha=0.5, label='second')
ax.set_xlabel('Kernel-phase')
ax.set_ylabel('Number')
ax.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>It is quite obvious that the second one is the PSF reference star while the first one is the binary star. Let's compare their Fourier phase data next:<jupyter_code># The complex visibility data is stored in the CVIS-DATA extension...
cvis1 = hdul1['CVIS-DATA'].data[0, 0, 0]+1.0j*hdul1['CVIS-DATA'].data[1, 0, 0]
cvis2 = hdul2['CVIS-DATA'].data[0, 0, 0]+1.0j*hdul2['CVIS-DATA'].data[1, 0, 0]

# The Fourier phase is simply the phase of the complex visibility...
fp1 = np.angle(cvis1)
fp2 = np.angle(cvis2)

# Plotting...
f = plt.figure()
ax = plt.gca()
ax.hist(fp1, bins=20, range=(-1.5, 1.5), alpha=0.5, label='first')
ax.hist(fp2, bins=20, range=(-1.5, 1.5), alpha=0.5, label='second')
ax.set_xlabel('Fourier phase')
ax.set_ylabel('Number')
ax.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Things are not that obvious anymore because the Fourier phase is affected by pupil plane phase errors, only the kernel-phase is more robust with respect to these errors. This is the advantage of using kernel-phase (or closure-phase)!  PROBLEMS Problem 1In the beginning of this notebook, we considered the object intensity distribution of a single star, which is a Dirac delta function. Let's know consider a single star which is 1 pixel off-center. First, create its object intensity distribution: Next, compute and plot the Fourier transform of its object intensity distribution: Explain what you can see!HINT: page 1 of this document (https://ethz.ch/content/dam/ethz/special-interest/baug/ibk/structural-mechanics-dam/education/identmeth/fourier.pdf) will help.What does this mean for Fourier plane imaging? Do you think centering the data is important before extracting the Fourier observables? Problem 2We have seen that in the high-contrast regime, the Fourier phase of a binary (i.e., two Dirac delta functions) is a sine wave with amplitude roughly equal to the contrast. Starting with the expression for the complex visibility of a binary, which is$V_\rm{bin}(u,v) = 1+c*\exp(-2*\pi*i*(a*u/\lambda+b*v/\lambda))$,where $(u,v)$ is the spatial frequency vector and $(a,b)$ is the separation vector and $\lambda$ is the observing wavelength, proof that the amplitude of the Fourier phase is roughly $c$ for $c \ll 1$.HINT: as you will see, we don't really care about the expression in the $\exp()$ function here. Simply substitute it with $x$ for better readability. Note that $\arctan(z) \approx z$ for $z \ll 1$. Problem 3The ultimate goal of the third problem is to show that (small) random errors in the pupil plane phase do not affect the kernel-phase observables. Let's start with the matrix A, which we assume to be given from our kernel-phase data set:<jupyter_code># The baseline mapping matrix A is stored in the BLM-MAT extension...
A = hdul1['BLM-MAT'].data[:, 1:] # ignore 1. row because we can normalize all phase offsets relative to 1. aperture<jupyter_output><empty_output><jupyter_text>Let's check out its shape:<jupyter_code>print(A.shape)<jupyter_output><empty_output><jupyter_text>This means that our pupil has 104 independent apertures that map onto 210 unique baselines in the Fourier plane! The apertures actually look like this:<jupyter_code>xxc = hdul1['APERTURE'].data['XXC']
yyc = hdul1['APERTURE'].data['YYC']

f = plt.figure()
ax = plt.gca()
ax.scatter(xxc, yyc)
ax.set_aspect('equal')
ax.set_xlabel('Size [m]')
ax.set_ylabel('Size [m]')
plt.show()<jupyter_output><empty_output><jupyter_text>Now compute the kernel of A using singular value decomposition and look at the singular values: There are 104 singular values, none of which is zero. However, there are 210 unique baselines, so what about the remaining 210-104 = 106 kernels? Note that since the matrix A is not square (it has more columns than rows), it has a trivial kernel. This trivial kernel are all the 106 remaining basis vectors. Hence, we append 106 zeros to the list of singular values: Now, compute the kernel of A similar as before: Run the test that K.A = 0: Now that we have the matrix K, the only thing that is missing from the kernel-phase equation is the diagonal redundancy matrix R. Luckily, it is also stored in the data:<jupyter_code># The redundancy matrix R is stored in the UV-PLANE extension
R = np.diag(hdul1['UV-PLANE'].data['RED'])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12831089882761615, 'fraction_numerical': 0.038355767839050516, 'mean_word_length': 5.318701417466849, 'pattern_counts': {'":': 0, '<': 74, '<?xml version=': 0, '>': 70, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '347597', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5121, 'n_tokens_neox': 4689, 'n_words': 1629, 'path': 'Class10/ami_kpi.ipynb', 'repo_name': 'anand0xff/FourierOpticsJHU2021', 'size': '13825'}