<jupyter_start><jupyter_text>Verification of the partial coherence phase screens<jupyter_code>%pylab
%matplotlib inline
import SimMLA.fftpack as simfft
import SimMLA.grids   as grids
import SimMLA.fields  as fields
from numpy.fft import fft, ifft, fftshift, ifftshift
from scipy.integrate import simps<jupyter_output>Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib<jupyter_text>I am simulating a Gaussian Schell model beam using a technique described in [Xiao and Voelz, "Wave optics simulation approach for partial spatially coherent beams," Opt. Express 14, 6986-6992 (2006)](https://www.osapublishing.org/oe/abstract.cfm?uri=oe-14-16-6986) and Chapter 9 of [Computational Fourier Optics: A MATLAB Tutorial](http://spie.org/Publications/Book/858456), which is also by Voelz. This workbook verifies that we can produce the same results for the partial spatial coherence simulation as presented on page 186 of Voelz's book.**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.**  Build the coordinate system<jupyter_code>focalLength = 250e-3

wavelength  = 650e-9

gridSize     = 251   # Number of grid (or lattice) sites
physicalSize = 50e-3 # The full extent of the grid

# dim = 1 makes the grid 1D
grid = grids.Grid(gridSize, physicalSize, wavelength, focalLength, dim = 1)<jupyter_output><empty_output><jupyter_text>The desired coherence length is set in `coherenceLength`. The rest of the beam is defined as described in another notebook.<jupyter_code>Z0              = 376.73 # Impedance of free space, Ohms
power           = 100  # mW
beamStd         = 1e-3
coherenceLength = 8e-3
fieldAmp = np.sqrt(power / 1000 * Z0 / beamStd / np.sqrt(np.pi)) # Factor of 1000 converts from mW to W

beam     = fields.GSMBeamRealization(fieldAmp, beamStd, coherenceLength, grid)<jupyter_output><empty_output><jupyter_text>Let's go ahead and plot one realization of the phase angle generated from this beam. We won't use the Gaussian envelope; only the phase.<jupyter_code>beamSample = beam(grid.px)

plt.plot(grid.px, np.angle(beamSample), linewidth = 2)
plt.xlabel(r'x-position')
plt.ylabel(r'Field phase, rad')
plt.grid(True)
plt.show()<jupyter_output><empty_output><jupyter_text>Computing the two-slit diffraction patternFirst, let's define the field. We have a plane wave incident on two slits symmetrically placed about the origin.<jupyter_code>slitWidth  = 2e-3 
slitOffset = 2.5e-3
beamPower  = 100  # mW

Z0    = 376.73 # Impedance of free space, Ohms
amp   = np.sqrt((beamPower / 1000) * Z0 / slitWidth / 2)
field = np.zeros(grid.px.size)
field[np.logical_and(grid.px >  (-slitOffset - slitWidth / 2), grid.px <= (-slitOffset + slitWidth / 2))] = amp
field[np.logical_and(grid.px >= ( slitOffset - slitWidth / 2), grid.px <  ( slitOffset + slitWidth / 2))] = amp
plt.plot(grid.px, field, linewidth = 2)
plt.xlim((-5e-3, 5e-3))
plt.xlabel(r'x-coordinate, $\mu m$')
plt.ylabel(r'Field amplitude, $V / \sqrt{\mu m}$')
plt.grid(True)
plt.show()<jupyter_output><empty_output><jupyter_text>Next, let's compute the diffraction pattern from this field distribution.<jupyter_code>dx = grid.px[1] - grid.px[0]
F  = fftshift(fft(ifftshift(field))) * dx
Irrad = np.abs(F)**2 / wavelength / focalLength
plt.plot(grid.pX, Irrad)
plt.xlim((-2e-4,2e-4))
plt.xlabel('x-position')
plt.ylabel('Irradiance, $W / m$')
plt.grid(True)
plt.show()
# Check power conservation
intPower = simps(Irrad, grid.pX)
print(intPower / Z0 * 1000)<jupyter_output>100.399980538<jupyter_text>Add partial spatial coherence to the incident beam.Now that we have a two slit interferometer to experiment with, we can begin adjust the spatial coherence of the beam.Let's first start by computing the random realization of the phase screen and apply it to the incident field.<jupyter_code># Generate a new phase screen
beamSample = beam(grid.px)

t = np.exp(1j * np.angle(beamSample))
newField = t * field
# Plot the real part of the field
plt.plot(grid.px, np.real(newField))
plt.show()
# Propagate and plot the new field
dx = grid.px[1] - grid.px[0]
F  = fftshift(fft(ifftshift(newField))) * dx
newIrrad = np.abs(F)**2 / wavelength / focalLength

plt.plot(grid.pX, newIrrad)
plt.plot(grid.pX, Irrad)
#plt.xlim((-2e-4,2e-4))
plt.xlabel('x-position')
plt.ylabel('Irradiance, $W / m$')
plt.grid(True)
plt.show()<jupyter_output><empty_output><jupyter_text>Compute the averaged irradiance for multiple realizations<jupyter_code>nIter = 1000
dx = grid.px[1] - grid.px[0]

finalIrrad = np.zeros(grid.px.size)
for ctr in range(nIter):
    # Create a new realization of the field
    t = np.exp(1j * np.angle(beam(grid.px)))
    newField = t * field
    
    F  = fft(ifftshift(newField)) * dx
    Irrad = np.abs(F)**2 / wavelength / focalLength
    
    finalIrrad = finalIrrad + Irrad
    
# Find the averaged irradiance pattern    
finalIrrad = fftshift(finalIrrad) / nIter
plt.plot(grid.pX, finalIrrad)
#plt.xlim((-2e-4,2e-4))
plt.xlabel('x-position')
plt.ylabel('Irradiance, $W / m$')
plt.grid(True)
plt.show()
# Check the integrated power
intPower = simps(newIrrad, grid.pX)
print(intPower / Z0 * 1000)<jupyter_output>100.3996561<jupyter_text>Compute the degree of coherence from this profileThe degree of coherence may be obtained by fitting the analytical results to this irradiance profile for \\( \left| \mu_{12} \right| \\):$$ I \left( x' \right) = \frac{2 a^2 E_0^2}{\lambda f} \text{sinc}^2 \left( \frac{a x'}{\lambda z} \right) \times \left[ 1 + \left| \mu_{12} \right| \cos \left( 2 \pi \frac{\Delta s}{\lambda f} x' \right) \right] $$Where \\( a \\) is the width of a single slit, \\( E_0 \\) is the field amplitude, \\( \Delta s \\) is the separation between the two slits, and \\( \text{sinc} \left( x \right) = \frac{\sin \pi x}{\pi x} \\). The analytical result only holds when the coherence length is on the order of or larger than the beam diameters.<jupyter_code>from scipy.optimize import curve_fit as cf

def sinc(x):
    if (x != 0):
        # Prevent divide-by-zero
        return np.sin(np.pi * x) / (np. pi * x)
    else:
        return 1
sinc = np.vectorize(sinc)

def theory(x, *p):
    mu = p
    
    return 2 * slitWidth**2 * amp**2 / wavelength / focalLength * sinc(x * slitWidth / wavelength / focalLength)**2 \
        * (1 + mu * np.cos(2 * np.pi * (slitOffset * 2) * x / wavelength / focalLength))
initGuess = 0.6
popt, _ = cf(theory, grid.pX, finalIrrad, p0 = initGuess)
mu = popt[0]

plt.plot(grid.pX, finalIrrad, linewidth = 2, label = 'FFT')
plt.plot(grid.pX, theory(grid.pX, mu), '--', linewidth = 3, label = 'Fit, d.o.c. = {0:.2f}'.format(mu))
#plt.xlim((-2e-4,2e-4))
plt.xlabel('x-position')
plt.ylabel('Irradiance, $W / m$')
plt.grid(True)
plt.legend()
plt.show()


numericCohLength = np.sqrt(- (2 * slitOffset)**2 / np.log(mu))
print('The fitted coherence length is: {0:.4f}'.format(numericCohLength))
print('The input coherence length is: {0:.4f}'.format(coherenceLength))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13139617292700212, 'fraction_numerical': 0.028348688873139617, 'mean_word_length': 4.2383073496659245, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 36, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, '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': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '608673', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2756, 'n_tokens_neox': 2501, 'n_words': 735, 'path': 'tests/GSM Verification.ipynb', 'repo_name': 'nolsman/simmla', 'size': '7055'}