[FEATURE] Implement UTD `diffraction_coefficients`
enhancementnice-to-have
## Terms
- [x] Checked the [existing issues](https://github.com/jeertmans/DiffeRT/issues?q=is%3Aissue+label%3Aenhancement) 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:
```python
@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:
```python
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.
1 条评论