Python version of Joseph Goodman's Mathematica code for simulating laser speckle. Both free space and in an optical system are combined into the one file.

Goodman_speckle_test.py

# my conversion of Goodman's code
import numpy as np
import matplotlib.pyplot as plt

# Program for simulating Speckle Formation by Free Space Propagation
# Translate from Mathematica
# Ref: Goodman, Joseph W. 2020. Speckle Phenomena in Optics: Theory and Applications, Second Edition. SPIE.


n = 1024   # nxn array size
k = 16  # samples per speckle

r_0 = n // k
# generate a n/k x n/k array of random phasors
start = np.exp(np.random.rand(n//k, n//k) * 2*np.pi*1j)
# pad with zeros
scatter_array = np.pad(start, ((0, n - r_0), (0, n - r_0)), 'constant')
# Calculate the speckle field in the observation plane
speckle_field = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(scatter_array)))
# Calculate the intensity
speckle_int = abs(speckle_field)**2
plt.figure(1)
plt.imshow(speckle_int)
plt.colorbar
plt.title('Speckle Intensity')

## Now simulate the Speckle formation in imaging

n = 1024   # nxn array size
k = 8  # samples per speckle
r_0 = n // k  # radius of lens pupil function in pixels


# create a simple image of rectangles
object_intensity = np.ones((n, n))
object_intensity[0:341, :] = 0.1
object_intensity[342:684, :] = 0.5

# generate the random phasors
random_field = np.exp(np.random.rand(n, n) * 2*np.pi*1j)
# calculate the object's amplitude
scatter_field = np.sqrt(object_intensity)*random_field

plt.figure(2)
plt.imshow(abs(scatter_field)**2, cmap='gray')
plt.colorbar()
plt.title('Intensity of the scattering surface')


# generate the circular pupil function of the lens
bandpass = np.zeros((n, n))
xv = np.arange(n)
xx, yy = np.meshgrid(xv, xv)
R = np.sqrt((xx-n/2) ** 2 + (yy-n/2) ** 2)/r_0
bandpass[R <= 1] = 1

plt.figure(3)
plt.imshow(bandpass, cmap='gray')
plt.colorbar()
plt.title('Pupil')

# Calculate the field transmitted by the lens pupil
pupil_field = bandpass*np.fft.fft2(scatter_field)
# Calculate the image field
image_field = np.fft.ifft2(pupil_field)
# Calculate the image intensity
image_intensity = np.abs(image_field)**2


plt.figure(4)
plt.imshow(image_intensity, cmap='gray')
plt.colorbar()
plt.title('Image Intensity')