diff --git a/meanas/fdfd/farfield.py b/meanas/fdfd/farfield.py
index 86ec0d7..0051cd0 100644
--- a/meanas/fdfd/farfield.py
+++ b/meanas/fdfd/farfield.py
@@ -78,15 +78,12 @@ def near_to_farfield(
     kx, ky = numpy.meshgrid(kxx, kyy, indexing='ij')
     kxy2 = kx * kx + ky * ky
     kxy = numpy.sqrt(kxy2)
-    kz = numpy.sqrt(k * k - kxy2)
+    kz = numpy.sqrt(numpy.maximum(0, k * k - kxy2))
 
-    sin_th = ky / kxy
-    cos_th = kx / kxy
+    sin_th = numpy.divide(ky, kxy, out=numpy.zeros_like(ky), where=kxy != 0)
+    cos_th = numpy.divide(kx, kxy, out=numpy.ones_like(kx), where=kxy != 0)
     cos_phi = kz / k
 
-    sin_th[numpy.logical_and(kx == 0, ky == 0)] = 0
-    cos_th[numpy.logical_and(kx == 0, ky == 0)] = 1
-
     # Normalized vector potentials N, L
     N = [-Hn_fft[1] * cos_phi * cos_th + Hn_fft[0] * cos_phi * sin_th,
           Hn_fft[1] * sin_th + Hn_fft[0] * cos_th]                      # noqa: E127
@@ -114,8 +111,8 @@ def near_to_farfield(
     outputs = {
         'E': E_far,
         'H': H_far,
-        'dkx': kx[1] - kx[0],
-        'dky': ky[1] - ky[0],
+        'dkx': float(kxx[1] - kxx[0]),
+        'dky': float(kyy[1] - kyy[0]),
         'kx': kx,
         'ky': ky,
         'theta': theta,
@@ -177,22 +174,19 @@ def far_to_nearfield(
     padded_shape = cast('Sequence[int]', padded_size)
 
     k = 2 * pi
-    kxs = fftshift(fftfreq(s[0], 1 / (s[0] * dkx)))
-    kys = fftshift(fftfreq(s[0], 1 / (s[1] * dky)))
+    kxs = dkx * fftshift(fftfreq(s[0], d=1 / s[0]))
+    kys = dky * fftshift(fftfreq(s[1], d=1 / s[1]))
 
     kx, ky = numpy.meshgrid(kxs, kys, indexing='ij')
     kxy2 = kx * kx + ky * ky
     kxy = numpy.sqrt(kxy2)
 
-    kz = numpy.sqrt(k * k - kxy2)
+    kz = numpy.sqrt(numpy.maximum(0, k * k - kxy2))
 
-    sin_th = ky / kxy
-    cos_th = kx / kxy
+    sin_th = numpy.divide(ky, kxy, out=numpy.zeros_like(ky), where=kxy != 0)
+    cos_th = numpy.divide(kx, kxy, out=numpy.ones_like(kx), where=kxy != 0)
     cos_phi = kz / k
 
-    sin_th[numpy.logical_and(kx == 0, ky == 0)] = 0
-    cos_th[numpy.logical_and(kx == 0, ky == 0)] = 1
-
     theta = numpy.arctan2(ky, kx)
     phi = numpy.arccos(cos_phi)
     theta[numpy.logical_and(kx == 0, ky == 0)] = 0
@@ -212,21 +206,41 @@ def far_to_nearfield(
     N = [L[1],
         -L[0]]  # noqa: E128
 
-    En_fft = [-( L[0] * sin_th + L[1] * cos_phi * cos_th) / cos_phi,
-              -(-L[0] * cos_th + L[1] * cos_phi * sin_th) / cos_phi]
+    En_fft = [
+        numpy.divide(
+            -(L[0] * sin_th + L[1] * cos_phi * cos_th),
+            cos_phi,
+            out=numpy.zeros_like(L[0]),
+            where=cos_phi != 0,
+            ),
+        numpy.divide(
+            -(-L[0] * cos_th + L[1] * cos_phi * sin_th),
+            cos_phi,
+            out=numpy.zeros_like(L[0]),
+            where=cos_phi != 0,
+            ),
+        ]
 
-    Hn_fft = [( N[0] * sin_th + N[1] * cos_phi * cos_th) / cos_phi,
-              (-N[0] * cos_th + N[1] * cos_phi * sin_th) / cos_phi]
-
-    for i in range(2):
-        En_fft[i][cos_phi == 0] = 0
-        Hn_fft[i][cos_phi == 0] = 0
+    Hn_fft = [
+        numpy.divide(
+            N[0] * sin_th + N[1] * cos_phi * cos_th,
+            cos_phi,
+            out=numpy.zeros_like(N[0]),
+            where=cos_phi != 0,
+            ),
+        numpy.divide(
+            -N[0] * cos_th + N[1] * cos_phi * sin_th,
+            cos_phi,
+            out=numpy.zeros_like(N[0]),
+            where=cos_phi != 0,
+            ),
+        ]
 
     E_near = [ifftshift(ifft2(ifftshift(Ei), s=padded_shape)) for Ei in En_fft]
     H_near = [ifftshift(ifft2(ifftshift(Hi), s=padded_shape)) for Hi in Hn_fft]
 
     dx = 2 * pi / (s[0] * dkx)
-    dy = 2 * pi / (s[0] * dky)
+    dy = 2 * pi / (s[1] * dky)
 
     outputs = {
         'E': E_near,
@@ -236,4 +250,3 @@ def far_to_nearfield(
     }
 
     return outputs
-