[FEATURE] Implement UTD `diffraction_coefficients`

[rwydaegh]

Terms

• Checked the existing issues to see if my suggestion has not already been suggested;

Description

diffraction_coefficients() in differt/em/_utd.py raises NotImplementedError, but almost everything is already there. The helpers _cot, _N, _a, L_i, and the transition function F are all implemented and tested. There's a dead-code draft of the actual computation after the raise (lines 260-303). The docstring literally asks for help:

This function is not yet implemented, as we are still thinking of the best API for it. If you want to get involved in the implementation of UTD coefficients, please reach out to us on GitHub!

I'd go with something like:

@jax.jit
def diffraction_coefficients(
    k: Float[ArrayLike, " *#batch"],           # wavenumber
    n: Float[ArrayLike, " *#batch"],           # exterior wedge parameter (angle / pi)
    phi_i: Float[ArrayLike, " *#batch"],       # incidence angle, edge-fixed coords
    phi_d: Float[ArrayLike, " *#batch"],       # diffraction angle, edge-fixed coords
    sin_beta_0: Float[ArrayLike, " *#batch"],  # sin(angle between incident ray and edge)
    L_i: Float[ArrayLike, " *#batch"],         # distance parameter (incident SB)
    L_r_n: Float[ArrayLike, " *#batch"] | None = None,
    L_r_o: Float[ArrayLike, " *#batch"] | None = None,
) -> tuple[Complex[Array, " *batch"], Complex[Array, " *batch"]]:
    """UTD diffraction coefficients D_s, D_h for a PEC wedge.
    McNamara et al. (1990), eq. 6.21-6.28. L_r_n/L_r_o default to L_i if not given.
    """

And a convenience function compute_diffraction_angles that extracts phi_i, phi_d, sin_beta_0 from 3D geometry (edge endpoints, source, observation point, face normals). Keeps the coefficient function clean by accepting pre-computed angles.

There's a bug in the existing _N helper

Lines 26-31 have identical implementations for "+" and "-" modes:

def _N(beta, n, mode):
    if mode == "+":
        return jnp.round((beta + jnp.pi) / (2 * n * jnp.pi))
    return jnp.round((beta + jnp.pi) / (2 * n * jnp.pi))  # same formula!

The "-" branch should use beta - jnp.pi based on the existing docstrings in _utd.py which reference McNamara et al. (1990) Introduction to the Uniform Geometrical Theory of Diffraction, eq. 6.22.

Scope

Fairly small. The math exists in dead code already. Clean up the API, fix the _N bug, remove the raise, test against McNamara's tables. PEC only for now. Lossy materials need Maliuzhinets functions which is a much bigger undertaking.

For validation: the existing docstrings reference McNamara et al. Table 6.1 for known coefficient values, so that's a natural first check. Also continuity across shadow boundaries, and knife-edge (n=2) against the closed-form Fresnel integral.

This needs #428 (edge detection) first for wedge_angles and edge geometry. Once both are in, #431 (transition_matrices) can add a diffraction branch that calls this.

[jeertmans]

Hi @rwydaegh! Indeed, the placeholder is very likely incorrect :-) It would be very nice to have proper tests (as you said, McNamara's book includes many examples that we can use), I agree!