<jupyter_start><jupyter_text>Fourier transforms of simple 2D functionsSo, in our application of optics, we will not be dealing with one dimensional plots but with two dimensional images instead. This means that our function $f$ in direct space is a function of **two** independent spatial variables, $f = f(x, y)$, and the FT of that will also have two independent spatial frequency variables, $F = F(s, r)$.The two-dimensional Fourier transform does exactly the same like the one-dimensional FT, decompose a signal in its sinusoidal components except that we now have to think all across a 2D surface as opposed to a single, 1D line.We have not touched yet any of the important basic theorems, but I though it is really important to get the visualizations for 2D transforms in, especially since this will set us up for jumping into Fourier Optics and from then on work with the optical application on our minds.The analytical definition of the 2D FT is:$$F(s, r) = \mathscr{F}\{f(x, y)\}(s, r) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x, y) e^{-i 2 \pi (x s+yr)} dx dy$$There are many functions $f(x,y)$ that have a well-defined analytical Fourier transform, but I will refrain from showing them here and focus only on the numerical implementation.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D   # this is for surface plots
from matplotlib import cm
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Like for the 1D case, we will define some convenience functions that will perform a 2D Fourier Transform with the respective array shifts needed. And I define some utility functions.<jupyter_code>def ft2d(func):
    ft = np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(func)))
    return ft

def ift2d(func):
    ift = np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(func)))
    return ift

def ft1d_freq(x):
    """Calculate the (spatial) frequency array based on the spatial array x."""
    s = np.fft.fftshift(np.fft.fftfreq(x.size, d=x[-1]-x[-2]))
    return s

def zoom(im, x, y, bb):
    """Cut out a square box from image im centered on (x,y) with half-box size bb."""
    return(im[y-bb:y+bb,x-bb:x+bb])

def padcplx(c, pad=5):
    """Puts a Complex array in the centre of a zero-filled Complex array.
    pad defines the padding multiplier for the output array."""
    (nx, ny) = c.shape
    bignx = nx * pad + 1
    bigny = ny * pad + 1
    big_c = np.zeros((bignx, bigny), dtype=complex)
    
    dx = int((nx * (pad-1)) / 2 + 1)
    dy = int((ny * (pad-1)) / 2 + 1)
    
    big_c[dx:dx+nx,dy:dy+ny] = c
    return(big_c)<jupyter_output><empty_output><jupyter_text>Setting up the gridsWe will first generate a grid for the independent variables $x$ and $y$. I chose the grid to have unity size and go from -0.5 to 0.5 so that I only ever have to define the sizes of my functions $f(x,y)$ with respect to the total grid.<jupyter_code># We first need to generate the independent variables,
# this time as a 2D grid.
npix = 512

lin = np.linspace(-0.5, 0.5, npix)
xx, yy = np.meshgrid(lin, lin)

print("Shape of xx: {}".format(xx.shape))
print("Shape of yy: {}".format(yy.shape))
print("xx[:3,0]: {}".format(xx[:3,0]))
print("xx[-3:,0]: {}".format(xx[-3:,0]))

plt.figure(figsize=(10,10))
plt.subplot(1, 2, 1)
plt.imshow(xx)
plt.title('xx')
plt.colorbar()
plt.subplot(1, 2, 2)
plt.imshow(yy)
plt.title('yy')
plt.colorbar()<jupyter_output><empty_output><jupyter_text>We will also set up the spatial frequency arrays for the variables $s$ and $r$.  Now, one thing that I don't want to get into in this notebook just yet is the subject of *padding* and *sampling*. I was able to get away with ignoring these things in the 1D cases because I just picked my different $f(x)$ in such a way that we were never undersampled and I never had to tweak the plots much in order to get nice displays. This gets more complicated in the 2D case and while I could still be smart about it in the same way, the plotting will need major adjustments, so I can just as well add the padding now to make my life easier.Again, I will come back to this at a later point, the short version is that we need to shove our $f(x,y)$ images into bigger arrays that pad the function before performing a Fourier transform. This means that the data points we are interested in only take up a small fraction of the arrays we're dealing with, and I will need to either adjust the plot ranges or specifically cut out the data parts that we want to inspect.<jupyter_code>pad = 5   # factor by how much do we pad our images before performing a FT

npix_pad = npix*pad+1   # figure out the padded big array sizes after the FT
xsf = np.linspace(-0.5, 0.5, npix_pad)
gf = ft1d_freq(xsf)
ss, rr = np.meshgrid(gf, gf)   # create the spatial frequency grids

print("Shape of ss: {}".format(ss.shape))
print("Shape of rr: {}".format(rr.shape))

plt.figure(figsize=(10,10))
plt.subplot(1, 2, 1)
plt.imshow(xx)
plt.title('ss')
plt.subplot(1, 2, 2)
plt.imshow(yy)
plt.title('rr')<jupyter_output><empty_output><jupyter_text>A 2D rectangle - a rectangular aperture Numerical representation of the rectangle functionWe can represent a rectangular aperture with a 2D function where the function value at each point on the grid is the amplitude of the aperture at that point. In the application of optics, simple shapes like squares, rectangles or circles are usually the functions used to describe pupil apertures, so the amplitude of a function like that is equivalent to the transmission of the aperture. In the simplest case, the aperture function will have non-zero values inside the aperture where it is transmissive to the (optical) signal and zeros outside of it, where we have no signal. Since the two extreme cases here are "all signal" vs. "no signal", we can represent that with zeros and ones. If we have an aperture that is partially transmissive in some areas, it will contain values between 0 and 1 there. In practical terms, this means that our $A$ will always have values between 0 and 1.<jupyter_code># Define a function for the rectangular aperture.
# The expression below returns values of True or
# False, so we have to convert that into floats.
def rect2d_noshift(size):
    """Rectangluar aperture. size is a tuple (x,y)."""
    rect = (np.abs(xx) <= (size[0]/2)) * (np.abs(yy) <= (size[1]/2))
    return rect.astype('float')
# Create a rectangular aperture
T = (0.4, 0.2)   # x and y size of the rectangular aperture
rect_ap = rect2d_noshift(T)

# Display the full function
fig = plt.figure(figsize=(15,7))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(xx, yy, rect_ap)
ax.set_title('$f(x, y)$ - rectangular aperture')
ax.set_xlabel('x')
ax.set_ylabel('y')<jupyter_output><empty_output><jupyter_text>Oftentimes it will not give us any extra information to display a 2D function like this, but instead the amplitude of the function at each point will be indicated by different colors in an instance of `plt.imshow()`.<jupyter_code># Display the rectangular aperture wtih imshow()
plt.figure(figsize=(5, 5))
plt.imshow(rect_ap)
plt.colorbar()<jupyter_output><empty_output><jupyter_text>We can see how the the transmissive part of this aperture is filled with ones, while the opaque parts around it are filled with zeros. FT of a rectangular apertureI will show the functions in the following both in a 3D surface plot as well as on a 2D image, just to give a feeling for what's going on. Also, remember that `plt.imshow()` cannot plot complex numbers, and as opposed to `plt.plot()` it will *not* default to plotting the real part, but will raise an error instead. We have to speficy what we want it to show us.<jupyter_code># Calculate the FT
rec_ft = ft2d(padcplx(rect_ap))   # we're zero-padding f(x,y) here, see note above

# Plot
fig = plt.figure(figsize=(15,15))
plt.suptitle('$F(s,r)$')

ax1 = fig.add_subplot(221, projection='3d')
ax1.plot_surface(ss, rr, np.real(rec_ft), cmap=cm.coolwarm)
ax1.set_title('Real')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.view_init(10, -60)
ax1.set_zlim3d(-5000,5000)

ax2 = fig.add_subplot(222, projection='3d')
ax2.plot_surface(ss, rr, np.imag(rec_ft), cmap=cm.coolwarm)
ax2.set_title('Imaginary')
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.view_init(30, -60)

ax3 = fig.add_subplot(223)
im1 = ax3.imshow(np.real(rec_ft))
ax3.set_title('Real')
fig.colorbar(im1)

ax4 = fig.add_subplot(224)
im2 = ax4.imshow(np.imag(rec_ft))
ax4.set_title('Imaginary')
fig.colorbar(im1)<jupyter_output><empty_output><jupyter_text>This might not look like much at first, but everything is ok here. All we have to do is zoom in a little in order to see the parts that are interesting to us (see also my note in the first part of the notebook where I create the grids).  Also, we can see that the imaginary part of this FT is just noise and not interesting to us right now.<jupyter_code># Zoom into the FT of the rectangular aperture
zoomfac = 20     # half-size of the zoom box will be 1/zoomfac of total image
box = int(npix_pad/zoomfac)

# This is a smaller data array wiht our region of interest.
rec_ft_zoom = zoom(rec_ft, int(npix_pad/2), int(npix_pad/2), box)

# I also have to adjust our s and r grids to match the zoomed data size
sz = zoom(ss, int(npix_pad/2), int(npix_pad/2), box)
rz = zoom(rr, int(npix_pad/2), int(npix_pad/2), box)

# Plot
fig = plt.figure(figsize=(20,7))
ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(sz, rz, np.real(rec_ft_zoom), cmap=cm.coolwarm)
ax1.set_title('$F(s,r)$')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.view_init(10, -60)
ax1.set_zlim3d(-5000,20000)

ax3 = fig.add_subplot(122)
im1 = ax3.imshow(np.real(rec_ft_zoom))
ax3.set_title('$F(s,r)$')
fig.colorbar(im1)<jupyter_output><empty_output><jupyter_text>Nice! We can see very well in the right hand plot what the FT of a rectangular aperture looks like. Note how the FT seems to be rotated with respect to the rectangle we transformed: while the rectangle has a longer extention along the x-axis, its FT is elongated along the y-axis. This is one of the main properties of the Fourier transform, where extended features in the direct plane become narrower, and narrow features get extended in the Fourier plane. We will cover this in more detail later.Let's move on to different 2D functions $f(x,y)$ that we want to transform. Circular apertureThe circular aperture is a classic one, since the standard case for a telescope is a circular entrance pupil.<jupyter_code># Define a function for a circular aperture
def circle_mask(im, xc, yc, rcirc):
    """Create a circular aperture centered on (xc, yc) with radius rcirc."""
    x, y = np.shape(im)
    newy, newx = np.mgrid[:y,:x]
    circ = (newx-xc)**2 + (newy-yc)**2 < rcirc**2
    return circ.astype('float')
# Create a circular aperture
rad = 0.7 * npix/2   # radius in pixels of the circular aperture
circ_ap = circle_mask(xx, int(npix/2), int(npix/2), rad)

# Plot
fig = plt.figure(figsize=(15,7))
plt.suptitle('$f(x, y)$ - circular aperture')

ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(xx, yy, circ_ap)
ax1.set_xlabel('x')
ax1.set_ylabel('y')

ax2 = fig.add_subplot(122)
im1 = ax2.imshow(circ_ap)
#fig.colorbar(im1)<jupyter_output><empty_output><jupyter_text>FT of a circular aperture<jupyter_code># Calculate the FT
circ_ft = ft2d(padcplx(circ_ap))
# Plot
zoomfac = 30     # half-size of the zoom box will be 1/zoomfac of total image
box = int(npix_pad/zoomfac)

# This is a smaller data array wiht our region of interest.
circ_ft_zoom = zoom(circ_ft, int(npix_pad/2), int(npix_pad/2), box)

# I also have to adjust our s and r grids to match the zoomed data size
sz = zoom(ss, int(npix_pad/2), int(npix_pad/2), box)
rz = zoom(rr, int(npix_pad/2), int(npix_pad/2), box)

# Plot
fig = plt.figure(figsize=(20,7))
ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(sz, rz, np.real(circ_ft_zoom), cmap=cm.coolwarm)
ax1.set_title('$F(s,r)$')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.view_init(10, -60)

ax3 = fig.add_subplot(122)
im1 = ax3.imshow(np.real(circ_ft_zoom))
ax3.set_title('$F(s,r)$')
fig.colorbar(im1)<jupyter_output><empty_output><jupyter_text>The FT of a circle yields a circular point in the middle of the image and additional concentric rings around it. It is such an important function that it has its own name: the **Airy funciton**. The circular core is also called the Airy core and the rings around it the Airy rings.  Inverse FTLet's see whether doing the inverse Fourier transform of the Fourier transform yields what we'd expect, like we did in the 1D case.<jupyter_code># Take the inverse FT
circ_ft_back = ift2d(circ_ft)
# Plot
circ_ft_back_zoom = zoom(circ_ft_back, int(npix_pad/2), int(npix_pad/2), int(npix/2))

fig = plt.figure(figsize=(15,7))
plt.suptitle('$f(x, y)$ - circular aperture')

ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(xx, yy, np.real(circ_ft_back_zoom))
ax1.set_xlabel('x')
ax1.set_ylabel('y')

ax2 = fig.add_subplot(122)
ax2.imshow(np.real(circ_ft_back_zoom))<jupyter_output><empty_output><jupyter_text>Looking pretty good, but how good is it really? Lets subtract the back-transformed aperture from the original one and look at the residuals.<jupyter_code>test = circ_ap - circ_ft_back_zoom

plt.figure(figsize=(20, 20))
plt.imshow(np.real(test))
plt.title('Back-transformed residual')
plt.colorbar()<jupyter_output><empty_output><jupyter_text>This is good enough to confirm that the invers FT of the FT gave us the original function $f(x,y)$, but we can also see that there is a small offset in the back-transformed function compared to the original one. We will keep this in mind and come back to it at a later point. Gaussian A 2D Gaussian is a possibility for creating tapered apertures, among other things. It can almost be viewed as a circular aperture with transmission values changing *gradually* from 0 to 1, with a Gaussian profile, instead of the hard jump from 0 to 1 in the standard circular aperture.<jupyter_code># Define a 2D Gaussian function
def gaussian2d(x, y, ampl, c):
    #func = ampl * np.exp(-np.square(x) / (2*np.square(c)))
    func = ampl * np.exp(-np.pi * (np.square(x) + np.square(y)) / (2 * c/2))
    return func
A = 1
T = 1
gauss = gaussian2d(xx, yy, A, T)

# Plot
fig = plt.figure(figsize=(15,7))
plt.suptitle('$f(x, y)$ - 2D Gaussian')

ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(xx, yy, gauss, cmap=cm.coolwarm)
ax1.set_xlabel('x')
ax1.set_ylabel('y')

ax2 = fig.add_subplot(122)
im1 = ax2.imshow(gauss)
fig.colorbar(im1)<jupyter_output><empty_output><jupyter_text>FT of a 2D Gaussian<jupyter_code># Perform the FT
gauss_ft = ft2d(padcplx(gauss))
# Plot
zoomfac = 30     # half-size of the zoom box will be 1/zoomfac of total image
box = int(npix_pad/zoomfac)

# This is a smaller data array wiht our region of interest.
gauss_ft_zoom = zoom(gauss_ft, int(npix_pad/2), int(npix_pad/2), box)

# The sz and rz arrays for this already exist.

# Plot
fig = plt.figure(figsize=(20,7))
ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(sz, rz, np.real(gauss_ft_zoom), cmap=cm.coolwarm)
ax1.set_title('$F(s,r)$')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.view_init(10, -60)

ax3 = fig.add_subplot(122)
im1 = ax3.imshow(np.real(gauss_ft_zoom))
ax3.set_title('$F(s,r)$')
fig.colorbar(im1)<jupyter_output><empty_output><jupyter_text>2D sines and cosinesSo far, I've been showing Fourier transforms of very simple functions $f(x,y)$. We are by no means restricted to work with binary functions like the rectangular or circular aperture. There's a particular set of functions that are of special intersest, as they have been in the 1D case - sinusoids.  Since we have two dimensions now, we could go ahead and make functions that have a sinusoid extending both ways, on the x- and on the y-axis. To keep things simple though, we will consider 2D sinusoid in only one dimension (but as a 2D function), while we keep the other dimension constant.From notebook 2, we remember that a sinusoidal function can generally be expressed as:$$f(x) = A\ cos(2\pi \nu x - \phi)$$Extendings this to two dimensions yields$$f(x,y) = A\ cos(2\pi \nu (x+y) - \phi)$$where $\nu$ is the frequency of $f(x,y)$: the number of complete oscillations the graph makes in a spatial interval of length 1; that is, the number of cycles per unit length, and $\phi$ is the phaselag: the value of $\omega x$ by which the function is shifted along $x$ (or $y$). The variables $x$ and $y$ will be represented by our grids `xx` and `yy`, respectively. Remember that a $sine$ and $cosine$ wave only differ by a phaselag of $\phi = 90^{\circ}$.<jupyter_code># Define a 2D sinusoid function
def sinusoid2d_norot(x, y, nu, phi, A):
    func = A * np.cos(2 * np.pi * nu * (x+y) - phi)
    return func
A = 1
phi = 0   # for phi = 0, we're doing a pure cosine
nu = 6

sinus_norot = sinusoid2d_norot(xx, yy, nu, phi, A)

plt.imshow(np.real(sinus_norot))<jupyter_output><empty_output><jupyter_text>This is looking really nice, but I would like to have this sinusoidal ripple to be along one axis only. I can achieve this by adding an axes rotation to my function, by the angle $\theta$.<jupyter_code>def sinusoid2d(x, y, nu, phi, A, theta=0):
    xr = x * np.cos(theta)
    yr = y * np.sin(theta)
    
    func = A * np.cos(2 * np.pi * nu * (xr+yr) - phi)
    return func
A = 1
phi = 0   # for phi = 0, we're doing a pure cosine
nu = 6

sinus = sinusoid2d(xx, yy, nu, phi, A)

# Plot
fig = plt.figure(figsize=(15,7))
plt.suptitle('$f(x, y)$ - 2D sinusoid')

ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(xx, yy, sinus, cmap=cm.coolwarm)
ax1.set_xlabel('x')
ax1.set_ylabel('y')

ax2 = fig.add_subplot(122)
im1 = ax2.imshow(sinus)
fig.colorbar(im1)<jupyter_output><empty_output><jupyter_text>This gives us a beautiful sinusoidal ripple that defaults to being along the x-direction with amplitude $A$. What does the Fourier transform of something like this look like? FT of a 2D sinusoid<jupyter_code># Perform the FT
sinus_ft = ft2d(padcplx(sinus))
# Plot
zoomfac = 30     # half-size of the zoom box will be 1/zoomfac of total image
box = int(npix_pad/zoomfac)

# This is a smaller data array wiht our region of interest.
sinus_ft_zoom = zoom(sinus_ft, int(npix_pad/2), int(npix_pad/2), box)

# The sz and rz arrays for this already exist.

# Plot
fig = plt.figure(figsize=(20,7))
ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(sz, rz, np.real(sinus_ft_zoom), cmap=cm.coolwarm)
ax1.set_title('$F(s,r)$')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.view_init(10, -60)

ax3 = fig.add_subplot(122)
im1 = ax3.imshow(np.real(sinus_ft_zoom))
ax3.set_title('$F(s,r)$')
fig.colorbar(im1)<jupyter_output><empty_output><jupyter_text>This is analogous to the 1D case! Remember fro notebook 4 how the transform of a 1D sinusoid yielded two delta functions - the 2D case of this is just the same!  Two circles We have gone through the most basic 2D functions, took their FT and displayed them. Going forward, we can make any sort of functions $f(x,y)$ and look at its FT $F(s,r)$.Looking at the FT of a sinusoidal function above, I can see that it yields someting that looks like two dots. Remember how we called a function $f(x,y)$ and its FT $F(s,r)$ a **Fourier pair**. We also know that we can go back and forth between the two members of a Fourier pair by doing a FT and its inverse.So, since the FT of the sinusoidal looks so much like two dots in the Fourier plane, how about we check what happens when we have two dots in real space instead, and then transform *that* into Fourier space? We can do this by creating two off-centered circular apertures and then combining them into one function.<jupyter_code># Make two circular apertures with different offset from center
rad_dots = 20
circ1 = circle_mask(xx, int(npix/2)-100, int(npix/2), rad_dots)
circ2 = circle_mask(xx, int(npix/2)+100, int(npix/2), rad_dots)

plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.imshow(circ1)
plt.title('circ1')
plt.subplot(122)
plt.imshow(circ2)
plt.title('circ2')
# Combine them into one array - this will be our f(x,y)
dots = circ1 + circ2

# Plot
fig = plt.figure(figsize=(15,7))
plt.suptitle('$f(x, y)$ - dots')

ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(xx, yy, dots)
ax1.set_xlabel('x')
ax1.set_ylabel('y')

ax2 = fig.add_subplot(122)
im1 = ax2.imshow(dots)<jupyter_output><empty_output><jupyter_text>FT of two dots<jupyter_code># Perform the FT
dots_ft = ft2d(padcplx(dots))
# Plot
zoomfac = 30     # half-size of the zoom box will be 1/zoomfac of total image
box = int(npix_pad/zoomfac)

# This is a smaller data array wiht our region of interest.
dots_ft_zoom = zoom(dots_ft, int(npix_pad/2), int(npix_pad/2), box)

# The sz and rz arrays for this already exist.

# Plot
fig = plt.figure(figsize=(20,7))
ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(sz, rz, np.real(dots_ft_zoom), cmap=cm.coolwarm)
ax1.set_title('$F(s,r)$')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.view_init(10, -60)

ax3 = fig.add_subplot(122)
im1 = ax3.imshow(np.real(dots_ft_zoom))
ax3.set_title('$F(s,r)$')
fig.colorbar(im1)<jupyter_output><empty_output><jupyter_text>This looks like a somewhat modulated sinusoid. That's because it is! We said that we can come back to the original function if we apply an inverse FT to its FT. This means that between a sinusoid and two dots, which are essentially broadened Delta functions, only ever lies a Fourier transform, no matter which of the two functions (sinusoid or dots) we define as our original funciton in real space! If the FT of a sinusoid are two dots, then the FT of two dots will be a sinusoid. What about the modulation though? Since we're not actually dealing with two Delta functions here, but with two circles with hard edges, the sinusoid in the Fourier plane will end up fuzzed out. This behavior of hard edges is yet another property of the Fourier transfor that plays a major role in astronomical optics.  Slits - two rectangles We can do the same thing with the dots but using two rectangles instead, effentively creating two slits for our function $f(x,y)$ that we then transform.<jupyter_code># Define a rectangle function in which you can define the shift of the rectangles
def rect2d(x, y, size, shift):
    """Rectangluar aperture. size and cen are a tuple (x,y)."""
    rect = (np.abs(x-shift[0]) <= (size[0]/2)) * (np.abs(y-shift[1]) <= (size[1]/2))
    return rect.astype('float')
# Create the two slits
# size and shift in terms of fraction of total array
T = (0.1, 0.3)

shift1 = (0.25, 0)
slit1 = rect2d(xx, yy, T, shift1)
shift2 = (-0.25, 0)
slit2 = rect2d(xx, yy, T, shift2)

plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.imshow(slit1)
plt.title('slit1')
plt.subplot(122)
plt.imshow(slit2)
plt.title('slit2')
# Combine them into one array - this will be our f(x,y)
slits = slit1 + slit2

# Plot
fig = plt.figure(figsize=(15,7))
plt.suptitle('$f(x, y)$ - slits')

ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(xx, yy, slits)
ax1.set_xlabel('x')
ax1.set_ylabel('y')

ax2 = fig.add_subplot(122)
im1 = ax2.imshow(slits)<jupyter_output><empty_output><jupyter_text>FT of slits<jupyter_code># Perform the FT
slits_ft = ft2d(padcplx(slits))
# Plot
zoomfac = 30     # half-size of the zoom box will be 1/zoomfac of total image
box = int(npix_pad/zoomfac)

# This is a smaller data array wiht our region of interest.
slits_ft_zoom = zoom(slits_ft, int(npix_pad/2), int(npix_pad/2), box)

# The sz and rz arrays for this already exist.

# Plot
fig = plt.figure(figsize=(20,7))
ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(sz, rz, np.real(slits_ft_zoom), cmap=cm.coolwarm)
ax1.set_title('$F(s,r)$')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.view_init(10, -60)

ax3 = fig.add_subplot(122)
im1 = ax3.imshow(np.real(slits_ft_zoom))
ax3.set_title('$F(s,r)$')
fig.colorbar(im1)<jupyter_output><empty_output><jupyter_text>RingFor the fun of it, I will also do the Fourier transform of a ring - because why not.<jupyter_code># A ring is essentially a big circle minus a smaller circle in the middle
rad_out = 0.7 * npix/2   # radius in pixels of the circular aperture
rad_in = 0.55 * npix/2
circ_out = circle_mask(xx, int(npix/2), int(npix/2), rad_out)
circ_in = circle_mask(xx, int(npix/2), int(npix/2), rad_in)

plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.imshow(circ_out)
plt.title('circ_out')
plt.subplot(122)
plt.imshow(circ_in)
plt.title('circ_in')
# Subtract the smaller circle off the bigger one to create ring
ring = circ_out - circ_in

# Plot
fig = plt.figure(figsize=(15,7))
plt.suptitle('$f(x, y)$ - ring')

ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(xx, yy, ring)
ax1.set_xlabel('x')
ax1.set_ylabel('y')

ax2 = fig.add_subplot(122)
im1 = ax2.imshow(ring)<jupyter_output><empty_output><jupyter_text>FT of a ring<jupyter_code># Perform the FT
ring_ft = ft2d(padcplx(ring))
# Plot
zoomfac = 30     # half-size of the zoom box will be 1/zoomfac of total image
box = int(npix_pad/zoomfac)

# This is a smaller data array wiht our region of interest.
ring_ft_zoom = zoom(ring_ft, int(npix_pad/2), int(npix_pad/2), box)

# The ss and rr arrays for this already exist.

# Plot
fig = plt.figure(figsize=(20,7))
ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(sz, rz, np.real(ring_ft_zoom), cmap=cm.coolwarm)
ax1.set_title('$F(s,r)$')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.view_init(10, -60)

ax3 = fig.add_subplot(122)
im1 = ax3.imshow(np.real(ring_ft_zoom))
ax3.set_title('$F(s,r)$')
fig.colorbar(im1)<jupyter_output><empty_output><jupyter_text>Comparing the different FTsWe have seen quite a number of examples for 2D functions and their Fourier transforms in this notebook. I will display the Fourier pairs we covered below to have a look at them all at once.Note how in all the figures of the FTs so far, we have been plotting their real parts, as in most cases, the imaginary parts are low-level noise. Anticipating what's coming in the next notebook, notebook 6 introducing light as an electric field, I will also plot the **absolute value squared** of the Fourier transforms, since that is a quantity we're very interested in - it's the **intensity** of the signal.<jupyter_code># Plot all 2D functions, the real parts of their FTs and their intensities
plt.figure(figsize=(15, 35))

plt.subplot(7, 3, 1)
plt.imshow(rect_ap)
plt.title('Rectangular aperture')

plt.subplot(7, 3, 2)
plt.imshow(np.real(rec_ft_zoom))
plt.title('Real part')

plt.subplot(7, 3, 3)
plt.imshow(np.abs(rec_ft_zoom)**2)
plt.title('Absolute value squared - intensity')

plt.subplot(7, 3, 4)
plt.imshow(circ_ap)
plt.title('Circular aperture')

plt.subplot(7, 3, 5)
plt.imshow(np.real(circ_ft_zoom))
plt.title('Real part')

plt.subplot(7, 3, 6)
plt.imshow(np.abs(circ_ft_zoom)**2)
plt.title('Absolute value squared - intensity')

plt.subplot(7, 3, 7)
plt.imshow(gauss)
plt.title('Gaussian')

plt.subplot(7, 3, 8)
plt.imshow(np.real(gauss_ft_zoom))
plt.title('Real part')

plt.subplot(7, 3, 9)
plt.imshow(np.abs(gauss_ft_zoom)**2)
plt.title('Absolute value squared - intensity')

plt.subplot(7, 3, 10)
plt.imshow(sinus)
plt.title('Sinusoid')

plt.subplot(7, 3, 11)
plt.imshow(np.real(sinus_ft_zoom))
plt.title('Real part')

plt.subplot(7, 3, 12)
plt.imshow(np.abs(sinus_ft_zoom)**2)
plt.title('Absolute value squared - intensity')

plt.subplot(7, 3, 13)
plt.imshow(dots)
plt.title('Dots')

plt.subplot(7, 3, 14)
plt.imshow(np.real(dots_ft_zoom))
plt.title('Real part')

plt.subplot(7, 3, 15)
plt.imshow(np.abs(dots_ft_zoom)**2)
plt.title('Absolute value squared - intensity')

plt.subplot(7, 3, 16)
plt.imshow(slits)
plt.title('Slits')

plt.subplot(7, 3, 17)
plt.imshow(np.real(slits_ft_zoom))
plt.title('Real part')

plt.subplot(7, 3, 18)
plt.imshow(np.abs(slits_ft_zoom)**2)
plt.title('Absolute value squared - intensity')

plt.subplot(7, 3, 19)
plt.imshow(ring)
plt.title('Ring')

plt.subplot(7, 3, 20)
plt.imshow(np.real(ring_ft_zoom))
plt.title('Real part')

plt.subplot(7, 3, 21)
plt.imshow(np.abs(ring_ft_zoom)**2)
plt.title('Absolute value squared - intensity')<jupyter_output><empty_output><jupyter_text>Usually, we would display the intensity on a logarithmic scale to be able to distinguish features with lower intensity, meaning the parts further away from the center, as well as high intensity features in the very center. And while not all of these functions are azimuthally symmetric, it can be helpful to plot an intensity profile, which we will also show on a log scale for the sake of the faint features.<jupyter_code># We will do the radial cut through the center of the images
# on the y-axis, having a radial plot along the x-axis.
cen = int(rec_ft_zoom.shape[0]/2)

# Plot all 2D functions, the log of their intensitiese and a radial cut of their intensities
plt.figure(figsize=(15, 35))

plt.subplot(7, 3, 1)
plt.imshow(rect_ap)
plt.title('Rectangular aperture')

plt.subplot(7, 3, 2)
plt.imshow(np.log10(np.abs(rec_ft_zoom)**2))
plt.title('Log intensity')

plt.subplot(7, 3, 3)
plt.plot(np.log10(np.abs(rec_ft_zoom)**2)[cen, :])
plt.title('Radial cut log intensity')

plt.subplot(7, 3, 4)
plt.imshow(circ_ap)
plt.title('Circular aperture')

plt.subplot(7, 3, 5)
plt.imshow(np.log10(np.abs(circ_ft_zoom)**2))
plt.title('Log intensity')

plt.subplot(7, 3, 6)
plt.plot(np.log10(np.abs(circ_ft_zoom)**2)[cen, :])
plt.title('Radial cut log intensity')

plt.subplot(7, 3, 7)
plt.imshow(gauss)
plt.title('Gaussian')

plt.subplot(7, 3, 8)
plt.imshow(np.log10(np.abs(gauss_ft_zoom)**2))
plt.title('Log intensity')

plt.subplot(7, 3, 9)
plt.plot(np.log10(np.abs(gauss_ft_zoom)**2)[cen, :])
plt.title('Radial cut log intensity')

plt.subplot(7, 3, 10)
plt.imshow(sinus)
plt.title('Sinusoid')

plt.subplot(7, 3, 11)
plt.imshow(np.log10(np.abs(sinus_ft_zoom)**2))
plt.title('Log intensity')

plt.subplot(7, 3, 12)
plt.plot(np.log10(np.abs(sinus_ft_zoom)**2)[cen, :])
plt.title('Radial cut log intensity')

plt.subplot(7, 3, 13)
plt.imshow(dots)
plt.title('Dots')

plt.subplot(7, 3, 14)
plt.imshow(np.log10(np.abs(dots_ft_zoom)**2))
plt.title('Log intensity')

plt.subplot(7, 3, 15)
plt.plot(np.log10(np.abs(dots_ft_zoom)**2)[cen, :])
plt.title('Radial cut log intensity')

plt.subplot(7, 3, 16)
plt.imshow(slits)
plt.title('Slits')

plt.subplot(7, 3, 17)
plt.imshow(np.log10(np.abs(slits_ft_zoom)**2))
plt.title('Log intensity')

plt.subplot(7, 3, 18)
plt.plot(np.log10(np.abs(slits_ft_zoom)**2)[cen, :])
plt.title('Radial cut log intensity')

plt.subplot(7, 3, 19)
plt.imshow(ring)
plt.title('Ring')

plt.subplot(7, 3, 20)
plt.imshow(np.log10(np.abs(ring_ft_zoom)**2))
plt.title('Log intensity')

plt.subplot(7, 3, 21)
plt.plot(np.log10(np.abs(ring_ft_zoom)**2)[cen, :])
plt.title('Radial cut log intensity')<jupyter_output><empty_output><jupyter_text>Good job, team, we've seen now quite a number of different Fourier transform pairs. In the next notebook, we will explain the connection of all of this to the application of optics in astronomy and look at a very simple telescope.You might have noticed that I mentioned a lot of Fourier properties on the fly and then said we'll have a closer look at them later. It is true that we still didn't treat a single one of the very basic Fourier theorems, but I really want to introduce the practical aspect of Fourier transforms in form of their application to telescopes, so I will have to defer further (super important!) mathematical aspects of the Fourier transform to afer notebook 6.<jupyter_code>#TODO: random apertures - hm, maybe better to put this in nb6, since that's where we get optics context:
# -> circular with struts and central obscuration, hexagonal, etc, etc, etc, ...

#TODO: more comments on the different apertures and their FT
#TODO: all the 3D plots kinda look like shit, I need to improve them
#TODO: I want to look through this:
#http://www.robots.ox.ac.uk/~az/lectures/ia/lect2.pdf<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1202130993418991, 'fraction_numerical': 0.02795361955499843, 'mean_word_length': 4.97920179876335, 'pattern_counts': {'":': 0, '<': 110, '<?xml version=': 0, '>': 106, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '26', 'dir': 'jupyter-structured-clean-dedup', 'id': '589144', 'license': 'MIT', 'n_tokens_mistral': 12199, 'n_tokens_neox': 11313, 'n_words': 3898, 'path': '05_2D Fourier transforms.ipynb', 'repo_name': 'ivalaginja/FouFourier', 'size': '31910'}