From be574d1a5bd453f2c2e6308758955e8422fe79c7 Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Wed, 28 Oct 2026 18:14:22 -0700 Subject: [PATCH] [farfield] preserve centered geometry and physical Fourier normalization --- meanas/fdfd/farfield.py | 78 ++++++++++++++++++++----------- meanas/test/test_fdfd_farfield.py | 52 +++++++++++++++++++++ 2 files changed, 102 insertions(+), 28 deletions(-) diff --git a/meanas/fdfd/farfield.py b/meanas/fdfd/farfield.py index 06f705b..2fece74 100644 --- a/meanas/fdfd/farfield.py +++ b/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_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) + padded_shape = _padded_shape(s, padded_size) - En_fft = [fftshift(fft2(fftshift(Eni), s=padded_shape)) for Eni in E_near] - Hn_fft = [fftshift(fft2(fftshift(Hni), s=padded_shape)) for Hni in H_near] + En_fft = [fftshift(fft2(ifftshift(_resize_centered(Eni, padded_shape)))) * dx * dy for Eni in E_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]), - 'dky': float(kyy[1] - kyy[0]), + 'dkx': float(2 * pi / (padded_shape[0] * dx)), + '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 - through several wavelengths of uniform medium. + Reconstruct the propagating near-field spectrum from a normalized farfield. - 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_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) + padded_shape = _padded_shape(s, 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[i][invalid_ind] = 0 - H_far[i][invalid_ind] = 0 + E_far = [numpy.where(invalid_ind, 0, field) for field in E_far] # 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] - H_near = [ifftshift(ifft2(ifftshift(Hi), s=padded_shape)) for Hi in Hn_fft] + dx = 2 * pi / (padded_shape[0] * dkx) + dy = 2 * pi / (padded_shape[1] * dky) - dx = 2 * pi / (s[0] * dkx) - dy = 2 * pi / (s[1] * dky) + E_near = [fftshift(ifft2(ifftshift(_resize_centered(Ei, padded_shape)))) / (dx * dy) for Ei in En_fft] + H_near = [fftshift(ifft2(ifftshift(_resize_centered(Hi, padded_shape)))) / (dx * dy) for Hi in Hn_fft] outputs = { 'E': E_near, diff --git a/meanas/test/test_fdfd_farfield.py b/meanas/test/test_fdfd_farfield.py index 5e2daab..2d6c0d3 100644 --- a/meanas/test/test_fdfd_farfield.py +++ b/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)