From 7b44eb0c2e4bad5f1b2bf185f65d3e189135291a Mon Sep 17 00:00:00 2001 From: Philipp Wagner Date: Sun, 11 Oct 2026 18:39:48 +0200 Subject: [PATCH] Reference cases for the path-analysis mathematics (#108) (#117) * test: pin the Maidenhead and great-circle maths of Location against definitions (#108) Subsquare centres from the Maidenhead definition, round trip through toMaidenhead (lower-case subsquare pinned as is), distance and bearing between locators against an independent spherical computation on R = 6371.009 km, consistency with PathGeometryUtils, and the NaN bearing for identical points pinned as current behaviour. Co-Authored-By: Claude Fable 5.1 * test: pin PathGeometryUtils against ITU-R P.525/P.526 and textbook values (#108) Great circle against the Veness worked example, Fresnel radius against P.526-15 eq. (2) and (3), earth bulge against the 8500 km effective radius, radio horizon against 4.12*sqrt(h), knife-edge loss J(v) against eq. (31), the v parameter against eq. (26), free-space loss pinned at the project's 32.44 constant with the P.525 and exact values documented, plus wattsToDbm, the sample-count rule and the k-factor sanitiser. Co-Authored-By: Claude Fable 5.1 * test: pin the path-analysis pipeline on a flat profile (#108) A flat 20 km profile with 10 m antennas at 144.3 MHz through an in-memory terrain provider: bulge, line-of-sight clearance, first Fresnel radius and intrusion, no knife edge and 0 dB diffraction. Pins that a partial Fresnel intrusion without a blocked line of sight contributes 0 dB to the link budget; P.526 would give 0 to 6 dB there, which is tracked separately. Fixes #108 Co-Authored-By: Claude Fable 5.1 --------- Co-authored-by: Claude Fable 5.1 --- .../LocationMaidenheadReferenceTest.java | 168 ++++++++++++ ...tryOnlyPathAnalysisServiceProfileTest.java | 197 ++++++++++++++ .../map/PathGeometryUtilsReferenceTest.java | 253 ++++++++++++++++++ 3 files changed, 618 insertions(+) create mode 100644 core/src/test/java/kst4contest/locatorUtils/LocationMaidenheadReferenceTest.java create mode 100644 core/src/test/java/kst4contest/view/map/GeometryOnlyPathAnalysisServiceProfileTest.java create mode 100644 core/src/test/java/kst4contest/view/map/PathGeometryUtilsReferenceTest.java diff --git a/core/src/test/java/kst4contest/locatorUtils/LocationMaidenheadReferenceTest.java b/core/src/test/java/kst4contest/locatorUtils/LocationMaidenheadReferenceTest.java new file mode 100644 index 00000000..88f1e1a1 --- /dev/null +++ b/core/src/test/java/kst4contest/locatorUtils/LocationMaidenheadReferenceTest.java @@ -0,0 +1,168 @@ +package kst4contest.locatorUtils; + +import kst4contest.view.map.PathGeometryUtils; +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertTrue; + +/** + * Reference cases for the Maidenhead and great-circle maths in {@link Location}. + * + *

None of the expected values below was derived from the code under test. The + * Maidenhead centres follow from the definition of the locator system, the distances + * and bearings from an independent spherical computation (haversine and atan2, see + * {@link #expectedDistanceAndBearingOfLocatorPairs}), both on the IUGG mean radius + * R = 6371.009 km that {@code Location} declares. Behaviour that is merely pinned + * (not endorsed) is marked as such in the test method's Javadoc.

+ * + *

Context: issue #108.

+ */ +class LocationMaidenheadReferenceTest { + + /** Tolerance for a subsquare centre: well below the 1e-7 degree print precision. */ + private static final double CENTRE_TOLERANCE_DEG = 1e-6; + + /** Tolerance for values that the code rounds to 0.01 km / 0.01 degrees. */ + private static final double ROUNDED_TOLERANCE = 0.005; + + /** + * Locator to subsquare centre, by definition. + * + *

Source: Wikipedia, "Maidenhead Locator System". A field spans 20 degrees of + * longitude and 10 degrees of latitude, a square 2 by 1 degrees, a subsquare 5 by + * 2.5 arc minutes. The centre of a subsquare is therefore

+ * 
+     *   lon = -180 + 20 * F1 + 2 * S1 + (5/60) * X1 + (5/60)/2
+     *   lat =  -90 + 10 * F2 +     S2 + (2.5/60) * X2 + (2.5/60)/2
+     *
+ *

with F = field letter index (A = 0), S = square digit, X = subsquare letter + * index. Tolerance: 1e-6 degrees.

+ */ + @Test + void locatorIsResolvedToTheCentreOfItsSubsquare() { + assertCentre("JO62QM", 52.5208333, 13.3750000); + assertCentre("JN49FK", 49.4375000, 8.4583333); + assertCentre("JO50JP", 50.6458333, 10.7916667); + assertCentre("AA00AA", -89.9791667, -179.9583333); + assertCentre("RR99XX", 89.9791667, 179.9583333); + } + + /** + * Round trip from the subsquare centre back to the locator. + * + *

The coordinates are the JO62QM centre of the previous test. The expected string + * uses a lower-case subsquare: that is what {@link Location#toMaidenhead(double, double)} + * produces today and is pinned here, not endorsed. The locator system itself is + * case-insensitive; the ON4KST convention is upper case.

+ */ + @Test + void centreCoordinatesRoundTripToTheLocatorWithLowerCaseSubsquare() { + assertEquals("JO62qm", Location.toMaidenhead(52.5208333, 13.375)); + } + + /** + * Distance and bearing between locator centres against an independent spherical + * computation. + * + *

Expected values were computed with Python 3 ({@code python3 -I}) on the subsquare + * centres of the previous test, R = 6371.009 km, using

+ * 
+     *   a     = sin^2(dphi/2) + cos(phi1) * cos(phi2) * sin^2(dlambda/2)
+     *   d     = 2 * R * atan2(sqrt(a), sqrt(1 - a))                        (haversine)
+     *   theta = atan2(sin(dlambda) * cos(phi2),
+     *                 cos(phi1) * sin(phi2) - sin(phi1) * cos(phi2) * cos(dlambda))
+     *   bearing = (deg(theta) + 360) mod 360
+     *
+ *

giving JO50JP to JO31AA 336.1207 km / 278.963 deg, JO50JP to JN49FK + * 214.0343 km / 232.018 deg, JO62QM to JN58TD 502.0447 km / 195.002 deg. + * {@link Location#getDistanceKmByTwoLocatorStrings} and {@link Location#getBearing} + * round to two decimals, so the rounded values are pinned with a tolerance of + * 0.005 km and 0.005 degrees.

+ */ + @Test + void expectedDistanceAndBearingOfLocatorPairs() { + assertDistanceAndBearing("JO50JP", "JO31AA", 336.12, 278.96); + assertDistanceAndBearing("JO50JP", "JN49FK", 214.03, 232.02); + assertDistanceAndBearing("JO62QM", "JN58TD", 502.04, 195.00); + } + + /** + * The two distance implementations in the project agree. + * + *

{@link Location} uses the spherical law of cosines, {@link PathGeometryUtils} + * the haversine formula; both declare R = 6371.009 km. On the same pair of subsquare + * centres the unrounded results must agree within 0.005 km (the law of cosines loses + * precision only for distances far below one kilometre).

+ */ + @Test + void lawOfCosinesAndHaversineAgreeOnTheSamePair() { + Location from = new Location("JO50JP"); + Location to = new Location("JN58TD"); + + double lawOfCosinesKm = from.getDistanceKm(to); + double haversineKm = PathGeometryUtils.calculateGreatCircleDistanceKm( + from.getLatitude().toDegrees(), + from.getLongitude().toDegrees(), + to.getLatitude().toDegrees(), + to.getLongitude().toDegrees() + ); + + assertEquals(haversineKm, lawOfCosinesKm, ROUNDED_TOLERANCE, + "law of cosines vs haversine on JO50JP -> JN58TD"); + } + + /** + * Bearing between identical points in {@link Location}: NaN. + * + *

Pinned, not endorsed. The bearing between a point and itself is undefined; + * {@code Location.getBearing} answers {@code Double.NaN}.

+ */ + @Test + void locationBearingOfIdenticalPointsIsNaN() { + Location point = new Location("JO62QM"); + + assertTrue(Double.isNaN(Location.getBearing(point, new Location("JO62QM"))), + "Location.getBearing for identical points"); + } + + /** + * Bearing between identical points in {@link PathGeometryUtils}: 0 degrees. + * + *

Pinned, not endorsed. For identical coordinates both atan2 arguments are zero + * and {@code atan2(0, 0)} is 0, so {@code calculateInitialBearingDeg} answers 0 where + * {@code Location.getBearing} answers NaN. The two implementations disagree on the + * undefined case; this test records the current answer.

+ */ + @Test + void pathGeometryBearingOfIdenticalPointsIsZero() { + double bearingDeg = PathGeometryUtils.calculateInitialBearingDeg( + 52.5208333, 13.375, 52.5208333, 13.375); + + assertEquals(0.0, bearingDeg, 1e-12, "PathGeometryUtils bearing for identical points"); + } + + private static void assertCentre(String locator, double expectedLatDeg, double expectedLonDeg) { + Location location = new Location(locator); + assertEquals(expectedLatDeg, location.getLatitude().toDegrees(), CENTRE_TOLERANCE_DEG, + locator + " latitude"); + assertEquals(expectedLonDeg, location.getLongitude().toDegrees(), CENTRE_TOLERANCE_DEG, + locator + " longitude"); + } + + private static void assertDistanceAndBearing(String fromLocator, + String toLocator, + double expectedRoundedKm, + double expectedRoundedBearingDeg) { + Location helper = new Location(); + + assertEquals(expectedRoundedKm, + helper.getDistanceKmByTwoLocatorStrings(fromLocator, toLocator), + ROUNDED_TOLERANCE, + fromLocator + " -> " + toLocator + " distance"); + assertEquals(expectedRoundedBearingDeg, + helper.getBearingOfTwoLocatorStrings(fromLocator, toLocator), + ROUNDED_TOLERANCE, + fromLocator + " -> " + toLocator + " bearing"); + } +} diff --git a/core/src/test/java/kst4contest/view/map/GeometryOnlyPathAnalysisServiceProfileTest.java b/core/src/test/java/kst4contest/view/map/GeometryOnlyPathAnalysisServiceProfileTest.java new file mode 100644 index 00000000..7f1f3394 --- /dev/null +++ b/core/src/test/java/kst4contest/view/map/GeometryOnlyPathAnalysisServiceProfileTest.java @@ -0,0 +1,197 @@ +package kst4contest.view.map; + +import org.junit.jupiter.api.Test; + +import java.util.ArrayList; +import java.util.List; + +import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertFalse; +import static org.junit.jupiter.api.Assertions.assertTrue; + +/** + * The path-analysis pipeline of {@link GeometryOnlyPathAnalysisService} on a flat + * terrain profile, with every expected value computed from the underlying formulas + * rather than from the code. + * + *

Scenario: 20.000 km due north from 50.0 N 10.0 E, both antennas 10 m above a + * terrain of 0 m everywhere, 144.3 MHz, k = 4/3, default link-budget settings. A path + * along a meridian makes the haversine distance exactly R times the latitude difference, + * so the 20 km are hit to floating-point precision. The in-memory terrain provider hands + * back the requested number of samples (121 for 20 km, so the midpoint is sample 60 at + * exactly 10.000 km) with elevation 0.

+ * + *

Expected midpoint values (R = 6 371 009 m, c = 299 792 458 m/s):

+ * 
+ *   bulge      = d1 d2 / (2 k R)        = 1e4 * 1e4 / (2 * 4/3 * 6371009) = 5.886 m
+ *   LOS height = 10 m (both endpoints have bulge 0 and antenna 10 m)
+ *   clearance  = 10 - 5.886                                              = 4.114 m
+ *   lambda     = c / 144.3e6                                             = 2.0776 m
+ *   R1         = sqrt(lambda d1 d2 / (d1 + d2)) = sqrt(2.0776 * 5000)   = 101.921 m
+ *   intrusion  = bulge - (LOS - R1) = 5.886 - (10 - 101.921)            = 97.807 m
+ *
+ *

Tolerances: 0.05 m for bulge and clearance, 0.1 m for Fresnel radius and intrusion, + * 0.01 dB for losses.

+ * + *

Context: issue #108.

+ */ +class GeometryOnlyPathAnalysisServiceProfileTest { + + private static final double EARTH_RADIUS_KM = 6371.009; + private static final double TOTAL_DISTANCE_KM = 20.0; + private static final double ANTENNA_HEIGHT_M = 10.0; + private static final double FREQUENCY_MHZ = 144.3; + private static final int MIDPOINT_INDEX = 60; + + private static final double EXPECTED_MIDPOINT_BULGE_M = 5.886; + private static final double EXPECTED_MIDPOINT_LOS_CLEARANCE_M = 4.114; + private static final double EXPECTED_MIDPOINT_FRESNEL_RADIUS_M = 101.921; + private static final double EXPECTED_MIDPOINT_INTRUSION_M = 97.807; + + /** + * Terrain provider that answers a flat profile of the requested sample count. + * Points are spaced uniformly along the path; latitude is interpolated linearly, + * which is exact for a meridian. + */ + private static final class FlatTerrainProfileProvider implements TerrainProfileProvider { + @Override + public TerrainProfileData loadProfile(TerrainProfileRequest request) { + int sampleCount = request.requestedSampleCount(); + List points = new ArrayList<>(sampleCount); + for (int i = 0; i < sampleCount; i++) { + double fraction = (double) i / (sampleCount - 1); + double distanceKm = request.totalDistanceKm() * fraction; + double latitudeDeg = request.fromLatitudeDeg() + + (request.toLatitudeDeg() - request.fromLatitudeDeg()) * fraction; + points.add(new PathProfilePoint(distanceKm, latitudeDeg, request.fromLongitudeDeg(), 0.0)); + } + return new TerrainProfileData(points, "Flat test terrain", false); + } + } + + private static PathAnalysisResult analyzeFlatPath() { + double fromLatitudeDeg = 50.0; + double longitudeDeg = 10.0; + double toLatitudeDeg = fromLatitudeDeg + Math.toDegrees(TOTAL_DISTANCE_KM / EARTH_RADIUS_KM); + + PathAnalysisRequest request = new PathAnalysisRequest( + "JO50AA", + fromLatitudeDeg, + longitudeDeg, + "DL0TEST", + "JO50AE", + toLatitudeDeg, + longitudeDeg, + FREQUENCY_MHZ, + ANTENNA_HEIGHT_M, + ANTENNA_HEIGHT_M, + PathGeometryUtils.DEFAULT_EFFECTIVE_EARTH_RADIUS_FACTOR, + PathLinkBudgetSettings.defaults() + ); + return new GeometryOnlyPathAnalysisService(new FlatTerrainProfileProvider()).analyze(request); + } + + /** + * The scenario is what the Javadoc says: 20.000 km, 121 samples, midpoint at 10 km. + * + *

Sample count: ceil(20 / 0.5) + 1 = 41, clamped to the minimum of 121. + * Tolerance: 1e-6 km for the distance.

+ */ + @Test + void scenarioIsTwentyKilometresWithTheMidpointAtSampleSixty() { + PathAnalysisResult result = analyzeFlatPath(); + + assertTrue(result.hasUsableProfile(), "profile usable"); + assertEquals(TOTAL_DISTANCE_KM, result.distanceKm(), 1e-6, "path distance"); + assertEquals(121, result.profilePoints().size(), "sample count"); + assertEquals(10.0, result.profilePoints().get(MIDPOINT_INDEX).distanceKm(), 1e-9, "midpoint distance"); + } + + /** + * Earth bulge at the midpoint: d1 d2 / (2 k R) = 5.886 m at k = 4/3. + * + *

On a flat profile the curvature-adjusted elevation is the bulge itself. + * Tolerance: 0.05 m.

+ */ + @Test + void midpointBulgeIsFivePointNineMetres() { + PathProfilePoint midpoint = analyzeFlatPath().profilePoints().get(MIDPOINT_INDEX); + + assertEquals(EXPECTED_MIDPOINT_BULGE_M, midpoint.curvatureAdjustedElevationMeters(), 0.05, + "curvature-adjusted elevation at 10 km"); + } + + /** + * Line of sight: 10 m antenna height minus 5.886 m bulge leaves 4.114 m, so the + * direct path is free and the minimum clearance of the whole profile is the midpoint + * value. Tolerance: 0.05 m. + */ + @Test + void lineOfSightIsFreeWithFourMetresClearanceAtTheMidpoint() { + PathAnalysisResult result = analyzeFlatPath(); + PathProfilePoint midpoint = result.profilePoints().get(MIDPOINT_INDEX); + + assertEquals(ANTENNA_HEIGHT_M, midpoint.lineOfSightHeightMeters(), 1e-9, "LOS height"); + assertEquals(EXPECTED_MIDPOINT_LOS_CLEARANCE_M, midpoint.lineOfSightClearanceMeters(), 0.05, + "LOS clearance at 10 km"); + assertTrue(result.lineOfSightClear(), "line of sight clear"); + assertEquals(EXPECTED_MIDPOINT_LOS_CLEARANCE_M, result.minimumLineOfSightClearanceMeters(), 0.05, + "minimum LOS clearance"); + } + + /** + * First Fresnel zone: R1 = sqrt(lambda d1 d2 / (d1 + d2)) = 101.921 m at the midpoint, + * far more than the 4.114 m of clearance, so the zone is intruded by + * 5.886 - (10 - 101.921) = 97.807 m at 10 km. Tolerance: 0.1 m. + */ + @Test + void fresnelZoneIsIntrudedByNinetyEightMetresAtTheMidpoint() { + PathAnalysisResult result = analyzeFlatPath(); + PathProfilePoint midpoint = result.profilePoints().get(MIDPOINT_INDEX); + + double localFresnelRadiusM = midpoint.lineOfSightHeightMeters() - midpoint.fresnelLowerHeightMeters(); + assertEquals(EXPECTED_MIDPOINT_FRESNEL_RADIUS_M, localFresnelRadiusM, 0.1, "R1 at 10 km"); + assertEquals(EXPECTED_MIDPOINT_INTRUSION_M, midpoint.fresnelIntrusionMeters(), 0.1, "intrusion at 10 km"); + + assertFalse(result.fresnelClear(), "Fresnel zone intruded"); + assertEquals(EXPECTED_MIDPOINT_INTRUSION_M, result.worstFresnelIntrusionMeters(), 0.1, "worst intrusion"); + assertEquals(10.0, result.worstFresnelDistanceKm(), 1e-9, "worst intrusion distance"); + assertEquals(MIDPOINT_INDEX, result.worstFresnelSampleIndex(), "worst intrusion sample"); + assertEquals(EXPECTED_MIDPOINT_INTRUSION_M / EXPECTED_MIDPOINT_FRESNEL_RADIUS_M, + result.worstFresnelIntrusionRatio(), 1e-3, "worst intrusion ratio"); + } + + /** + * No knife edge: the obstruction search only considers samples that rise above the + * line of sight (h greater than 0). On a flat profile with free LOS there is none, so + * there is no dominant obstruction and no knife-edge loss estimate. + */ + @Test + void noKnifeEdgeCandidateOnAFreeLineOfSight() { + PathObstructionSummary obstruction = analyzeFlatPath().obstructionSummary(); + + assertFalse(obstruction.hasDominantLosObstruction(), "dominant obstruction"); + assertEquals(-1, obstruction.dominantObstructionSampleIndex(), "obstruction sample index"); + assertTrue(Double.isNaN(obstruction.estimatedKnifeEdgeLossDb()), "knife-edge loss is NaN"); + assertEquals(EXPECTED_MIDPOINT_INTRUSION_M, obstruction.worstFresnelIntrusionMeters(), 0.1, + "obstruction summary carries the worst Fresnel intrusion"); + } + + /** + * Pinned, not endorsed: a partial Fresnel intrusion without LOS blockage adds 0 dB to + * the link budget. + * + *

The link budget takes its diffraction term only from a dominant LOS obstruction. + * Here the first Fresnel zone is intruded by 96 percent of its radius yet the direct + * path is open, so the code books 0 dB. ITU-R P.526 eq. (31) would still give a loss + * for negative nu down to -0.78 (here nu = -0.057, about 5.5 dB; in general 0 .. 6 dB + * for an open but Fresnel-obstructed path). Tracked as a separate issue; this test + * records the current behaviour. Tolerance: 0.01 dB.

+ */ + @Test + void partialFresnelIntrusionWithoutLosBlockageContributesZeroDiffractionLoss() { + PathLinkBudgetSummary budget = analyzeFlatPath().linkBudgetSummary(); + + assertEquals(0.0, budget.diffractionLossDb(), 0.01, "diffraction loss in the link budget"); + } +} diff --git a/core/src/test/java/kst4contest/view/map/PathGeometryUtilsReferenceTest.java b/core/src/test/java/kst4contest/view/map/PathGeometryUtilsReferenceTest.java new file mode 100644 index 00000000..03141249 --- /dev/null +++ b/core/src/test/java/kst4contest/view/map/PathGeometryUtilsReferenceTest.java @@ -0,0 +1,253 @@ +package kst4contest.view.map; + +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +/** + * Reference cases for {@link PathGeometryUtils} against ITU-R P.525 / P.526 and + * textbook values. + * + *

Every expected value is stated in the test method's Javadoc together with its + * source, the equation and the tolerance. Where the project deliberately deviates from + * the standard (the 32.44 constant of the free-space loss) the code's value is pinned + * and the standard's value is documented next to it. Nothing here was derived by + * running the code under test.

+ * + *

Context: issue #108.

+ */ +class PathGeometryUtilsReferenceTest { + + private static final double SPEED_OF_LIGHT_M_PER_S = 299_792_458.0; + + /** + * Great-circle distance and initial bearing, Veness worked example. + * + *

Source: Chris Veness, "Calculate distance, bearing and more between + * Latitude/Longitude points" (movable-type.co.uk/scripts/latlong.html): from + * 50 03 59 N, 005 42 53 W to 58 38 38 N, 003 04 12 W the haversine distance is + * 968.9 km on R = 6371 km and the initial bearing is 9.12 degrees. Recomputed on the + * project's R = 6371.009 km the distance is 968.855 km. Tolerance: 0.1 km, 0.01 deg.

+ */ + @Test + void greatCircleDistanceAndBearingMatchTheVenessWorkedExample() { + double fromLat = 50.0 + 3.0 / 60.0 + 59.0 / 3600.0; + double fromLon = -(5.0 + 42.0 / 60.0 + 53.0 / 3600.0); + double toLat = 58.0 + 38.0 / 60.0 + 38.0 / 3600.0; + double toLon = -(3.0 + 4.0 / 60.0 + 12.0 / 3600.0); + + assertEquals(968.855, + PathGeometryUtils.calculateGreatCircleDistanceKm(fromLat, fromLon, toLat, toLon), + 0.1, "Veness distance on R = 6371.009 km"); + assertEquals(9.1198, + PathGeometryUtils.calculateInitialBearingDeg(fromLat, fromLon, toLat, toLon), + 0.01, "Veness initial bearing"); + } + + /** + * First Fresnel zone radius, ITU-R P.526-15 section 2.1 eq. (2). + * + *

R_n = sqrt(n * lambda * d1 * d2 / (d1 + d2)) with n = 1 and lambda = c / f, + * c = 299 792 458 m/s. For d1 = d2 = 50 km at 144.3 MHz this is 227.9015 m, for + * d1 = 10 km, d2 = 40 km at 432.2 MHz it is 74.4926 m. Tolerance: 1e-3 m against the + * literal, relative 1e-6 against the equation evaluated in the test.

+ */ + @Test + void firstFresnelRadiusMatchesP526Equation2() { + double symmetric = PathGeometryUtils.calculateFirstFresnelRadiusMeters(50.0, 100.0, 144.3); + double asymmetric = PathGeometryUtils.calculateFirstFresnelRadiusMeters(10.0, 50.0, 432.2); + + assertEquals(227.9015, symmetric, 1e-3, "R1 for d1 = d2 = 50 km at 144.3 MHz"); + assertEquals(74.4926, asymmetric, 1e-3, "R1 for d1 = 10 km, d2 = 40 km at 432.2 MHz"); + + double expectedSymmetric = fresnelRadiusEq2(50_000.0, 50_000.0, 144.3); + double expectedAsymmetric = fresnelRadiusEq2(10_000.0, 40_000.0, 432.2); + assertEquals(expectedSymmetric, symmetric, expectedSymmetric * 1e-6, "R1 vs eq. (2), symmetric"); + assertEquals(expectedAsymmetric, asymmetric, expectedAsymmetric * 1e-6, "R1 vs eq. (2), asymmetric"); + } + + /** + * First Fresnel zone radius against the practical form, ITU-R P.526-15 eq. (3). + * + *

R_1 = 550 * sqrt(d1 * d2 / ((d1 + d2) * f)) with d in km and f in MHz. The + * constant 550 carries the rounding of sqrt(c) in these units, so the two forms + * differ by about 0.45 percent (228.928 m vs 227.902 m, 74.828 m vs 74.493 m). + * Tolerance: 0.5 percent.

+ */ + @Test + void firstFresnelRadiusAgreesWithP526Equation3WithinHalfAPercent() { + assertRelative(fresnelRadiusEq3(50.0, 50.0, 144.3), + PathGeometryUtils.calculateFirstFresnelRadiusMeters(50.0, 100.0, 144.3), + 0.005, "R1 vs eq. (3), symmetric"); + assertRelative(fresnelRadiusEq3(10.0, 40.0, 432.2), + PathGeometryUtils.calculateFirstFresnelRadiusMeters(10.0, 50.0, 432.2), + 0.005, "R1 vs eq. (3), asymmetric"); + } + + /** + * Recommended clearance is 60 percent of the first Fresnel radius. + * + *

Source: the usual planning rule (ITU-R P.530, "0.6 F1 clearance"). Expected + * 0.6 * 227.9015 m = 136.7409 m. Tolerance: 1e-3 m.

+ */ + @Test + void recommendedClearanceIsSixtyPercentOfTheFirstFresnelRadius() { + assertEquals(136.7409, + PathGeometryUtils.calculateRecommendedFresnelClearanceMeters(50.0, 100.0, 144.3), + 1e-3, "0.6 * R1 for d1 = d2 = 50 km at 144.3 MHz"); + } + + /** + * Earth bulge below the straight line between the endpoints. + * + *

h = d1 * d2 / (2 * k * R) with R = 6 371 009 m. For d1 = d2 = 50 km: 147.1509 m + * at k = 4/3 and 196.2013 m at k = 1. Tolerance: 1e-3 m. The k = 4/3 value is also + * checked against the effective Earth radius of 8500 km that ITU-R P.834 and the IEV + * quote for the standard atmosphere (147.0588 m), within 0.1 percent.

+ */ + @Test + void earthBulgeMatchesTheGeometricFormulaAndThe8500KmEffectiveRadius() { + double bulgeK43 = PathGeometryUtils.calculateEarthBulgeMeters(50.0, 100.0, 4.0 / 3.0); + double bulgeK1 = PathGeometryUtils.calculateEarthBulgeMeters(50.0, 100.0, 1.0); + + assertEquals(147.1509, bulgeK43, 1e-3, "bulge at k = 4/3"); + assertEquals(196.2013, bulgeK1, 1e-3, "bulge at k = 1"); + assertEquals(bulgeK43, PathGeometryUtils.calculateEarthBulgeMeters(50.0, 100.0), 1e-9, + "two-argument overload uses k = 4/3"); + + double bulge8500Km = 50_000.0 * 50_000.0 / (2.0 * 8_500_000.0); + assertRelative(bulge8500Km, bulgeK43, 0.001, "bulge vs 8500 km effective radius"); + } + + /** + * Radio horizon distance for an antenna height of 100 m. + * + *

d = sqrt(2 * k * R * h) with R = 6 371 009 m: 41.2181 km at k = 4/3 and + * 35.6960 km at k = 1. Tolerance: 1e-3 km. The textbook approximations + * d = 4.12 * sqrt(h) (radio, k = 4/3) and d = 3.57 * sqrt(h) (optical, k = 1), d in + * km and h in m, are met within 0.1 percent.

+ */ + @Test + void radioHorizonMatchesSqrtTwoKRHAndTheTextbookApproximations() { + double horizonK43 = PathGeometryUtils.calculateRadioHorizonDistanceKm(100.0, 4.0 / 3.0); + double horizonK1 = PathGeometryUtils.calculateRadioHorizonDistanceKm(100.0, 1.0); + + assertEquals(41.2181, horizonK43, 1e-3, "horizon at k = 4/3"); + assertEquals(35.6960, horizonK1, 1e-3, "horizon at k = 1"); + assertRelative(4.12 * Math.sqrt(100.0), horizonK43, 0.001, "horizon vs 4.12 sqrt(h)"); + assertRelative(3.57 * Math.sqrt(100.0), horizonK1, 0.001, "horizon vs 3.57 sqrt(h)"); + } + + /** + * Single knife-edge diffraction loss, ITU-R P.526-15 eq. (31). + * + *

J(nu) = 6.9 + 20 log10(sqrt((nu - 0.1)^2 + 1) + nu - 0.1) for nu greater than + * -0.78, otherwise 0. Expected: J(0) = 6.033, J(1) = 13.926, J(2.4) = 20.539, + * J(-0.5) = 1.959, J(-1) = 0. Tolerance: 0.01 dB.

+ * + *

At exactly nu = -0.78 the formula gives 0.004 dB; the code compares with + * {@code v <= -0.78} and answers 0 there. That is pinned here, not endorsed: the step + * is 0.004 dB and of no practical consequence, but the boundary belongs to the + * zero branch in this implementation.

+ */ + @Test + void knifeEdgeLossMatchesP526Equation31() { + assertEquals(6.033, PathGeometryUtils.calculateSingleKnifeEdgeLossDb(0.0), 0.01, "J(0)"); + assertEquals(13.926, PathGeometryUtils.calculateSingleKnifeEdgeLossDb(1.0), 0.01, "J(1)"); + assertEquals(20.539, PathGeometryUtils.calculateSingleKnifeEdgeLossDb(2.4), 0.01, "J(2.4)"); + assertEquals(1.959, PathGeometryUtils.calculateSingleKnifeEdgeLossDb(-0.5), 0.01, "J(-0.5)"); + assertEquals(0.0, PathGeometryUtils.calculateSingleKnifeEdgeLossDb(-1.0), 0.0, "J(-1)"); + // Pinned: the boundary nu = -0.78 falls into the zero branch (v <= -0.78). + assertEquals(0.0, PathGeometryUtils.calculateSingleKnifeEdgeLossDb(-0.78), 0.0, "J(-0.78)"); + } + + /** + * Fresnel-Kirchhoff diffraction parameter, ITU-R P.526-15 eq. (26) in SI units. + * + *

nu = h * sqrt(2 * (d1 + d2) / (lambda * d1 * d2)). For h = 50 m, d1 = d2 = 25 km + * at 144.3 MHz (lambda = 2.0776 m): nu = 0.4388. Tolerance: 1e-4.

+ */ + @Test + void knifeEdgeVParameterMatchesP526Equation26() { + assertEquals(0.4388, + PathGeometryUtils.calculateKnifeEdgeVParameter(50.0, 25.0, 25.0, 144.3), + 1e-4, "nu for h = 50 m, d1 = d2 = 25 km, 144.3 MHz"); + } + + /** + * Free-space basic transmission loss for 100 km at 144.3 MHz. + * + *

The code uses L = 32.44 + 20 log10(d/km) + 20 log10(f/MHz), which gives + * 115.625 dB; that value is pinned (tolerance 0.001 dB) as decided in #108.

+ * + *

For the record: ITU-R P.525-5 eq. (6) writes the constant as 32.4 (115.585 dB + * here), and the exact constant 20 log10(4 pi * 1e9 / c) = 32.4478 gives 115.633 dB. + * The project keeps 32.44; the spread of 0.05 dB is far below the uncertainty of + * the rest of the link budget.

+ */ + @Test + void freeSpacePathLossUsesTheProjectConstant32Point44() { + double lossDb = PathGeometryUtils.calculateFreeSpacePathLossDb(100.0, 144.3); + + assertEquals(115.625, lossDb, 0.001, "FSPL 100 km / 144.3 MHz with 32.44"); + + // Documented, not asserted as the project's answer: the standard's variants. + double p525Variant = 32.4 + 20.0 * Math.log10(100.0) + 20.0 * Math.log10(144.3); + double exactConstant = 20.0 * Math.log10(4.0 * Math.PI * 1e9 / SPEED_OF_LIGHT_M_PER_S); + double exactVariant = exactConstant + 20.0 * Math.log10(100.0) + 20.0 * Math.log10(144.3); + assertEquals(115.585, p525Variant, 0.001, "P.525-5 eq. (6) with 32.4"); + assertEquals(115.633, exactVariant, 0.001, "exact constant 32.4478"); + } + + /** + * Watts to dBm: P_dBm = 10 log10(P_W * 1000). + * + *

750 W = 58.751 dBm, 100 W = 50.000 dBm. Tolerance: 1e-3 dB.

+ */ + @Test + void wattsToDbmMatchesTheDefinition() { + assertEquals(58.751, PathGeometryUtils.wattsToDbm(750.0), 1e-3, "750 W"); + assertEquals(50.0, PathGeometryUtils.wattsToDbm(100.0), 1e-3, "100 W"); + } + + /** + * Adaptive sample count: ceil(d / 0.5 km) + 1, clamped to 121 .. 1201. + * + *

0 km (non-positive) yields the minimum 121, 50 km computes 101 and is clamped + * to 121, 100 km gives 201, 1000 km computes 2001 and is clamped to 1201.

+ */ + @Test + void adaptiveSampleCountFollowsTheHalfKilometreRuleWithClamping() { + assertEquals(121, PathGeometryUtils.resolveAdaptiveSampleCount(0.0), "0 km"); + assertEquals(121, PathGeometryUtils.resolveAdaptiveSampleCount(50.0), "50 km"); + assertEquals(201, PathGeometryUtils.resolveAdaptiveSampleCount(100.0), "100 km"); + assertEquals(1201, PathGeometryUtils.resolveAdaptiveSampleCount(1000.0), "1000 km"); + } + + /** + * Effective Earth radius factor sanitising: values outside 0.5 .. 10 fall back to 4/3. + * + *

0.4 and 11 are replaced by 4/3, 1.0 is kept. Tolerance: 1e-12.

+ */ + @Test + void effectiveEarthRadiusFactorOutsideTheAcceptedRangeFallsBackToFourThirds() { + assertEquals(4.0 / 3.0, PathGeometryUtils.sanitizeEffectiveEarthRadiusFactor(0.4), 1e-12, "0.4"); + assertEquals(4.0 / 3.0, PathGeometryUtils.sanitizeEffectiveEarthRadiusFactor(11.0), 1e-12, "11"); + assertEquals(1.0, PathGeometryUtils.sanitizeEffectiveEarthRadiusFactor(1.0), 1e-12, "1.0"); + } + + /** ITU-R P.526-15 eq. (2) with n = 1, distances in metres. */ + private static double fresnelRadiusEq2(double d1Meters, double d2Meters, double frequencyMHz) { + double wavelengthMeters = SPEED_OF_LIGHT_M_PER_S / (frequencyMHz * 1e6); + return Math.sqrt(wavelengthMeters * d1Meters * d2Meters / (d1Meters + d2Meters)); + } + + /** ITU-R P.526-15 eq. (3): distances in km, frequency in MHz, result in metres. */ + private static double fresnelRadiusEq3(double d1Km, double d2Km, double frequencyMHz) { + return 550.0 * Math.sqrt(d1Km * d2Km / ((d1Km + d2Km) * frequencyMHz)); + } + + private static void assertRelative(double expected, double actual, double relativeTolerance, String message) { + assertEquals(expected, actual, Math.abs(expected) * relativeTolerance, message); + } +}