diff --git a/api/app/core/__pycache__/diffraction.cpython-313.pyc b/api/app/core/__pycache__/diffraction.cpython-313.pyc index 83212c8..fba477b 100644 Binary files a/api/app/core/__pycache__/diffraction.cpython-313.pyc and b/api/app/core/__pycache__/diffraction.cpython-313.pyc differ diff --git a/api/app/core/diffraction.py b/api/app/core/diffraction.py index 2f97fd3..d0e314a 100644 --- a/api/app/core/diffraction.py +++ b/api/app/core/diffraction.py @@ -19,53 +19,43 @@ def knife_edge_v(h: float, d1: float, d2: float, freq_hz: float) -> float: return h * sqrt(2 * (d1 + d2) / (lmbda * d1 * d2)) +def bullington_loss( + profile: SurfaceProfile, + tx_height_agl: float, + rx_height_agl: float, + freq_hz: float, +) -> float: + """Bullington-style equivalent edge loss for a terrain profile. + + The current implementation uses the dominant obstacle relative to the TX-RX + chord as the equivalent Bullington edge, then applies ITU-R P.526 J(v). + """ + if len(profile.samples) < 3: + return 0.0 + + total_distance = profile.distance_m + if total_distance <= 0: + return 0.0 + + tx_elevation = profile.samples[0].ground_m + tx_height_agl + rx_elevation = profile.samples[-1].ground_m + rx_height_agl + max_v = float("-inf") + + for sample in profile.samples[1:-1]: + d1 = sample.distance_m + d2 = total_distance - d1 + path_height = tx_elevation + (rx_elevation - tx_elevation) * (d1 / total_distance) + h = sample.surface_m - path_height + max_v = max(max_v, knife_edge_v(h, d1, d2, freq_hz)) + + return knife_edge_loss(max_v) + + def deygout( profile: SurfaceProfile, tx_height_agl: float, rx_height_agl: float, freq_hz: float, ) -> float: - """Recursive Deygout diffraction loss using the dominant edge and subprofiles.""" - if len(profile.samples) < 3: - return 0.0 - - endpoint_heights = { - 0: profile.samples[0].ground_m + tx_height_agl, - len(profile.samples) - 1: profile.samples[-1].ground_m + rx_height_agl, - } - - def sample_height(index: int) -> float: - return endpoint_heights.get(index, profile.samples[index].surface_m) - - def solve(left: int, right: int) -> float: - if right - left < 2: - return 0.0 - - left_sample = profile.samples[left] - right_sample = profile.samples[right] - span_m = right_sample.distance_m - left_sample.distance_m - left_height = sample_height(left) - right_height = sample_height(right) - max_v = float("-inf") - max_index: int | None = None - - for index in range(left + 1, right): - sample = profile.samples[index] - d1 = sample.distance_m - left_sample.distance_m - d2 = right_sample.distance_m - sample.distance_m - path_height = left_height + (right_height - left_height) * (d1 / span_m) - h = sample.surface_m - path_height - v = knife_edge_v(h, d1, d2, freq_hz) - if v > max_v: - max_v = v - max_index = index - - if max_index is None: - return 0.0 - - main_loss = knife_edge_loss(max_v) - if main_loss == 0.0: - return 0.0 - return main_loss + solve(left, max_index) + solve(max_index, right) - - return solve(0, len(profile.samples) - 1) + """Compatibility wrapper; use Bullington equivalent loss for multi-edge profiles.""" + return bullington_loss(profile, tx_height_agl, rx_height_agl, freq_hz) diff --git a/api/tests/__pycache__/test_core.cpython-313-pytest-9.0.3.pyc b/api/tests/__pycache__/test_core.cpython-313-pytest-9.0.3.pyc index 9aa5319..b4912c4 100644 Binary files a/api/tests/__pycache__/test_core.cpython-313-pytest-9.0.3.pyc and b/api/tests/__pycache__/test_core.cpython-313-pytest-9.0.3.pyc differ diff --git a/api/tests/test_core.py b/api/tests/test_core.py index 3dff795..a97ebe8 100644 --- a/api/tests/test_core.py +++ b/api/tests/test_core.py @@ -3,9 +3,10 @@ from math import isclose import pytest from app.core.antenna import AntennaPattern, gain -from app.core.diffraction import knife_edge_loss +from app.core.diffraction import bullington_loss, knife_edge_loss, knife_edge_v from app.core.fresnel import earth_bulge, fresnel_radius, wavelength from app.core.propagation import fspl +from app.core.surface import SurfaceProfile, SurfaceSample from app.core.vegetation import P833Coefficients, p833_attenuation @@ -34,6 +35,66 @@ def test_knife_edge_loss_reference_values(v: float, expected: float) -> None: assert isclose(knife_edge_loss(v), expected, abs_tol=0.05) +def test_bullington_loss_uses_dominant_equivalent_edge() -> None: + profile = SurfaceProfile( + distance_m=10_000, + samples=[ + SurfaceSample( + i=0, + lat=0, + lon=0, + distance_m=0, + ground_m=0, + building_m=0, + canopy_m=0, + surface_m=0, + ), + SurfaceSample( + i=1, + lat=0, + lon=0, + distance_m=5_000, + ground_m=20, + building_m=0, + canopy_m=0, + surface_m=20, + ), + SurfaceSample( + i=2, + lat=0, + lon=0, + distance_m=10_000, + ground_m=0, + building_m=0, + canopy_m=0, + surface_m=0, + ), + ], + ) + v = knife_edge_v(20, 5_000, 5_000, 433_000_000) + + assert isclose(bullington_loss(profile, 0, 0, 433_000_000), knife_edge_loss(v)) + + +def test_bullington_loss_does_not_sum_every_dem_sample() -> None: + samples = [ + SurfaceSample( + i=i, + lat=0, + lon=0, + distance_m=i * 100, + ground_m=5 if 20 <= i <= 80 else 0, + building_m=0, + canopy_m=0, + surface_m=5 if 20 <= i <= 80 else 0, + ) + for i in range(101) + ] + profile = SurfaceProfile(distance_m=10_000, samples=samples) + + assert bullington_loss(profile, 0, 0, 433_000_000) < 30 + + def test_p833_formula_with_supplied_coefficients() -> None: attenuation = p833_attenuation( depth_m=100,