mobility: Error in computation of elevation angle in class GeocentricConstantPositionMobilityModel

Inside file geocentric-constant-position-mobility-model.cc, in method GetElevationAngle, the following code is present:

// DoGetGeocentricPosition:
// https://www.nsnam.org/doxygen/d3/d1f/classns3_1_1_geocentric_constant_position_mobility_model.html#abe513712e9b30269a7f0db501e4641ea

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);

However, since the Vector objects me and them contain the position in geocentric coordinates, the component z is not related to the altitude but represents instead the distance from the Equatorial plane.

With the current code, any geostationary satellite (which by definition orbits on the Equatorial plane) has component z equal to 0 and is take by the code as the terminal with the lowest altitude if the other point is located in the northern hemisphere. This causes the elevation angle to be computed incorrectly.

The correct code would use the length of the vector as proxy for the altitude

Vector& a = (me.GetLength() < them.GetLength() ? me : them);
Vector& b = (me.GetLength() < them.GetLength() ? them : me);

Another error is the application of std::abs, both in the computation of the numerator

double numerator = std::abs(a * bMinusA);

and of the final elevation angle

// x:
// https://www.nsnam.org/doxygen/d2/d34/cairo-wideint_8c.html#a18faf656404f0e3e7e86226b22cff4bb

double elevAngle = std::abs((180.0 * M_1_PI) * asin(x));

Taking the absolute value forces the elevation angle to be positive, returning a positive value even for a situation where the satellite is below the horizon of the other point. In this case, instead, the returned elevation angle must be negative.

Edited by Eduardo Almeida