[farfield] preserve centered geometry and physical Fourier normalization

meanas/fdfd/farfield.py

@@ -1,8 +1,9 @@
 """
 Functions for performing near-to-farfield transformation (and the reverse).
 """
-from typing import Any, cast
+from typing import Any
 from collections.abc import Sequence
+import operator
 import numpy
 from numpy.fft import fft2, fftshift, fftfreq, ifft2, ifftshift
 from numpy import pi
@@ -13,6 +14,34 @@ type farfield_slice = NDArray[complexfloating]
 type transverse_slice_pair = Sequence[farfield_slice]
 
 
+def _padded_shape(shape: tuple[int, ...], padded_size: Sequence[int] | int | None) -> tuple[int, int]:
+    if len(shape) != 2 or any(n == 0 for n in shape):
+        raise ValueError('fields must be nonempty two-dimensional arrays')
+    if padded_size is None:
+        sizes = tuple(1 << (n - 1).bit_length() for n in shape)
+    else:
+        raw_sizes = (padded_size, padded_size) if numpy.isscalar(padded_size) else padded_size
+        sizes = tuple(operator.index(n) for n in raw_sizes)
+    if len(sizes) != 2 or any(n <= 0 for n in sizes):
+        raise ValueError('padded_size must contain two positive integer sizes')
+    return sizes[0], sizes[1]
+
+
+def _resize_centered(field: farfield_slice, shape: tuple[int, int]) -> farfield_slice:
+    """Pad or crop while preserving the sample at index size // 2 as the origin."""
+    if field.shape == shape:
+        return field
+    result = numpy.zeros(shape, dtype=field.dtype)
+    source, target = [], []
+    for old, new in zip(field.shape, shape, strict=True):
+        start_old, start_new = max(0, old // 2 - new // 2), max(0, new // 2 - old // 2)
+        count = min(old - start_old, new - start_new)
+        source.append(slice(start_old, start_old + count))
+        target.append(slice(start_new, start_new + count))
+    result[tuple(target)] = field[tuple(source)]
+    return result
+
+
 def near_to_farfield(
         E_near: transverse_slice_pair,
         H_near: transverse_slice_pair,
@@ -24,7 +53,9 @@ def near_to_farfield(
     Compute the farfield, i.e. the distribution of the fields after propagation
       through several wavelengths of uniform medium.
 
-    The input fields should be complex phasors.
+    The input fields should be complex phasors. The spatial origin is the
+    sample at ``(nx // 2, ny // 2)``. Padding/cropping preserves that origin.
+    Fourier integrals include the cell area ``dx * dy``.
 
     Args:
         E_near: List of 2 ndarrays containing the 2D phasor field slices for the transverse
@@ -40,9 +71,9 @@ def near_to_farfield(
     Returns:
         Dict with keys
 
-        -   `E_far`: Normalized E-field farfield; multiply by
+        -   `E`: Normalized E-field farfield; multiply by
                 (i k exp(-i k r) / (4 pi r)) to get the actual field value.
-        -   `H_far`: Normalized H-field farfield; multiply by
+        -   `H`: Normalized H-field farfield; multiply by
                 (i k exp(-i k r) / (4 pi r)) to get the actual field value.
         -   `kx`, `ky`: Wavevector values corresponding to the x- and y- axes in E_far and H_far,
                 normalized to wavelength (dimensionless).
@@ -62,14 +93,10 @@ def near_to_farfield(
     if not all(s == f.shape for f in [*E_near, *H_near]):
         raise Exception('All fields must be the same shape!')
 
-    if padded_size is None:
+    padded_shape = _padded_shape(s, padded_size)
-        padded_size = (2**numpy.ceil(numpy.log2(s))).astype(int)
-    if not hasattr(padded_size, '__len__'):
-        padded_size = (padded_size, padded_size)            # type: ignore  # checked if sequence
-    padded_shape = cast('Sequence[int]', padded_size)
 
-    En_fft = [fftshift(fft2(fftshift(Eni), s=padded_shape)) for Eni in E_near]
+    En_fft = [fftshift(fft2(ifftshift(_resize_centered(Eni, padded_shape)))) * dx * dy for Eni in E_near]
-    Hn_fft = [fftshift(fft2(fftshift(Hni), s=padded_shape)) for Hni in H_near]
+    Hn_fft = [fftshift(fft2(ifftshift(_resize_centered(Hni, padded_shape)))) * dx * dy for Hni in H_near]
 
     # Propagation vectors kx, ky
     k  = 2 * pi
@@ -112,8 +139,8 @@ def near_to_farfield(
     outputs = {
         'E': E_far,
         'H': H_far,
-        'dkx': float(kxx[1] - kxx[0]),
+        'dkx': float(2 * pi / (padded_shape[0] * dx)),
-        'dky': float(kyy[1] - kyy[0]),
+        'dky': float(2 * pi / (padded_shape[1] * dy)),
         'kx': kx,
         'ky': ky,
         'theta': theta,
@@ -131,10 +158,11 @@ def far_to_nearfield(
         padded_size: Sequence[int] | int | None = None
         ) -> dict[str, Any]:
     """
-    Compute the farfield, i.e. the distribution of the fields after propagation
+    Reconstruct the propagating near-field spectrum from a normalized farfield.
-      through several wavelengths of uniform medium.
 
-    The input fields should be complex phasors.
+    Inputs are not modified. Spectral padding/cropping preserves zero frequency
+    at index ``size // 2``. The output spacings follow its resized shape, and
+    inverse Fourier quadrature restores the cell-area normalization.
 
     Args:
         E_far: List of 2 ndarrays containing the 2D phasor field slices for the transverse
@@ -168,11 +196,7 @@ def far_to_nearfield(
     if not all(s == f.shape for f in [*E_far, *H_far]):
         raise Exception('All fields must be the same shape!')
 
-    if padded_size is None:
+    padded_shape = _padded_shape(s, padded_size)
-        padded_size = (2 ** numpy.ceil(numpy.log2(s))).astype(int)
-    if not hasattr(padded_size, '__len__'):
-        padded_size = (padded_size, padded_size)            # type: ignore  # checked if sequence
-    padded_shape = cast('Sequence[int]', padded_size)
 
     k = 2 * pi
     kxs = dkx * fftshift(fftfreq(s[0], d=1 / s[0]))
@@ -197,9 +221,7 @@ def far_to_nearfield(
     invalid_ind = kxy2 >= k * k
     theta[invalid_ind] = 0
     phi[invalid_ind] = 0
-    for i in range(2):
+    E_far = [numpy.where(invalid_ind, 0, field) for field in E_far]
-        E_far[i][invalid_ind] = 0
-        H_far[i][invalid_ind] = 0
 
     # Normalized vector potentials N, L
     L = [0.5 * E_far[1],
@@ -237,11 +259,11 @@ def far_to_nearfield(
             ),
         ]
 
-    E_near = [ifftshift(ifft2(ifftshift(Ei), s=padded_shape)) for Ei in En_fft]
+    dx = 2 * pi / (padded_shape[0] * dkx)
-    H_near = [ifftshift(ifft2(ifftshift(Hi), s=padded_shape)) for Hi in Hn_fft]
+    dy = 2 * pi / (padded_shape[1] * dky)
 
-    dx = 2 * pi / (s[0] * dkx)
+    E_near = [fftshift(ifft2(ifftshift(_resize_centered(Ei, padded_shape)))) / (dx * dy) for Ei in En_fft]
-    dy = 2 * pi / (s[1] * dky)
+    H_near = [fftshift(ifft2(ifftshift(_resize_centered(Hi, padded_shape)))) / (dx * dy) for Hi in Hn_fft]
 
     outputs = {
         'E': E_near,

meanas/test/test_fdfd_farfield.py

@@ -98,3 +98,55 @@ def test_farfield_roundtrip_supports_rectangular_arrays() -> None:
     assert restored['dy'] == pytest.approx(0.3)
     assert numpy.isfinite(restored['E'][0]).all()
     assert numpy.isfinite(restored['H'][0]).all()
+
+
+@pytest.mark.parametrize('shape', [(4, 4), (4, 5), (5, 4), (5, 7)])
+def test_farfield_padding_preserves_aperture_origin(shape) -> None:
+    rng = numpy.random.default_rng(7)
+    fields = [rng.normal(size=shape) + 1j * rng.normal(size=shape) for _ in range(4)]
+    target = (8, 10)
+    padding = [(new // 2 - old // 2, new - old - (new // 2 - old // 2)) for old, new in zip(shape, target, strict=True)]
+    padded = [numpy.pad(field, padding) for field in fields]
+    automatic = farfield.near_to_farfield(fields[:2], fields[2:], 0.5, 0.7, padded_size=target)
+    manual = farfield.near_to_farfield(padded[:2], padded[2:], 0.5, 0.7, padded_size=target)
+    for key in ('E', 'H'):
+        numpy.testing.assert_allclose(automatic[key], manual[key], atol=1e-13)
+
+
+@pytest.mark.parametrize('size', [4, 8])
+def test_farfield_amplitude_is_independent_of_aperture_sampling(size) -> None:
+    # A constant 2 by 2 wavelength aperture has integral 4 on either mesh.
+    ones = numpy.ones((size, size), dtype=complex)
+    zeros = numpy.zeros_like(ones)
+    result = farfield.near_to_farfield([ones, zeros], [zeros, ones], 2 / size, 2 / size, padded_size=size)
+    assert result['E'][0][size // 2, size // 2] == pytest.approx(8)
+
+
+@pytest.mark.parametrize('shape', [(4, 4), (16, 12), (7, 9)])
+def test_inverse_resizing_matches_plane_wave_quadrature(shape) -> None:
+    spectrum_shape = (8, 6)
+    dkx, dky = 0.4, 0.3
+    amplitude = 1.7 + 0.2j
+    e = [numpy.zeros(spectrum_shape, dtype=complex) for _ in range(2)]
+    e[0][spectrum_shape[0] // 2 + 1, spectrum_shape[1] // 2] = 2 * amplitude
+    h = [-e[1], e[0]]
+    result = farfield.far_to_nearfield(e, h, dkx, dky, padded_size=shape)
+    dx, dy = 2 * numpy.pi / (shape[0] * dkx), 2 * numpy.pi / (shape[1] * dky)
+    assert result['dx'] == pytest.approx(dx)
+    assert result['dy'] == pytest.approx(dy)
+    x = (numpy.arange(shape[0]) - shape[0] // 2) * dx
+    expected = numpy.broadcast_to(amplitude * dkx * dky / (2 * numpy.pi) ** 2 * numpy.exp(1j * dkx * x[:, None]), shape)
+    numpy.testing.assert_allclose(result['E'][0], expected, atol=1e-14)
+    numpy.testing.assert_allclose(result['H'][1], expected / numpy.sqrt(1 - (dkx / (2 * numpy.pi)) ** 2), atol=1e-14)
+
+
+def test_inverse_accepts_readonly_spectra_without_modifying_them() -> None:
+    rng = numpy.random.default_rng(9)
+    fields = [rng.normal(size=(8, 8)).astype(complex) for _ in range(4)]
+    originals = [field.copy() for field in fields]
+    for field in fields:
+        field.flags.writeable = False
+    result = farfield.far_to_nearfield(fields[:2], fields[2:], dkx=3, dky=3)
+    assert numpy.isfinite(result['E']).all()
+    for field, original in zip(fields, originals, strict=True):
+        numpy.testing.assert_array_equal(field, original)