diff --git a/Migration.md b/Migration.md index 5305eaa59..5c81932ab 100644 --- a/Migration.md +++ b/Migration.md @@ -24,6 +24,10 @@ release will remove the deprecated code. + All mutator functions that lacked a `Set` prefix have been deprecated and replaced by version with a `Set` prefix. +1. **SphericalCoordinates.hh** + + ***Deprecation:*** public: static double Distance(const Angle&, const Angle&, const Angle&, const Angle&) + + ***Replacement:*** public: static double DistanceWGS84(const Angle&, const Angle&, const Angle&, const Angle&) + 1. **Matrix3.hh** + ***Deprecation:*** public: void Axes(const Vector3 &, const Vector3 &, const Vector3 &) + ***Replacement:*** public: void SetAxes(const Vector3 &, const Vector3 &, const Vector3 &) diff --git a/include/gz/math/SphericalCoordinates.hh b/include/gz/math/SphericalCoordinates.hh index 7ea183a6b..546939646 100644 --- a/include/gz/math/SphericalCoordinates.hh +++ b/include/gz/math/SphericalCoordinates.hh @@ -41,7 +41,15 @@ namespace gz { /// \brief Model of reference ellipsoid for earth, based on /// WGS 84 standard. see wikipedia: World_Geodetic_System - EARTH_WGS84 = 1 + EARTH_WGS84 = 1, + + /// \brief Model of the moon, based on the Selenographic + /// coordinate system, see wikipedia: Selenographic + /// Coordinate System. + MOON_SCS = 2, + + /// \brief Custom surface type + CUSTOM_SURFACE = 10 }; /// \enum CoordinateType @@ -73,6 +81,16 @@ namespace gz /// \param[in] _type SurfaceType specification. public: explicit SphericalCoordinates(const SurfaceType _type); + /// \brief Constructor with surface type input and properties + /// input. To be used for CUSTOM_SURFACE. + /// \param[in] _type SurfaceType specification. + /// \param[in] _axisEquatorial Semi major axis of the surface. + /// \param[in] _axisPolar Semi minor axis of the surface. + public: SphericalCoordinates( + const SurfaceType _type, + const double _axisEquatorial, + const double _axisPolar); + /// \brief Constructor with surface type, angle, and elevation inputs. /// \param[in] _type SurfaceType specification. /// \param[in] _latitude Reference latitude. @@ -85,7 +103,6 @@ namespace gz const double _elevation, const gz::math::Angle &_heading); - /// \brief Convert a Cartesian position vector to geodetic coordinates. /// This performs a `PositionTransform` from LOCAL to SPHERICAL. /// @@ -136,15 +153,66 @@ namespace gz /// \param[in] _latB Latitude of point B. /// \param[in] _lonB Longitude of point B. /// \return Distance in meters. - public: static double Distance(const gz::math::Angle &_latA, + /// \deprecated Use DistanceWGS84 instead. + public: GZ_DEPRECATED(7) static double Distance( + const gz::math::Angle &_latA, + const gz::math::Angle &_lonA, + const gz::math::Angle &_latB, + const gz::math::Angle &_lonB); + + /// \brief Get the distance between two points expressed in geographic + /// latitude and longitude. It assumes that both points are at sea level. + /// Example: _latA = 38.0016667 and _lonA = -123.0016667) represents + /// the point with latitude 38d 0'6.00"N and longitude 123d 0'6.00"W. + /// This method assumes that the surface model is EARTH_WGS84. + /// \param[in] _latA Latitude of point A. + /// \param[in] _lonA Longitude of point A. + /// \param[in] _latB Latitude of point B. + /// \param[in] _lonB Longitude of point B. + /// \return Distance in meters. + public: static double DistanceWGS84( + const gz::math::Angle &_latA, const gz::math::Angle &_lonA, const gz::math::Angle &_latB, const gz::math::Angle &_lonB); + /// \brief Get the distance between two points expressed in geographic + /// latitude and longitude. It assumes that both points are at sea level. + /// Example: _latA = 38.0016667 and _lonA = -123.0016667) represents + /// the point with latitude 38d 0'6.00"N and longitude 123d 0'6.00"W. + /// This is different from the deprecated static Distance() method as it + /// takes into account the set surface's radius. + /// \param[in] _latA Latitude of point A. + /// \param[in] _lonA Longitude of point A. + /// \param[in] _latB Latitude of point B. + /// \param[in] _lonB Longitude of point B. + /// \return Distance in meters. + public: double DistanceBetweenPoints( + const gz::math::Angle &_latA, + const gz::math::Angle &_lonA, + const gz::math::Angle &_latB, + const gz::math::Angle &_lonB); + /// \brief Get SurfaceType currently in use. /// \return Current SurfaceType value. public: SurfaceType Surface() const; + /// \brief Get the radius of the surface. + /// \return radius of the surface in use. + public: double SurfaceRadius(); + + /// \brief Get the major axis of the surface. + /// \return Equatorial axis of the surface in use. + public: double SurfaceAxisEquatorial(); + + /// \brief Get the minor axis of the surface. + /// \return Polar axis of the surface in use. + public: double SurfaceAxisPolar(); + + /// \brief Get the flattening of the surface. + /// \return Flattening parameter of the surface in use. + public: double SurfaceFlattening(); + /// \brief Get reference geodetic latitude. /// \return Reference geodetic latitude. public: gz::math::Angle LatitudeReference() const; @@ -167,6 +235,16 @@ namespace gz /// \param[in] _type SurfaceType value. public: void SetSurface(const SurfaceType &_type); + /// \brief Set SurfaceType for planetary surface model with + /// custom ellipsoid properties. + /// \param[in] _type SurfaceType value. + /// \param[in] _axisEquatorial Equatorial axis of the surface. + /// \param[in] _axisPolar Polar axis of the surface. + public: void SetSurface( + const SurfaceType &_type, + const double _axisEquatorial, + const double _axisPolar); + /// \brief Set reference geodetic latitude. /// \param[in] _angle Reference geodetic latitude. public: void SetLatitudeReference(const gz::math::Angle &_angle); diff --git a/src/SphericalCoordinates.cc b/src/SphericalCoordinates.cc index 9c36a1d3c..87934ee36 100644 --- a/src/SphericalCoordinates.cc +++ b/src/SphericalCoordinates.cc @@ -38,12 +38,31 @@ const double g_EarthWGS84Flattening = 1.0/298.257223563; // Radius of the Earth (meters). const double g_EarthRadius = 6371000.0; +// Radius of the Moon (meters). +// Source: https://lunar.gsfc.nasa.gov/library/451-SCI-000958.pdf +const double g_MoonRadius = 1737400.0; + +// a: Equatorial radius of the Moon. +// Source : https://nssdc.gsfc.nasa.gov/planetary/factsheet/moonfact.html +const double g_MoonAxisEquatorial = 1738100.0; + +// b: Polar radius of the Moon. +// Source : https://nssdc.gsfc.nasa.gov/planetary/factsheet/moonfact.html +const double g_MoonAxisPolar = 1736000.0; + +// if: Unitless flattening parameter for the Moon. +// Source : https://nssdc.gsfc.nasa.gov/planetary/factsheet/moonfact.html +const double g_MoonFlattening = 0.0012; + // Private data for the SphericalCoordinates class. class gz::math::SphericalCoordinates::Implementation { /// \brief Type of surface being used. public: SphericalCoordinates::SurfaceType surfaceType; + /// \brief Radius of the given SurfaceType. + public: double surfaceRadius = 0; + /// \brief Latitude of reference point. public: gz::math::Angle latitudeReference; @@ -94,6 +113,10 @@ SphericalCoordinates::SurfaceType SphericalCoordinates::Convert( { if ("EARTH_WGS84" == _str) return EARTH_WGS84; + else if ("MOON_SCS" == _str) + return MOON_SCS; + else if ("CUSTOM_SURFACE" == _str) + return CUSTOM_SURFACE; std::cerr << "SurfaceType string not recognized, " << "EARTH_WGS84 returned by default" << std::endl; @@ -106,6 +129,10 @@ std::string SphericalCoordinates::Convert( { if (_type == EARTH_WGS84) return "EARTH_WGS84"; + else if (_type == MOON_SCS) + return "MOON_SCS"; + else if (_type == CUSTOM_SURFACE) + return "CUSTOM_SURFACE"; std::cerr << "SurfaceType not recognized, " << "EARTH_WGS84 returned by default" << std::endl; @@ -128,6 +155,20 @@ SphericalCoordinates::SphericalCoordinates(const SurfaceType _type) this->SetElevationReference(0.0); } +////////////////////////////////////////////////// +SphericalCoordinates::SphericalCoordinates( + const SurfaceType _type, + const double _axisEquatorial, + const double _axisPolar) + : SphericalCoordinates() +{ + // Set properties + this->SetSurface(_type, _axisEquatorial, + _axisPolar); + + this->SetElevationReference(0.0); +} + ////////////////////////////////////////////////// SphericalCoordinates::SphericalCoordinates(const SurfaceType _type, const gz::math::Angle &_latitude, @@ -208,6 +249,42 @@ void SphericalCoordinates::SetSurface(const SurfaceType &_type) std::pow(this->dataPtr->ellA, 2) / std::pow(this->dataPtr->ellB, 2) - 1.0); + // Set the radius of the surface. + this->dataPtr->surfaceRadius = g_EarthRadius; + + break; + } + case MOON_SCS: + { + // Set the semi-major axis + this->dataPtr->ellA = g_MoonAxisEquatorial; + + // Set the semi-minor axis + this->dataPtr->ellB = g_MoonAxisPolar; + + // Set the flattening parameter + this->dataPtr->ellF = g_MoonFlattening; + + // Set the first eccentricity ellipse parameter + // https://en.wikipedia.org/wiki/Eccentricity_(mathematics)#Ellipses + this->dataPtr->ellE = sqrt(1.0 - + std::pow(this->dataPtr->ellB, 2) / std::pow(this->dataPtr->ellA, 2)); + + // Set the second eccentricity ellipse parameter + // https://en.wikipedia.org/wiki/Eccentricity_(mathematics)#Ellipses + this->dataPtr->ellP = sqrt( + std::pow(this->dataPtr->ellA, 2) / std::pow(this->dataPtr->ellB, 2) - + 1.0); + + // Set the radius of the surface. + this->dataPtr->surfaceRadius = g_MoonRadius; + + break; + } + case CUSTOM_SURFACE: + { + std::cerr << "For custom surfaces, use SetSurface(type, radius," + "axisEquatorial, axisPolar)" << std::endl; break; } default: @@ -219,6 +296,53 @@ void SphericalCoordinates::SetSurface(const SurfaceType &_type) } } +////////////////////////////////////////////////// +void SphericalCoordinates::SetSurface( + const SurfaceType &_type, + const double _axisEquatorial, + const double _axisPolar) +{ + if ((_type != EARTH_WGS84) && + (_type != MOON_SCS) && + (_type != CUSTOM_SURFACE)) + { + std::cerr << "Unknown surface type[" + << _type << "]\n"; + return; + } + + this->dataPtr->surfaceType = _type; + + if ((_axisEquatorial > 0) + && (_axisPolar > 0) + && (_axisPolar <= _axisEquatorial)) + { + this->dataPtr->ellA = _axisEquatorial; + this->dataPtr->ellB = _axisPolar; + this->dataPtr->ellF = + (_axisEquatorial - _axisPolar) / _axisEquatorial; + // Arithmetic mean radius + this->dataPtr->surfaceRadius = + (2 * _axisEquatorial + _axisPolar) / 3.0; + } + else + { + std::cerr << "Invalid parameters found, defaulting to " + "Earth's parameters" << std::endl; + + this->dataPtr->ellA = g_EarthWGS84AxisEquatorial; + this->dataPtr->ellB = g_EarthWGS84AxisPolar; + this->dataPtr->ellF = g_EarthWGS84Flattening; + this->dataPtr->surfaceRadius = g_EarthRadius; + } + + this->dataPtr->ellE = sqrt(1.0 - + std::pow(this->dataPtr->ellB, 2) / std::pow(this->dataPtr->ellA, 2)); + this->dataPtr->ellP = sqrt( + std::pow(this->dataPtr->ellA, 2) / std::pow(this->dataPtr->ellB, 2) - + 1.0); +} + ////////////////////////////////////////////////// void SphericalCoordinates::SetLatitudeReference( const gz::math::Angle &_angle) @@ -286,7 +410,7 @@ gz::math::Vector3d SphericalCoordinates::LocalFromGlobalVelocity( ////////////////////////////////////////////////// /// Based on Haversine formula (http://en.wikipedia.org/wiki/Haversine_formula). -double SphericalCoordinates::Distance(const gz::math::Angle &_latA, +double SphericalCoordinates::DistanceWGS84(const gz::math::Angle &_latA, const gz::math::Angle &_lonA, const gz::math::Angle &_latB, const gz::math::Angle &_lonB) @@ -303,6 +427,62 @@ double SphericalCoordinates::Distance(const gz::math::Angle &_latA, return d; } +////////////////////////////////////////////////// +/// Based on Haversine formula (http://en.wikipedia.org/wiki/Haversine_formula). +double SphericalCoordinates::Distance(const gz::math::Angle &_latA, + const gz::math::Angle &_lonA, + const gz::math::Angle &_latB, + const gz::math::Angle &_lonB) +{ + return gz::math::SphericalCoordinates::DistanceWGS84( + _latA, _lonA, _latB, _lonB); +} + +////////////////////////////////////////////////// +/// Based on Haversine formula (http://en.wikipedia.org/wiki/Haversine_formula). +/// This takes into account the surface type. +double SphericalCoordinates::DistanceBetweenPoints( + const gz::math::Angle &_latA, + const gz::math::Angle &_lonA, + const gz::math::Angle &_latB, + const gz::math::Angle &_lonB) +{ + gz::math::Angle dLat = _latB - _latA; + gz::math::Angle dLon = _lonB - _lonA; + + double a = sin(dLat.Radian() / 2) * sin(dLat.Radian() / 2) + + sin(dLon.Radian() / 2) * sin(dLon.Radian() / 2) * + cos(_latA.Radian()) * cos(_latB.Radian()); + + double c = 2 * atan2(sqrt(a), sqrt(1 - a)); + double d = this->dataPtr->surfaceRadius * c; + return d; +} + +////////////////////////////////////////////////// +double SphericalCoordinates::SurfaceRadius() +{ + return this->dataPtr->surfaceRadius; +} + +////////////////////////////////////////////////// +double SphericalCoordinates::SurfaceAxisEquatorial() +{ + return this->dataPtr->ellA; +} + +////////////////////////////////////////////////// +double SphericalCoordinates::SurfaceAxisPolar() +{ + return this->dataPtr->ellB; +} + +////////////////////////////////////////////////// +double SphericalCoordinates::SurfaceFlattening() +{ + return this->dataPtr->ellF; +} + ////////////////////////////////////////////////// void SphericalCoordinates::UpdateTransformationMatrix() { diff --git a/src/SphericalCoordinates_TEST.cc b/src/SphericalCoordinates_TEST.cc index 3468da276..26eee5552 100644 --- a/src/SphericalCoordinates_TEST.cc +++ b/src/SphericalCoordinates_TEST.cc @@ -63,6 +63,11 @@ TEST(SphericalCoordinatesTest, Constructor) math::SphericalCoordinates sc2(sc); EXPECT_EQ(sc, sc2); } + + // Bad surface type, this should throw an error + math::SphericalCoordinates invalidSC( + static_cast(3)); + EXPECT_EQ(invalidSC.Surface(), 3); } ////////////////////////////////////////////////// @@ -80,7 +85,17 @@ TEST(SphericalCoordinatesTest, Convert) EXPECT_EQ("EARTH_WGS84", math::SphericalCoordinates::Convert(st)); EXPECT_EQ("EARTH_WGS84", math::SphericalCoordinates::Convert( - static_cast(2))); + static_cast(3))); + + // For the Moon surface type + st = math::SphericalCoordinates::MOON_SCS; + EXPECT_EQ(math::SphericalCoordinates::Convert("MOON_SCS"), st); + EXPECT_EQ("MOON_SCS", math::SphericalCoordinates::Convert(st)); + + // For the custom surface type + st = math::SphericalCoordinates::CUSTOM_SURFACE; + EXPECT_EQ(math::SphericalCoordinates::Convert("CUSTOM_SURFACE"), st); + EXPECT_EQ("CUSTOM_SURFACE", math::SphericalCoordinates::Convert(st)); } ////////////////////////////////////////////////// @@ -98,6 +113,14 @@ TEST(SphericalCoordinatesTest, SetFunctions) EXPECT_EQ(sc.LongitudeReference(), gz::math::Angle()); EXPECT_EQ(sc.HeadingOffset(), gz::math::Angle()); EXPECT_NEAR(sc.ElevationReference(), 0.0, 1e-6); + EXPECT_NEAR(sc.SurfaceRadius(), + 6371000.0, 1e-3); + EXPECT_NEAR(sc.SurfaceAxisEquatorial(), + 6378137.0, 1e-3); + EXPECT_NEAR(sc.SurfaceAxisPolar(), + 6356752.314245, 1e-3); + EXPECT_NEAR(sc.SurfaceFlattening(), + 1.0/298.257223563, 1e-5); { gz::math::Angle lat(0.3), lon(-1.2), heading(0.5); @@ -113,7 +136,71 @@ TEST(SphericalCoordinatesTest, SetFunctions) EXPECT_EQ(sc.LongitudeReference(), lon); EXPECT_EQ(sc.HeadingOffset(), heading); EXPECT_NEAR(sc.ElevationReference(), elev, 1e-6); + EXPECT_NEAR(sc.SurfaceRadius(), + 6371000.0, 1e-3); + EXPECT_NEAR(sc.SurfaceAxisEquatorial(), + 6378137.0, 1e-3); + EXPECT_NEAR(sc.SurfaceAxisPolar(), + 6356752.314245, 1e-3); + EXPECT_NEAR(sc.SurfaceFlattening(), + 1.0/298.257223563, 1e-5); } + + // Moon surface type + st = math::SphericalCoordinates::MOON_SCS; + math::SphericalCoordinates moonSC(st); + moonSC.SetSurface(st); + EXPECT_EQ(moonSC.Surface(), st); + EXPECT_NEAR(moonSC.SurfaceRadius(), + 1737400.0, 1e-3); + EXPECT_NEAR(moonSC.SurfaceAxisEquatorial(), + 1738100.0, 1e-3); + EXPECT_NEAR(moonSC.SurfaceAxisPolar(), + 1736000.0, 1e-3); + EXPECT_NEAR(moonSC.SurfaceFlattening(), + 0.0012, 1e-5); +} + +////////////////////////////////////////////////// +/// Test invalid parameters for custom surface +TEST(SphericalCoordinatesTest, InvalidParameters) +{ + // Earth's constants + double g_EarthWGS84AxisEquatorial = 6378137.0; + double g_EarthWGS84AxisPolar = 6356752.314245; + double g_EarthWGS84Flattening = 1.0/298.257223563; + double g_EarthRadius = 6371000.0; + + // Create a custom surface with invalid parameters. + math::SphericalCoordinates scInvalid( + math::SphericalCoordinates::CUSTOM_SURFACE, + -1, -1); + + // These should be rejected and default to Earth's + // parameters. + EXPECT_NEAR(scInvalid.SurfaceRadius(), g_EarthRadius, + 1e-3); + EXPECT_NEAR(scInvalid.SurfaceAxisEquatorial(), + g_EarthWGS84AxisEquatorial, 1e-3); + EXPECT_NEAR(scInvalid.SurfaceAxisPolar(), + g_EarthWGS84AxisPolar, 1e-3); + EXPECT_NEAR(scInvalid.SurfaceFlattening(), + g_EarthWGS84Flattening, 1e-3); + + // Create a custom surface with valid parameters. + math::SphericalCoordinates scValid( + math::SphericalCoordinates::CUSTOM_SURFACE, + 100, 100); + + // These should be accepted + EXPECT_NEAR(scValid.SurfaceRadius(), 100, + 1e-3); + EXPECT_NEAR(scValid.SurfaceAxisEquatorial(), + 100, 1e-3); + EXPECT_NEAR(scValid.SurfaceAxisPolar(), + 100, 1e-3); + EXPECT_NEAR(scValid.SurfaceFlattening(), + 0, 1e-3); } ////////////////////////////////////////////////// @@ -311,17 +398,51 @@ TEST(SphericalCoordinatesTest, Distance) longA.SetDegree(-122.249972); latB.SetDegree(46.124953); longB.SetDegree(-122.251683); - double d = math::SphericalCoordinates::Distance(latA, longA, latB, longB); - EXPECT_NEAR(14002, d, 20); + // Calculating distance using the static method. + double d1 = math::SphericalCoordinates::DistanceWGS84( + latA, longA, latB, longB); + EXPECT_NEAR(14002, d1, 20); + + // Using the non static method. The default surface type is EARTH_WGS84. + auto earthSC = math::SphericalCoordinates(); + double d2 = earthSC.DistanceBetweenPoints(latA, longA, latB, longB); + EXPECT_NEAR(d1, d2, 0.1); + + earthSC = math::SphericalCoordinates( + math::SphericalCoordinates::SurfaceType::EARTH_WGS84); + double d3 = earthSC.DistanceBetweenPoints(latA, longA, latB, longB); + EXPECT_NEAR(d2, d3, 0.1); + + // Setting the surface type as Moon. + auto moonSC = math::SphericalCoordinates( + math::SphericalCoordinates::SurfaceType::MOON_SCS); + double d4 = moonSC.DistanceBetweenPoints(latA, longA, latB, longB); + EXPECT_NEAR(3820, d4, 5); + + // Using a custom surface. + // For custom surfaces, the surface properties need to be set. + // This one will throw an error. + auto invalidCustomSC = math::SphericalCoordinates( + math::SphericalCoordinates::CUSTOM_SURFACE); + // This one should be accepted. + auto customSC = math::SphericalCoordinates( + math::SphericalCoordinates::SurfaceType::CUSTOM_SURFACE, + 6378137.0, + 6356752.314245); + + EXPECT_NEAR(customSC.DistanceBetweenPoints(latA, longA, latB, longB), + d1, 0.1); } ////////////////////////////////////////////////// TEST(SphericalCoordinatesTest, BadSetSurface) { math::SphericalCoordinates sc; - sc.SetSurface(static_cast(2)); - EXPECT_EQ(sc.Surface(), 2); + sc.SetSurface(static_cast(3), + 10, 10); + sc.SetSurface(static_cast(3)); + EXPECT_EQ(sc.Surface(), 3); } ////////////////////////////////////////////////// diff --git a/src/python_pybind11/src/SphericalCoordinates.cc b/src/python_pybind11/src/SphericalCoordinates.cc index 69e5046e5..f2d27ed2e 100644 --- a/src/python_pybind11/src/SphericalCoordinates.cc +++ b/src/python_pybind11/src/SphericalCoordinates.cc @@ -42,6 +42,8 @@ void defineMathSphericalCoordinates(py::module &m, const std::string &typestr) .def(py::init()) + .def(py::init()) .def(py::self != py::self) .def(py::self == py::self) .def("spherical_from_local_position", @@ -57,8 +59,15 @@ void defineMathSphericalCoordinates(py::module &m, const std::string &typestr) .def("convert", py::overload_cast(&Class::Convert), "Convert a SurfaceType to a string.") - .def("distance", - &Class::Distance, + .def("distance_WGS84", + &Class::DistanceWGS84, + "Get the distance between two points expressed in geographic " + "latitude and longitude. It assumes that both points are at sea level." + " Example: _latA = 38.0016667 and _lonA = -123.0016667) represents " + "the point with latitude 38d 0'6.00\"N and longitude 123d 0'6.00\"W." + " This function assumes the surface is EARTH_WGS84.") + .def("distance_between_points", + &Class::DistanceBetweenPoints, "Get the distance between two points expressed in geographic " "latitude and longitude. It assumes that both points are at sea level." " Example: _latA = 38.0016667 and _lonA = -123.0016667) represents " @@ -66,6 +75,18 @@ void defineMathSphericalCoordinates(py::module &m, const std::string &typestr) .def("surface", &Class::Surface, "Get SurfaceType currently in use.") + .def("surface_radius", + &Class::SurfaceRadius, + "Get the radius of the surface.") + .def("surface_axis_equatorial", + &Class::SurfaceAxisEquatorial, + "Get the major of the surface.") + .def("surface_axis_polar", + &Class::SurfaceAxisPolar, + "Get the minor axis of the surface.") + .def("surface_flattening", + &Class::SurfaceFlattening, + "Get the flattening parameter of the surface.") .def("latitude_reference", &Class::LatitudeReference, "Get reference geodetic latitude.") @@ -81,7 +102,12 @@ void defineMathSphericalCoordinates(py::module &m, const std::string &typestr) "angle from East to x-axis, or equivalently " "from North to y-axis.") .def("set_surface", - &Class::SetSurface, + py::overload_cast(&Class::SetSurface), + "Set SurfaceType for planetary surface model.") + .def("set_surface", + py::overload_cast(&Class::SetSurface), "Set SurfaceType for planetary surface model.") .def("set_latitude_reference", &Class::SetLatitudeReference, @@ -125,6 +151,8 @@ void defineMathSphericalCoordinates(py::module &m, const std::string &typestr) .export_values(); py::enum_(sphericalCoordinates, "SurfaceType") .value("EARTH_WGS84", Class::SurfaceType::EARTH_WGS84) + .value("MOON_SCS", Class::SurfaceType::MOON_SCS) + .value("CUSTOM_SURFACE", Class::SurfaceType::CUSTOM_SURFACE) .export_values(); } } // namespace python diff --git a/src/python_pybind11/test/SphericalCoordinates_TEST.py b/src/python_pybind11/test/SphericalCoordinates_TEST.py index 50a52ef96..72f8d8a53 100644 --- a/src/python_pybind11/test/SphericalCoordinates_TEST.py +++ b/src/python_pybind11/test/SphericalCoordinates_TEST.py @@ -67,6 +67,18 @@ def test_convert(self): self.assertEqual("EARTH_WGS84", SphericalCoordinates.convert(st)) + # For Moon surface type + st = SphericalCoordinates.MOON_SCS + self.assertEqual(SphericalCoordinates.convert("MOON_SCS"), st) + self.assertEqual("MOON_SCS", SphericalCoordinates.convert(st)) + + # For custom surface type + st = SphericalCoordinates.CUSTOM_SURFACE + self.assertEqual(SphericalCoordinates.convert("CUSTOM_SURFACE"), + st) + self.assertEqual("CUSTOM_SURFACE", + SphericalCoordinates.convert(st)) + def test_set_functions(self): # Default surface type st = SphericalCoordinates.EARTH_WGS84 @@ -78,6 +90,13 @@ def test_set_functions(self): self.assertEqual(sc.longitude_reference(), Angle()) self.assertEqual(sc.heading_offset(), Angle()) self.assertAlmostEqual(sc.elevation_reference(), 0.0, delta=1e-6) + self.assertAlmostEqual(sc.surface_radius(), 6371000.0, delta=1e-3) + self.assertAlmostEqual(sc.surface_axis_equatorial(), + 6378137.0, delta=1e-3) + self.assertAlmostEqual(sc.surface_axis_polar(), + 6356752.314245, delta=1e-3) + self.assertAlmostEqual(sc.surface_flattening(), + 1.0/298.257223563, delta=1e-5) lat = Angle(0.3) lon = Angle(-1.2) @@ -94,6 +113,63 @@ def test_set_functions(self): self.assertEqual(sc.longitude_reference(), lon) self.assertEqual(sc.heading_offset(), heading) self.assertAlmostEqual(sc.elevation_reference(), elev, delta=1e-6) + self.assertAlmostEqual(sc.surface_radius(), 6371000.0, delta=1e-3) + self.assertAlmostEqual(sc.surface_axis_equatorial(), + 6378137.0, delta=1e-3) + self.assertAlmostEqual(sc.surface_axis_polar(), + 6356752.314245, delta=1e-3) + self.assertAlmostEqual(sc.surface_flattening(), + 1.0/298.257223563, delta=1e-5) + + # Moon surface type + st = SphericalCoordinates.MOON_SCS + sc = SphericalCoordinates(st) + sc.set_surface(st) + self.assertAlmostEqual(sc.surface_radius(), 1737400.0, + delta=1e-3) + self.assertAlmostEqual(sc.surface_axis_equatorial(), + 1738100.0, delta=1e-3) + self.assertAlmostEqual(sc.surface_axis_polar(), + 1736000.0, delta=1e-3) + self.assertAlmostEqual(sc.surface_flattening(), + 0.0012, delta=1e-5) + + def test_invalid_parameters(self): + # Earth's constants + g_EarthWGS84AxisEquatorial = 6378137.0; + g_EarthWGS84AxisPolar = 6356752.314245; + g_EarthWGS84Flattening = 1.0/298.257223563; + g_EarthRadius = 6371000.0; + + # Create a custom surface with invalid parameters + sc_invalid = SphericalCoordinates( + SphericalCoordinates.CUSTOM_SURFACE, + -1, -1) + + # These should be rejected and default to Earth's parameters + self.assertAlmostEqual(sc_invalid.surface_radius(), + g_EarthRadius, delta=1e-3) + self.assertAlmostEqual(sc_invalid.surface_axis_equatorial(), + g_EarthWGS84AxisEquatorial, delta=1e-3) + self.assertAlmostEqual(sc_invalid.surface_axis_polar(), + g_EarthWGS84AxisPolar, delta=1e-3) + self.assertAlmostEqual(sc_invalid.surface_flattening(), + g_EarthWGS84Flattening, delta=1e-3) + + # Creating a custom surface with valid parameters + sc_valid = SphericalCoordinates( + SphericalCoordinates.CUSTOM_SURFACE, + 100, 100) + + # These should be accepted + self.assertAlmostEqual(sc_valid.surface_radius(), + 100, delta=1e-3) + self.assertAlmostEqual(sc_valid.surface_axis_equatorial(), + 100, delta=1e-3) + self.assertAlmostEqual(sc_valid.surface_axis_polar(), + 100, delta=1e-3) + self.assertAlmostEqual(sc_valid.surface_flattening(), + 0, delta=1e-3) def test_coordinate_transforms(self): # Default surface type @@ -263,14 +339,42 @@ def test_distance(self): longA.set_degree(-122.249972) latB.set_degree(46.124953) longB.set_degree(-122.251683) - d = SphericalCoordinates.distance(latA, longA, latB, longB) - - self.assertAlmostEqual(14002, d, delta=20) + d1 = SphericalCoordinates.distance_WGS84(latA, longA, latB, longB) + + self.assertAlmostEqual(14002, d1, delta=20) + + # Using the non static method. The default surface is EARTH_WGS84. + earth_sc = SphericalCoordinates() + d2 = earth_sc.distance_between_points(latA, longA, latB, longB) + self.assertAlmostEqual(d1, d2, delta=0.1) + + earth_sc = SphericalCoordinates(SphericalCoordinates.EARTH_WGS84) + d3 = earth_sc.distance_between_points(latA, longA, latB, longB) + self.assertAlmostEqual(d2, d3, delta=0.1) + + # Using the surface type as Moon. + moon_sc = SphericalCoordinates(SphericalCoordinates.MOON_SCS) + d4 = moon_sc.distance_between_points(latA, longA, latB, longB) + self.assertAlmostEqual(3820, d4, delta=5) + + # Using a custom surface + # For custom surfaces, the surface properties need to be set. + # This line should throw an error. + invalid_custom_sc = SphericalCoordinates( + SphericalCoordinates.CUSTOM_SURFACE) + # This one should be accepted. + valid_custom_sc = SphericalCoordinates( + SphericalCoordinates.CUSTOM_SURFACE, + 6378137.0, + 6356752.314245); + + self.assertAlmostEqual(valid_custom_sc.distance_between_points(latA, longA, latB, longB), + d1, delta=0.1) def test_bad_set_surface(self): sc = SphericalCoordinates() - sc.set_surface(SphericalCoordinates.SurfaceType(2)) - self.assertEqual(sc.surface(), SphericalCoordinates.SurfaceType(2)) + sc.set_surface(SphericalCoordinates.SurfaceType(3)) + self.assertEqual(sc.surface(), SphericalCoordinates.SurfaceType(3)) def test_transform(self): sc = SphericalCoordinates()