# SPDX-License-Identifier: AGPL-3.0-or-later
# Copyright (C) 2025 SWGY, Inc
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
"""
Validation ladder for the BTM DC edge-diffraction coefficient (Stage 1).
The defensibility chain, each rung an independent cross-check:
(a) the exact static rigid-wedge Green's function (a Legendre-Q modal series)
reproduces exact image sums at wedge angles pi/n -> the oracle is trusted;
(b) the BTM coefficient in :mod:`diffraction` reproduces that oracle in the
deep-shadow zone (where the total field is purely diffracted), across wedge
angles including the screen (nu=1/2, the helmet-rim case);
(c) the diffracted field's jump across the shadow boundary equals the incident
direct field -- the continuity condition that also pins KAPPA=1/(16 pi^2).
These are pure NumPy/SciPy/mpmath (no Warp), so they gate the physics before any
engine work.
"""
import numpy as np
import mpmath as mp
import pytest
import diffraction
mp.mp.dps = 30
# --- The exact oracle: static rigid-wedge Green's function -------------------
def _legenq(deg, x):
# Q_deg(x), x>1, real degree -> real; mpmath returns mpc on this branch.
return float(mp.legenq(deg, 0, x, type=3).real)
def g_wedge(rho, theta, z, rho0, theta0, z0, theta_w, mmax=2000):
"""Exact Laplace Green's function for a Neumann wedge (fluid angle theta_w).
Modal series with the k-integral done in closed form
(int cos(kZ) I_mu K_mu dk = Q_{mu-1/2}(chi)/(2 sqrt(rho rho0))).
"""
Z = z - z0
chi = (rho**2 + rho0**2 + Z**2) / (2.0 * rho * rho0)
nu = np.pi / theta_w
s, small = 0.0, 0
for m in range(mmax + 1):
eps = 1.0 if m == 0 else 2.0
mu = m * nu
q = _legenq(mu - 0.5, chi)
s += eps * np.cos(mu * theta) * np.cos(mu * theta0) * q
small = small + 1 if eps * abs(q) < 1e-15 * max(abs(s), 1e-300) else 0
if small >= 8:
break
return s / (2.0 * np.pi * theta_w * np.sqrt(rho * rho0))
def g_images(rho, theta, z, rho0, theta0, z0, theta_w, n):
"""Exact field for theta_w = pi/n via 2n rigid (Neumann) images."""
Z = z - z0
total = 0.0
for m in range(n):
for sgn in (+1.0, -1.0):
alpha = 2.0 * m * theta_w + sgn * theta0
d2 = rho**2 + rho0**2 - 2*rho*rho0*np.cos(theta - alpha) + Z**2
total += 1.0 / (4.0 * np.pi * np.sqrt(d2))
return total
# --- (a) the oracle itself is exact ------------------------------------------
@pytest.mark.parametrize("n", [1, 2, 3, 4, 6])
def test_oracle_matches_image_sums(n):
tw = np.pi / n
rho, z, rho0, z0 = 1.15, 0.3, 0.8, -0.1
theta, theta0 = 0.35 * tw, 0.72 * tw
gw = g_wedge(rho, theta, z, rho0, theta0, z0, tw)
gi = g_images(rho, theta, z, rho0, theta0, z0, tw, n)
assert abs(gw - gi) <= 1e-9 * abs(gi), \
"modal wedge oracle disagrees with the exact image sum"
# --- (b) BTM reproduces the exact oracle in the shadow ------------------------
# (theta_w, theta_S, deep-shadow theta_R fractions of theta_w)
_SHADOW_CASES = [
(1.5 * np.pi, 0.15 * np.pi, (1.25, 1.30, 1.35, 1.40)), # 270 deg wedge
(2.0 * np.pi, 0.15 * np.pi, (1.30, 1.50, 1.70, 1.90)), # screen (rim case)
]
@pytest.mark.parametrize("theta_w,theta_S,thR_fracs", _SHADOW_CASES)
def test_btm_matches_oracle_in_shadow(theta_w, theta_S, thR_fracs):
r_S = 1.0
for f in thR_fracs:
theta_R = f * np.pi
for r_R, z_R in ((1.7, 0.0), (0.6, 0.5), (1.3, -0.9), (2.2, 0.4)):
btm = diffraction.diffracted_dc_field(
theta_w, r_S, theta_S, 0.0, r_R, theta_R, z_R)
exact = g_wedge(r_R, theta_R, z_R, r_S, theta_S, 0.0, theta_w)
assert abs(btm - exact) <= 2e-3 * abs(exact), (
f"BTM {btm:.6e} vs exact {exact:.6e} at "
f"theta_w={np.degrees(theta_w):.0f} theta_R={np.degrees(theta_R):.0f}")
# --- (c) shadow-boundary continuity (jump = incident, pins KAPPA) -------------
def test_shadow_boundary_jump_equals_incident():
# screen; incident shadow boundary at theta_R = theta_S + pi
theta_w, theta_S, r_S = 2.0 * np.pi, np.radians(50.0), 1.0
r_R, z_R = 1.4, 0.0
sb = theta_S + np.pi
eps = np.radians(0.4)
lit = diffraction.diffracted_dc_field(theta_w, r_S, theta_S, 0.0, r_R, sb - eps, z_R)
sha = diffraction.diffracted_dc_field(theta_w, r_S, theta_S, 0.0, r_R, sb + eps, z_R)
# incident direct field of a unit Laplace source at the boundary
S = np.array([r_S, 0.0, 0.0])
R = np.array([r_R * np.cos(np.pi), r_R * np.sin(np.pi), z_R]) # relative angle sb-theta_S=pi
R_SR = np.linalg.norm(S - R)
incident = 1.0 / (4.0 * np.pi * R_SR)
assert abs((sha - lit) - incident) <= 1e-2 * incident, \
"diffracted jump across the shadow boundary must equal the incident field"
# --- per-segment quadrature (kernel path) matches the reference integral -------
def test_segment_sum_matches_integral():
"""Midpoint per-segment sum -> the adaptive edge integral (the kernel path)."""
theta_w, nu = 2.0 * np.pi, 0.5
r_S, theta_S = 1.0, np.radians(40.0)
r_R, theta_R, z_R = 1.3, np.radians(300.0), 0.2 # deep shadow, off-apex
ref = diffraction.diffracted_dc_field(theta_w, r_S, theta_S, 0.0, r_R, theta_R, z_R)
zs = np.linspace(-40, 40, 40001)
dz = zs[1] - zs[0]
total = 0.0
for z in zs:
total += diffraction.segment_dc_gain(
nu, r_S, theta_S, 0.0 - z, r_R, theta_R, z_R - z, dz)
assert abs(total - ref) <= 2e-3 * abs(ref), \
"per-segment quadrature disagrees with the reference edge integral"
# --- (Stage 2) the Warp kernel reproduces the engine formula in NumPy ---------
def _engine_diff_numpy(edges, S, R, s_fan, p0, k):
"""Independent NumPy evaluation of the _diffract per-segment sum."""
R_SR = np.linalg.norm(S - R)
total = 0.0
for i in range(edges.apex.shape[0]):
rS, thS, zS = diffraction.to_edge_frame(S, edges.apex[i], edges.tangent[i],
edges.n_ref[i])
rR, thR, zR = diffraction.to_edge_frame(R, edges.apex[i], edges.tangent[i],
edges.n_ref[i])
if rS < 1e-6 or rR < 1e-6:
continue
m, l = np.hypot(rS, zS), np.hypot(rR, zR)
cosh_eta = (m * l + zS * zR) / (rS * rR)
beta = diffraction._beta_sum(edges.nu[i], thS, thR, cosh_eta)
total += ((-edges.nu[i] / (4.0 * np.pi)) * s_fan * (p0 * k / R_SR)
* (beta / (m * l)) * edges.seg_len[i])
return total
def test_warp_diffract_matches_numpy():
wp = pytest.importorskip("warp")
wp.config.cpu_compiler_flags = "-mcpu=x86-64" # container CPU misdetection
import warp_core
device = warp_core.configure("cpu")
# a synthetic straight edge (along x), off to one side; sensor at origin
n = 60
xs = np.linspace(-0.5, 0.5, n)
seg = xs[1] - xs[0]
apex = np.stack([xs, np.full(n, 0.4), np.full(n, -0.3)], axis=1)
tangent = np.tile([1.0, 0.0, 0.0], (n, 1))
n_ref = np.tile([0.0, 1.0, 0.0], (n, 1))
edges = diffraction.EdgeArrays(apex, tangent, n_ref,
np.full(n, 0.5), np.full(n, seg))
S = np.array([0.2, -1.5, 0.1])
R = np.zeros(3)
s_fan, p0, k = 3.2e-4, 50000.0, 3.679e-3
ref = _engine_diff_numpy(edges, S, R, s_fan, p0, k)
out = wp.zeros(1, dtype=wp.float64, device=device)
wp.launch(warp_core._diffract, dim=(1, n), device=device, inputs=[
wp.array(S.reshape(1, 3).astype(np.float32), dtype=wp.vec3, device=device),
wp.array(apex.astype(np.float32), dtype=wp.vec3, device=device),
wp.array(tangent.astype(np.float32), dtype=wp.vec3, device=device),
wp.array(n_ref.astype(np.float32), dtype=wp.vec3, device=device),
wp.array(np.full(n, 0.5, np.float32), dtype=wp.float32, device=device),
wp.array(np.full(n, seg, np.float32), dtype=wp.float32, device=device),
wp.array(np.array([s_fan]), dtype=wp.float64, device=device),
wp.uint64(0), wp.uint64(0), 0, # no meshes -> no occlusion
wp.vec3(0.0, 0.0, 0.0), float(p0), float(k), float(warp_core._T_MIN), out,
])
wp.synchronize()
got = float(out.numpy()[0])
assert abs(got - ref) <= 1e-4 * abs(ref) + 1e-12, \
f"Warp _diffract {got:.6e} != NumPy reference {ref:.6e}"
# --- (Stage 3) end-to-end: diffraction on the real helmet ---------------------
_E2E_OVERRIDES = [
"blast.cube_segments=5", "blast.cube_extents_m=12",
"blast.cube_center_m=[0,0,0]", "blast.min_dist_m=1.0",
"rays.resolution_arcmin=25", "rays.max_bounces=1",
]
def test_diffraction_end_to_end():
wp = pytest.importorskip("warp")
wp.config.cpu_compiler_flags = "-mcpu=x86-64"
import config as config_mod
import simulate_blasts
off = np.array([r[3:] for r in simulate_blasts.run(
config_mod.load_config("config.yaml", _E2E_OVERRIDES))])
on = np.array([r[3:] for r in simulate_blasts.run(
config_mod.load_config("config.yaml", _E2E_OVERRIDES + [
"diffraction.enabled=true", "diffraction.segment_length_m=0.03"]))])
assert np.isfinite(on).all() and (on >= 0.0).all(), \
"diffracted totals must be finite and clamped nonnegative"
assert np.any(on != off), "enabling diffraction must change some results"
# shadow fill: some cell the geometric model nulled becomes nonzero
helmet_j = 2 # CONFIGS index of helmet_only
assert np.any((off[:, helmet_j] == 0.0) & (on[:, helmet_j] > 0.0)), \
"edge diffraction should fill at least one geometric-shadow cell"