<jupyter_start><jupyter_text>Simulation of plane wave incident on dual microlens arrays (MLA's) for laser beam array generation.This notebook is used to verify that the results of the code match Fig. 7b in [<NAME>, "Wave optical analysis of light-emitting diode beam shaping using microlens arrays", SPIE, 2002](http://opticalengineering.spiedigitallibrary.org/article.aspx?articleid=1099186).<jupyter_code>%pylab
%matplotlib inline
import SimMLA.fftpack as simfft
import SimMLA.grids   as grids<jupyter_output>Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib<jupyter_text>Setup the MLA geometryWe'll start by defining the geometry of the MLA's.+ Prior to the objective there are two MLA's with the same parameters.+ They are spaced by one focal length such that the second MLA is in the focal plane of the first.+ Each lenslet has the same focal length as all the other lenslets. This is true for lenslets in both MLA's.**Note: This notebook contains LaTeX that may not be visible when viewed from GitHub. Try downloading it and opening it with the Jupyter Notebook application.**<jupyter_code>numLenslets = 11    # Must be odd; corresponds to the number of lenslets in one dimension
lensletSize = 200   # microns
focalLength = lensletSize / 2 / 0.11 # microns

wavelength  = 0.520 # microns<jupyter_output><empty_output><jupyter_text>Now we'll establish a GridArray to represent the microlenses. Such an array is a coordinate system onto which the individual lenslets are mapped. The GridArray is primarily a convenience class because it manages all unit conversions when performing Fourier transforms.The most important parameter here is the subgridSize. This is the number of discrete lattice sites into which a lenslet will be divided. The spatial sampling period for each lenslet will therefore be$$\delta = \frac{p}{N_{\text{sub}}}$$where $p$ is the MLA pitch in microns and $ N_{\text{sub}} $ is the subgrid size, i.e. the number of discrete lattice sites along one dimension of the MLA. For example, a MLA pitch of $p = 500 \, \mu m$ and a subgrid size of $N_{\text{sub}} = 51$ will result in a spatial sampling period of $\frac{500 \, \mu m}{51} \approx 9.8 \mu m$.A good sampling rate should be smaller than the fastest oscillations in the amplitude and phase of the Fourier transform of the input field to avoid aliasing. Unfortunately, it is impossible to completely avoid aliasing when computing transforms of real apertures, which are not band-limited. We can however minimize their effects by choosing sampling periods as small as possible (or, equivalently, sampling frequencies as large as possible).We also need to zero pad the input by a large enough amount to minimize rippling effects in the solutions. If we don't zero pad, the numerical Fourier transforms of the input fields are convolved with a sinc whose width is equal to the width of the field's spectrum. If this occurs, the sinc-based interpolation of the input field's spectrum results in additional artifacts. The spatial sampling period used below is$$ \delta = \frac{p}{N_{\text{sub}}} = \frac{200 \, \mu m}{100001} $$<jupyter_code>print('Spatial sampling period: {0:.4f} microns'.format(200 / 100001.0))
subgridSize  = 100001                     # Number of grid (or lattice) sites for a single lenslet
physicalSize = numLenslets * lensletSize # The full extent of the MLA

# dim = 1 makes the grid 1D
# zeroPad = N makes the computational grid N times larger than the MLA
grid = grids.GridArray(numLenslets, subgridSize, physicalSize, wavelength, focalLength, dim = 1, zeroPad = 5)<jupyter_output><empty_output><jupyter_text>Define the input fieldNext, we'll define an input laser beam that will impinge upon the first MLA. Because this is a 1D geometry, the input field \\( u \left( x \right) \\) is in units of \\( \frac{\text{Volts}}{\sqrt{\mu m}} \\). With \\( a \\) denoting the total width of the array (the variable `physicalSize` above), the total power input into the MLA's is$$ P = \frac{\int_{-a/2}^{a/2} \left| u \left( x \right) \right|^2 dx}{Z_{0}} $$which has units of \\( \frac{\text{Volts}^2}{\Omega} \\). \\( Z_0 \\) is the [impedance of free space](https://en.wikipedia.org/wiki/Impedance_of_free_space), \\( Z_0 \approx 376.73 \, \Omega \\).<jupyter_code>powerIn = 100 # mW
Z0      = 376.73 # Impedance of free space, Ohms

# Determine the field amplitude
fieldAmp = np.sqrt((powerIn / 1000) * Z0 / physicalSize) # Volts / sqrt(microns)

# Plane wave
uIn = lambda x: fieldAmp
uIn = np.vectorize(uIn)

# Plot the input field over the array
plt.plot(grid.px, uIn(grid.px), linewidth = 2)
plt.xlabel(r'x-position, $\mu m$')
plt.ylabel(r'Field amplitude, $V / \sqrt{\mu m}$')
plt.grid(True)
plt.show()<jupyter_output><empty_output><jupyter_text>Find the field immediately after the second MLAThe field immediately *after* the second MLA is a summation of the Fourier transforms of all the fields sampled by each lenslet in the first array. Each transform is centered around the axis of its corresponding lenslet. Therefore, we have to compute a Fourier transform for each lenslet separately and then shift the origin of its coordinate axes onto the axis for the corresponding lenslet. Mathematically, the field immediately after the second MLA is$$ u_{1}' \left( x' \right) = \sum_{m = -\frac{N - 1}{2}}^{\frac{N-1}{2}} \frac{1}{\sqrt{j \lambda f_{LA}}} \int_{-\infty}^{\infty} u \left( x \right) \text{rect} \left( \frac{x - mp}{p}\right) \exp \left[ - j \frac{2 \pi}{\lambda f_{LA}} \left( x - mp \right) \left( x' + mp \right) \right] dx$$with \\(N \\) odd, \\( f_{LA} \\) representing the lenslet focal length, and \\( \text{rect} \left( x \right) \\) equal to one between \\(-1/2 \\) and \\(+ 1/2 \\) and zero elsewhere. The multiplication of the input field by the \\( \text{rect} \\) function models a single lenslet sampling a region of the field. The sampled field is then Fourier transformed by the integral with the Fourier transform centered on the center coordinates for the respective lenslet. The transformed fields are then summed over all \\( N \\) lenslets.The reason that there is no quadratic phase curvature preceding the transforms is because the second MLA acts as a field lens, effectively canceling the quadratic phase terms. In doing so, we ignore regions of the diffracted field that leak onto neighboring lenslets because1. they carry little energy relative to the center of the diffracted field, and2. the tails of the diffracted fields contain aliasing artifacts anyway.To compute this spatially-parallelized Fourier transform, we use SimMLA's fftSubgrid routine. It will return two lists of interpolated fields, one for the magnitude and one for the phase. We will then resample this field onto a new grid representing the coordinate axes immediately after the second MLA.<jupyter_code># Compute the interpolated fields
# Linear interpolation is used for speed
%time interpMag, interpPhase = simfft.fftSubgrid(uIn, grid)<jupyter_output>CPU times: user 5min 4s, sys: 6.45 s, total: 5min 11s
Wall time: 5min 10s<jupyter_text>Now that we have the interpolated fields, we'll define a new grid for sampling them. We'll give the grid the same physical extent as before, but we may choose to increase the spatial sampling rate slightly. We will also use the focal length of the objective because it will serve as the final Fourier transforming lens in the system.*If we increase the sampling rate for the new grid, it won't necessarily fix any subsampling that occurred in the previous step. It will only provide better resolution at the possibly aliased results.*<jupyter_code># Check that we've successfully computed the Fourier transform of a lenslet.
m = 6 # lenslet index, integers between [0, numLenslets - 1]
plt.plot(grid.px, np.abs(interpMag[m](grid.px) * np.exp(1j * interpPhase[m](grid.px))))
plt.xlabel(r'x-position, $\mu m$')
plt.ylabel(r'Field amplitude, $V / \sqrt{\mu m}$')
plt.xlim((0,400))
plt.grid(True)
plt.show()
# Check that we've successfully computed the Fourier transform of a lenslet.
m = 6 # lenslet index, integers between [0, numLenslets - 1]
plt.plot(grid.px, np.angle(interpMag[m](grid.px) * np.exp(1j * interpPhase[m](grid.px))),'.')
plt.xlabel(r'x-position, $\mu m$')
plt.ylabel(r'Phase, radians')
plt.xlim((0,400))
plt.grid(True)
plt.show()<jupyter_output><empty_output><jupyter_text>Most values for the phase should be 0 or \\( \pi \\). In this example, any phase values not equal to either of these are probably a result of aliasing, which simply cannot be avoided completely using the FFT.<jupyter_code>%%time
fObj        = 100000 # microns
newGridSize = subgridSize * numLenslets # microns

# Upsample grid by an odd factor to reduce rippling
newGrid = grids.Grid(5*newGridSize, 5*physicalSize, wavelength, fObj, dim = 1)
field   = np.zeros(newGrid.gridSize)

# For each interpolated magnitude and phase corresponding to a lenslet
# 1) Compute the full complex field
# 2) Sum it with the other complex fields
for currMag, currPhase in zip(interpMag, interpPhase):
    fieldMag   = currMag(newGrid.px)
    fieldPhase = currPhase(newGrid.px)
    
    currField = fieldMag * np.exp(1j * fieldPhase)
    field     = field + currField<jupyter_output>CPU times: user 12.7 s, sys: 2.44 s, total: 15.2 s
Wall time: 15.2 s<jupyter_text>Let's plot the field amplitude and phase just immediately after the second MLA. Note that there are 11 focused spots corresponding the focus of each microlens array.<jupyter_code>fig, (ax0, ax1) = plt.subplots(nrows = 1,
                               ncols = 2,
                               sharex = True,
                               sharey = False,
                               figsize = (10,6))

ax0.plot(newGrid.px, np.abs(field))
ax0.set_xlim((-7500, 7500))
ax0.set_xlabel(r'x-position, $\mu m$')
ax0.set_ylabel(r'Amplitude, $V / \sqrt{\mu m}$')
ax0.grid()

ax1.plot(newGrid.px, np.angle(field), '.')
ax1.set_xlabel(r'x-position, $\mu m$')
ax1.set_ylabel('Phase, radians')
ax1.grid()

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>We can also check that energy has been conserved, at least to within numerical round off errors. The power delivered by the diffracted field \\( u_{1}' \left( x' \right) \\) just past the second MLA is$$ P_{1}' = \frac{\int_{-\infty}^{\infty} \left| u \left( x' \right) \right|^2 dx'}{Z_{0}} $$Note that the diffracted field carries a prefactor of \\( \frac{1}{\sqrt{\lambda f_{LA}}} \\) which is taken care of by SimMLA.<jupyter_code>from scipy.integrate import simps
# Factor of 1000 converts from Watts to milliWatts
power1Prime = simps(np.abs(field)**2, newGrid.px) / Z0 * 1000 

print('Input power is: {0:.4f} mW'.format(powerIn))
print('Diffracted power is: {0:.4f} mW'.format(power1Prime))<jupyter_output>Input power is: 100.0000 mW
Diffracted power is: 100.4267 mW<jupyter_text>Find the field on the sample Compute the field in the objective's back focal planeThe final step is to compute the field formed in the sample plane by Fourier transforming the field just behind the second MLA.  Compute the field on the sample<jupyter_code># Sample plane irradiance including limiting aperture in the BFP
scalingFactor = newGrid.physicalSize / (newGrid.gridSize - 1) / np.sqrt(newGrid.wavelength * newGrid.focalLength)
F             = scalingFactor * np.fft.fftshift(np.fft.fft(np.fft.ifftshift(field)))
plt.plot(newGrid.pX, np.abs(F)**2 / Z0 * 1000)
plt.xlim((-15000, 15000))
plt.xlabel(r'x-position, $\mu m$')
plt.ylabel(r'Irradiance, $mW / \mu m$')
plt.grid(True)
plt.show()

powerSample = simps(np.abs(F)**2, newGrid.pX) / Z0 * 1000
print('Power delivered to sample: {0:.4f} mW'.format(powerSample))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10200953678474114, 'fraction_numerical': 0.02111716621253406, 'mean_word_length': 4.288158487167943, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 38, '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': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '360505', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3794, 'n_tokens_neox': 3480, 'n_words': 1480, 'path': 'tests/Plane Wave Incidence Verification.ipynb', 'repo_name': 'nolsman/simmla', 'size': '11757'}