BPM and WPM fixes for odd grids

diffractio/scalar_fields_X.py

@@ -1802,7 +1802,7 @@ def PWD_kernel(u: NDArrayComplex, n: NDArrayComplex, k0: float, k_perp2: NDArray
     Ek = fftshift(fft(u))
     H = np.exp(1j * dz * csqrt(n**2 * k0**2 - k_perp2))
 
-    return ifft(fftshift(H * Ek))
+    return ifft(np.fft.ifftshift(H * Ek))
 
 
 def WPM_schmidt_kernel(u, n: NDArrayComplex, k0: float, k_perp2: NDArrayComplex, dz: float):

diffractio/scalar_fields_XY.py

@@ -3234,7 +3234,7 @@ def PWD_kernel(u: NDArrayComplex, n: float, k0: float, k_perp2: NDArrayComplex,
     Ek = fftshift(fft2(u))
     H = np.exp(1j * dz * csqrt(n**2 * k0**2 - k_perp2) - absorption)
 
-    result = (ifft2(fftshift(H * Ek)))
+    result = ifft2(np.fft.ifftshift(H * Ek))
     return result
 
 

diffractio/scalar_fields_XYZ.py

@@ -889,21 +889,19 @@ def BPM(self, has_edges: bool = True, pow_edge: int = 80, verbose: bool = False)
         numy = len(self.y)
         numz = len(self.z)
 
-        deltax = self.x[-1] - self.x[0]
-        deltay = self.y[-1] - self.y[0]
+        deltax = self.x[1] - self.x[0]
+        deltay = self.y[1] - self.y[0]
         deltaz = self.z[1] - self.z[0]
         # pixelx = np.linspace(-int(numx/2), int(numx/2), numx)
         # pixely = np.linspace(-numy/2, numy/2, numy)
 
         modo = self.u0.u
 
-        kx1 = np.linspace(0, int(numx/2) + 1, int(numx/2))
-        kx2 = np.linspace(-int(numx/2), -1, int(numx/2))
-        kx = (2 * np.pi / deltax) * np.concatenate((kx1, kx2))
+        # Use fftfreq to guarantee kx has exactly numx samples.
+        kx = 2 * np.pi * np.fft.fftfreq(numx, d=deltax)
 
-        ky1 = np.linspace(0, numy/2 + 1, int(numy/2))
-        ky2 = np.linspace(-numy/2, -1, int(numy/2))
-        ky = (2 * np.pi / deltay) * np.concatenate((ky1, ky2))
+        # Use fftfreq to guarantee ky has exactly numy samples.
+        ky = 2 * np.pi * np.fft.fftfreq(numy, d=deltay)
 
         KX, KY = np.meshgrid(kx, ky)
 

diffractio/scalar_fields_XZ.py

@@ -779,17 +779,15 @@ def __BPM__(self,
         numz = len(self.z)  # distance en z
         numx = len(self.x)  # distance en x
         deltaz = self.z[1] - self.z[0]  # Tamaño del sampling
-        rangox = self.x[-1] - self.x[0]
+        deltax = self.x[1] - self.x[0]
 
         pixelx = np.linspace(-int(numx/2), int(numx/2), numx)
         # initial field
         field_z = self.u0.u
 
         # Calculo de la phase 1 normalizada -------------------
-        kx1 = np.linspace(0, int(numx/2) + 1, int(numx/2))
-        kx2 = np.linspace(-int(numx/2), -1, int(numx/2))
-        # Número de ondas del material en una dimensión
-        kx = (2 * np.pi / rangox) * np.concatenate((kx1, kx2))
+        # Use fftfreq to guarantee kx has exactly numx samples (odd/even-safe).
+        kx = 2 * np.pi * np.fft.fftfreq(numx, d=deltax)
         # Función de transferencia para la propagación que es identica
         # a la respuesta de frecuencia espacial en óptica de Fourier
         # incorporando el termino (-j k0 z).