mobility: (fixes #1308) Fix elevation angle computation

GeocentricConstantPositionMobilityModel::DoGetElevationAngle selected the "lower" terminal by ECEF z-coordinate, which is the distance from the equatorial plane rather than the altitude, so for many geometries the wrong terminal was chosen as the horizon vertex; it also took std::abs of the dot-product numerator and of the final angle, so a below-horizon terminal reported a positive elevation.

Select the lower terminal by distance from the Earth's centre and keep the signs, so the elevation is correct, negative below the horizon, and independent of argument order. Add a unit test covering zenith satellites across a latitude sweep, a GEO satellite seen from 45 N and its southern-hemisphere mirror, oblique north/south mirror geometries at high latitude (which the z-based vertex selection broke), longitude rotational invariance, an exact-horizon (0 degree) terminal, argument-order independence, and a below-horizon terminal.

Author: Gabrielcarvfer

src/mobility/CMakeLists.txt

@@ -61,6 +61,7 @@ build_lib(
     test/box-line-intersection-test.cc
     test/geo-to-cartesian-test.cc
     test/geocentric-topocentric-conversion-test.cc
+    test/geocentric-elevation-angle-test.cc
     test/mobility-test-suite.cc
     test/mobility-trace-test-suite.cc
     test/ns2-mobility-helper-test-suite.cc

src/mobility/model/geocentric-constant-position-mobility-model.cc

@@ -195,21 +195,29 @@ GeocentricConstantPositionMobilityModel::DoGetElevationAngle(
     Vector me = this->DoGetGeocentricPosition();
     Vector them = other->DoGetGeocentricPosition();
 
-    // a is assumed to be the terminal with the lowest altitude
-    Vector& a = (me.z < them.z ? me : them);
-    Vector& b = (me.z < them.z ? them : me);
+    // a is the terminal with the lowest altitude, used as the vertex whose local
+    // vertical (the radial direction a/|a|) defines the horizon. The altitude is
+    // approximated by the distance from the Earth's centre, not by the ECEF
+    // z-coordinate, which is the distance from the equatorial plane
+    // (@issueid{1308}).
+    Vector& a = (me.GetLength() < them.GetLength() ? me : them);
+    Vector& b = (me.GetLength() < them.GetLength() ? them : me);
 
     Vector bMinusA = b - a;
-    double numerator = std::abs(a * bMinusA);
+    // sin(elevation) = (a . (b - a)) / (|a| |b - a|). The dot product is signed:
+    // it is negative when the other terminal is below the local horizon, so the
+    // sign must be preserved.
+    double numerator = a * bMinusA;
     double denominator = a.GetLength() * bMinusA.GetLength();
     double x = numerator / denominator;
 
     // Enforce the appropriate domain range ([-1, 1]) for the argument of std::asin.
     x = std::min(x, 1.0);
     x = std::max(x, -1.0);
 
-    // asin returns radians, we convert to degrees
-    double elevAngle = std::abs((180.0 * M_1_PI) * asin(x));
+    // asin returns radians, we convert to degrees. A negative result means the
+    // other terminal is below the local horizon.
+    double elevAngle = (180.0 * M_1_PI) * asin(x);
 
     return elevAngle;
 }

src/mobility/model/geocentric-constant-position-mobility-model.h

@@ -50,8 +50,9 @@ class GeocentricConstantPositionMobilityModel : public MobilityModel
      * After calculating the plane perpendicular to a cartesian position vector,
      * the elevation angle is calculated using
      * https://www.w3schools.blog/angle-between-a-line-and-a-plane.
-     * The altitude of the position passed as a parameter must be higher than that of the reference
-     * point.
+     * The elevation is measured at the lower of the two terminals (the one closer
+     * to the Earth's centre); the result is independent of the argument order. A
+     * negative value indicates that the other terminal is below the local horizon.
      * @param other pointer to the HAPS/satellite mobility model
      * @return the elevation angle
      * in degrees

src/mobility/test/geocentric-elevation-angle-test.cc

@@ -0,0 +1,156 @@
+/*
+ * Copyright (c) 2026 Centre Tecnologic de Telecomunicacions de Catalunya (CTTC)
+ *
+ * SPDX-License-Identifier: GPL-2.0-only
+ *
+ * Author: Gabriel Ferreira <gabrielcarvfer@gmail.com>
+ */
+
+#include "ns3/geocentric-constant-position-mobility-model.h"
+#include "ns3/geographic-positions.h"
+#include "ns3/test.h"
+
+#include <cmath>
+
+using namespace ns3;
+
+/**
+ * @ingroup mobility-test
+ *
+ * @brief Tests GeocentricConstantPositionMobilityModel::GetElevationAngle (@issueid{1308}).
+ *
+ * The elevation angle is measured at the lower terminal, with the local
+ * vertical given by the geocentric radial direction (the model uses a perfect
+ * sphere). The angle must be signed (negative when the other terminal is
+ * below the local horizon), independent of the argument order, and identical
+ * for mirrored geometries in either hemisphere. This test pins that behavior.
+ */
+class GeocentricElevationAngleTestCase : public TestCase
+{
+  public:
+    GeocentricElevationAngleTestCase()
+        : TestCase("GeocentricConstantPositionMobilityModel elevation angle")
+    {
+    }
+
+  private:
+    void DoRun() override;
+
+    /// Build a model at the given geographic position.
+    /// @param lat latitude (deg)
+    /// @param lon longitude (deg)
+    /// @param alt altitude (m)
+    /// @return the configured mobility model
+    static Ptr<GeocentricConstantPositionMobilityModel> Make(double lat, double lon, double alt)
+    {
+        auto m = CreateObject<GeocentricConstantPositionMobilityModel>();
+        m->SetGeographicPosition(Vector(lat, lon, alt));
+        return m;
+    }
+};
+
+void
+GeocentricElevationAngleTestCase::DoRun()
+{
+    constexpr double TOL = 1e-2; // degrees
+    const double geoAlt = 35786e3;
+    const double leoAlt = 600e3;
+
+    // A satellite directly above the ground station is at zenith: 90 degrees.
+    auto gs = Make(45.0, 0.0, 0.0);
+    auto overhead = Make(45.0, 0.0, geoAlt);
+    NS_TEST_ASSERT_MSG_EQ_TOL(gs->GetElevationAngle(overhead),
+                              90.0,
+                              TOL,
+                              "Satellite at zenith should have 90 degrees elevation");
+
+    // GEO satellite over the equator seen from a ground station at 45 N.
+    auto geoSat = Make(0.0, 0.0, geoAlt);
+    double elevFromGs = gs->GetElevationAngle(geoSat);
+    NS_TEST_ASSERT_MSG_EQ_TOL(elevFromGs, 38.1771, TOL, "Wrong elevation for GEO over the equator");
+    // The result must be independent of the argument order.
+    NS_TEST_ASSERT_MSG_EQ_TOL(geoSat->GetElevationAngle(gs),
+                              elevFromGs,
+                              TOL,
+                              "Elevation angle must not depend on the argument order");
+
+    // A terminal below the local horizon must yield a negative elevation.
+    auto equatorGs = Make(0.0, 0.0, 0.0);
+    auto farSat = Make(0.0, 120.0, leoAlt);
+    NS_TEST_ASSERT_MSG_EQ_TOL(equatorGs->GetElevationAngle(farSat),
+                              -58.5127,
+                              TOL,
+                              "Below-horizon terminal should have a negative elevation");
+
+    // Southern-hemisphere mirror of the GEO case: the geometry is symmetric,
+    // so the elevation must match the northern result.
+    auto southGs = Make(-45.0, 0.0, 0.0);
+    NS_TEST_ASSERT_MSG_EQ_TOL(southGs->GetElevationAngle(geoSat),
+                              elevFromGs,
+                              TOL,
+                              "Southern-hemisphere GEO elevation must mirror the northern one");
+
+    // Oblique LEO geometry at high latitude and its southern mirror: the
+    // endpoints have genuinely different local verticals, so hemisphere
+    // symmetry is a strong discriminator for the vertex selection.
+    auto northHighGs = Make(60.0, 0.0, 0.0);
+    auto northLeo = Make(50.0, 0.0, leoAlt);
+    double northElev = northHighGs->GetElevationAngle(northLeo);
+    NS_TEST_ASSERT_MSG_EQ_TOL(northElev, 22.2040, TOL, "Wrong elevation for oblique LEO geometry");
+    auto southHighGs = Make(-60.0, 0.0, 0.0);
+    auto southLeo = Make(-50.0, 0.0, leoAlt);
+    NS_TEST_ASSERT_MSG_EQ_TOL(southHighGs->GetElevationAngle(southLeo),
+                              northElev,
+                              TOL,
+                              "North and south mirror geometries must have equal elevation");
+
+    // Longitude rotational invariance: shifting both terminals by the same
+    // longitude must not change the (purely relative) geometry.
+    auto rotatedGs = Make(60.0, 130.0, 0.0);
+    auto rotatedLeo = Make(50.0, 130.0, leoAlt);
+    NS_TEST_ASSERT_MSG_EQ_TOL(rotatedGs->GetElevationAngle(rotatedLeo),
+                              northElev,
+                              TOL,
+                              "Elevation must be invariant under a common longitude shift");
+
+    // Exact horizon: place the satellite so the line of sight grazes the
+    // sphere and the ground station sees it at exactly 0 degrees, pinning the
+    // sign boundary between visible and below-horizon terminals.
+    const double horizonLon = std::acos(GeographicPositions::EARTH_SPHERE_RADIUS /
+                                        (GeographicPositions::EARTH_SPHERE_RADIUS + leoAlt)) *
+                              180.0 / M_PI;
+    auto horizonSat = Make(0.0, horizonLon, leoAlt);
+    NS_TEST_ASSERT_MSG_EQ_TOL(equatorGs->GetElevationAngle(horizonSat),
+                              0.0,
+                              TOL,
+                              "Satellite at the geometric horizon should have 0 degrees elevation");
+
+    // Zenith sanity sweep across latitudes. Zenith is degenerate (both radial
+    // directions coincide), so this is an invariant check.
+    for (double lat : {-80.0, -30.0, 0.0, 30.0, 80.0})
+    {
+        NS_TEST_ASSERT_MSG_EQ_TOL(Make(lat, 0.0, 0.0)->GetElevationAngle(Make(lat, 0.0, leoAlt)),
+                                  90.0,
+                                  TOL,
+                                  "Zenith satellite should have 90 degrees elevation at latitude "
+                                      << lat);
+    }
+}
+
+/**
+ * @ingroup mobility-test
+ *
+ * @brief Geocentric elevation angle test suite.
+ */
+class GeocentricElevationAngleTestSuite : public TestSuite
+{
+  public:
+    GeocentricElevationAngleTestSuite()
+        : TestSuite("geocentric-elevation-angle", Type::UNIT)
+    {
+        AddTestCase(new GeocentricElevationAngleTestCase, TestCase::Duration::QUICK);
+    }
+};
+
+static GeocentricElevationAngleTestSuite
+    g_geocentricElevationAngleTestSuite; //!< Static variable for test initialization