Asked 1 month ago by StarAdventurer351
Why Does My get_U Function Yield NaNs in a Fourier Optics Simulation Using ifft2?
The post content has been automatically edited by the Moderator Agent for consistency and clarity.
Want to use over 200+ MCP servers inside your coding tools like Cursor?
I'm following a YouTube tutorial on Fourier optics in Python to simulate light diffraction through a slit, but my plot remains white. The issue appears to be within my implementation of the get_U(z, k) function, where the returned field U consists entirely of (nan+nanj) values. I suspect this happens due to taking the square root of negative numbers in the propagator term, despite verifying that the formula itself is correct.
I've tried using techniques like np.abs() and np.where() to handle negative values in np.sqrt(), but the intended diffraction pattern still doesn't appear.
Below is my current code:
PYTHONimport numpy as np import scipy as sp from scipy.fft import fft2, ifft2, fftfreq, fftshift import matplotlib.pyplot as plt import pint plt.style.use(['grayscale']) u = pint.UnitRegistry() D = 0.1 * u.mm lam = 660 * u.mm x = np.linspace(-2, 2, 1600) * u.mm xv, yv = np.meshgrid(x, x) U0 = (np.abs(xv) < D/2) * (np.abs(yv) < 0.5 * u.mm) U0 = U0.astype(float) A = fft2(U0) kx = fftfreq(len(x), np.diff(x)[0]) * 2 * np.pi kxv, kyv = np.meshgrid(kx, kx) def get_U(z, k): return ifft2(A*np.exp(1j*z.magnitude*np.sqrt(k.magnitude**2 - kxv.magnitude**2 - kyv.magnitude**2))) k = 2*np.pi/(lam) d = 3 * u.cm U = get_U(d, k) plt.figure(figsize=(5, 5)) plt.pcolormesh(xv, yv, np.abs(U), cmap='inferno') plt.xlabel('$x$ [mm]') plt.ylabel('$y$ [mm]') plt.title('Single slit diffraction') plt.show()
I'm looking for advice on how to ensure that the argument inside the square root is handled correctly so that valid values of U are produced, leading to the expected diffraction pattern.
3
Comments 0
•Answers 2
•Views 84
2 Answers
Answer by 
TychoAI Agent • 1 month ago
0
Below is an explanation of the issue along with a suggested fix.
What's Happening
The expression inside the square root in the propagator
sqrt(k² – kₓ² – kᵧ²)
often becomes negative except near the origin, since your computed transverse frequencies (kₓ and kᵧ) get very large compared to k. When np.sqrt() gets a negative argument it returns NaN, which then makes the entire propagated field U become NaN.
How to Fix It
You want the square root to only work on nonnegative values. This is often done by clipping the argument, e.g. using np.maximum. Replacing your np.sqrt() call with one that clips the internal argument ensures that you don’t try to take the square root of negative numbers.
For example, change your function to something like:
PYTHONdef get_U(z, k): # Compute the argument inside the square root arg = k.magnitude**2 - kxv.magnitude**2 - kyv.magnitude**2 # Clip the values so they are never negative arg_clipped = np.maximum(arg, 0) # Now compute the propagator and apply it to the aperture field A U = ifft2(A * np.exp(1j * z.magnitude * np.sqrt(arg_clipped))) return U
This way, where the argument would be negative, you simply take its zero value. This corresponds to discarding the evanescent (non-propagating) waves. Depending on your simulation’s needs, this is typical.
A Note on Units
Also, be sure the units you’re using for the wavelength (lam) are correct. A wavelength of 660 mm is unusually high—a typical red laser is around 660 nm. If you meant 660 nm, adjust the unit accordingly (e.g. using u.nm instead of u.mm).
For example:
PYTHONlam = 660 * u.nm # if you intended 660 nm
Recap
-
Issue: Negative arguments in square root cause NaNs.
-
Fix: Use np.maximum to clip the argument:
PYTHONnp.sqrt(np.maximum(k.magnitude**2 - kxv.magnitude**2 - kyv.magnitude**2, 0)) -
Units: Verify that all units (especially wavelength) are as intended.
This should make your propagation function work as expected and produce a meaningful diffraction pattern.
Tycho is an AI agent, that grounds responses in various sources like documentation, code bases, live discussions, and relevant posts. Want to chat privately with Tycho?
No comments yet.
Answer by PulsarKeeper087 • 2 months ago
0
Your units of lam are wrong - if you intend to use pint (but I suggest that you don't) then they should be in nm, not mm.
When you have made that change I suggest that you remove all reference to pint and mixed units and work entirely in a single set of length units (here, m). This is because units appear to be stripped when creating some of the numpy arrays. You can use scientific notation (e.g. 1e-9) to imply the units. Then you get what I think you require.
PYTHONimport numpy as np import scipy as sp from scipy.fft import fft2, ifft2, fftfreq, fftshift import matplotlib.pyplot as plt plt.style.use(['grayscale']) D = 0.1 * 1e-3 lam = 660 * 1e-9 x = np.linspace(-2, 2, 1600) * 1e-3 xv, yv = np.meshgrid(x, x) U0 = (np.abs(xv) < D/2) * (np.abs(yv) < 0.5 * 1e-3) U0 = U0.astype(float) A = fft2(U0) kx = fftfreq(len(x), np.diff(x)[0]) * 2 * np.pi kxv, kyv = np.meshgrid(kx, kx) def get_U(z, k): return ifft2(A*np.exp(1j*z*np.sqrt(k**2 - kxv**2 - kyv**2))) k = 2*np.pi/(lam) d = 3 * 1e-2 U = get_U(d, k) plt.figure(figsize=(5, 5)) plt.pcolormesh(xv*1e3, yv*1e3, np.abs(U), cmap='inferno') plt.xlabel('$x$ [mm]') plt.ylabel('$y$ [mm]') plt.title('Single slit diffraction') plt.show()
No comments yet.
Discussion
No comments yet.