Make greatCircleCourse private before moving to eckit

ecmwf · Dec 12, 2023 · 1212789 · 1212789
1 parent 09af7e7
commit 1212789
Show file tree

Hide file tree

Showing 5 changed files with 72 additions and 42 deletions.
diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc
@@ -26,7 +26,7 @@
 #include "atlas/runtime/Exception.h"
 #include "atlas/runtime/Trace.h"
 #include "atlas/util/Constants.h"
-#include "atlas/util/UnitSphere.h"
+#include "atlas/util/Geometry.h"
 
 #include "eckit/linalg/Triplet.h"
 
@@ -161,13 +161,15 @@ void SphericalVector::do_setup(const FunctionSpace& source,
   const auto sourceLonLats = array::make_view<double, 2>(source_.lonlat());
   const auto targetLonLats = array::make_view<double, 2>(target_.lonlat());
 
+  geometry::UnitSphere unitSphere;
+
   const auto setWeights = [&](auto i, auto j, const auto& baseWeight) {
     const auto sourceLonLat =
         PointLonLat(sourceLonLats(j, 0), sourceLonLats(j, 1));
     const auto targetLonLat =
         PointLonLat(targetLonLats(i, 0), targetLonLats(i, 1));
 
-    const auto alpha = util::greatCircleCourse(sourceLonLat, targetLonLat);
+    const auto alpha = unitSphere.greatCircleCourse(sourceLonLat, targetLonLat);
 
     const auto deltaAlpha =
         (alpha.first - alpha.second) * util::Constants::degreesToRadians();

diff --git a/src/atlas/util/Geometry.cc b/src/atlas/util/Geometry.cc
@@ -8,12 +8,14 @@
  * nor does it submit to any jurisdiction.
  */
 
+#include "atlas/util/Geometry.h"
+
 #include "eckit/geometry/Point2.h"
 #include "eckit/geometry/Point3.h"
 
 #include "atlas/library/config.h"
 #include "atlas/runtime/Exception.h"
-#include "atlas/util/Geometry.h"
+#include "atlas/util/Constants.h"
 
 namespace atlas {
 
@@ -30,6 +32,44 @@ void GeometrySphere::xyz2lonlat(const Point3& xyz, Point2& lonlat) const {
     Sphere::convertCartesianToSpherical(radius_, xyz, lonlat);
 }
 
+
+/// @brief   Calculate great-cricle course between points
+///
+/// @details Calculates the direction (clockwise from north) of a great-circle
+///          arc between lonLat1 and lonLat2. Returns the direction of the arc
+///          at lonLat1 (first) and lonLat2 (second). Angle is normalised to the
+///          range of atan2 (usually (-180, 180]). All input and output values
+///          are in units of degrees.
+/// @ref     https://en.wikipedia.org/wiki/Great-circle_navigation
+///
+std::pair<double, double> greatCircleCourse(const Point2& lonLat1,
+                                                   const Point2& lonLat2) {
+
+  const auto lambda1 = lonLat1[0] * util::Constants::degreesToRadians();
+  const auto lambda2 = lonLat2[0] * util::Constants::degreesToRadians();
+  const auto phi1 = lonLat1[1] * util::Constants::degreesToRadians();
+  const auto phi2 = lonLat2[1] * util::Constants::degreesToRadians();
+
+  const auto sinLambda12 = std::sin(lambda2 - lambda1);
+  const auto cosLambda12 = std::cos(lambda2 - lambda1);
+  const auto sinPhi1 = std::sin(phi1);
+  const auto sinPhi2 = std::sin(phi2);
+  const auto cosPhi1 = std::cos(phi1);
+  const auto cosPhi2 = std::cos(phi2);
+
+  const auto alpha1 =
+      std::atan2(cosPhi2 * sinLambda12,
+                 cosPhi1 * sinPhi2 - sinPhi1 * cosPhi2 * cosLambda12);
+
+  const auto alpha2 =
+      std::atan2(cosPhi1 * sinLambda12,
+                 -cosPhi2 * sinPhi1 + sinPhi2 * cosPhi1 * cosLambda12);
+
+  return std::make_pair(alpha1 * util::Constants::radiansToDegrees(),
+                        alpha2 * util::Constants::radiansToDegrees());
+};
+
+
 } // namespace detail
 } // namespace geometry
 

diff --git a/src/atlas/util/Geometry.h b/src/atlas/util/Geometry.h
@@ -27,6 +27,20 @@ namespace geometry {
 
 namespace detail {
 
+// TODO: move greatCircleCourse to eckit::geometry::Sphere
+
+/// @brief   Calculate great-cricle course between points
+///
+/// @details Calculates the direction (clockwise from north) of a great-circle
+///          arc between lonLat1 and lonLat2. Returns the direction of the arc
+///          at lonLat1 (first) and lonLat2 (second). Angle is normalised to the
+///          range of atan2 (usually (-180, 180]). All input and output values
+///          are in units of degrees.
+/// @ref     https://en.wikipedia.org/wiki/Great-circle_navigation
+std::pair<double, double> greatCircleCourse(const Point2& lonLat1, const Point2& lonLat2);
+
+//------------------------------------------------------------------------------------------------------
+
 class GeometryBase : public util::Object {
 public:
     virtual ~GeometryBase()                                           = default;
@@ -36,6 +50,7 @@ class GeometryBase : public util::Object {
     virtual double distance(const Point3& p1, const Point3& p2) const = 0;
     virtual double radius() const                                     = 0;
     virtual double area() const                                       = 0;
+    virtual std::pair<double, double> greatCircleCourse(const Point2& p1, const Point2& p2) = 0;
 
     Point3 xyz(const Point2& lonlat) const {
         Point3 xyz;
@@ -64,6 +79,9 @@ class GeometrySphereT : public GeometryBase {
     double distance(const Point3& p1, const Point3& p2) const override { return SphereT::distance(p1, p2); }
     double radius() const override { return SphereT::radius(); }
     double area() const override { return SphereT::area(); }
+    std::pair<double, double> greatCircleCourse(const Point2& p1, const Point2& p2) override {
+        return atlas::geometry::detail::greatCircleCourse(p1, p2);
+    }
 };
 
 class GeometrySphere : public GeometryBase {
@@ -77,6 +95,9 @@ class GeometrySphere : public GeometryBase {
     double distance(const Point3& p1, const Point3& p2) const override { return Sphere::distance(radius_, p1, p2); }
     double radius() const override { return radius_; }
     double area() const override { return Sphere::area(radius_); }
+    std::pair<double, double> greatCircleCourse(const Point2& p1, const Point2& p2) override {
+        return atlas::geometry::detail::greatCircleCourse(p1, p2);
+    }
 
 private:
     double radius_;
@@ -107,6 +128,7 @@ class Geometry : public util::ObjectHandle<geometry::detail::GeometryBase> {
     double distance(const Point3& p1, const Point3& p2) const { return get()->distance(p1, p2); }
     double radius() const { return get()->radius(); }
     double area() const { return get()->area(); }
+    std::pair<double, double> greatCircleCourse(const Point2& p1, const Point2& p2) { return get()->greatCircleCourse(p1, p2); }
 
 protected:
     template <typename GeometryT, typename... Args>

diff --git a/src/atlas/util/UnitSphere.h b/src/atlas/util/UnitSphere.h
@@ -27,42 +27,6 @@ namespace util {
 
 using eckit::geometry::UnitSphere;
 
-/// @brief   Calculate great-cricle course between points
-///
-/// @details Calculates the direction (clockwise from north) of a great-circle
-///          arc between lonLat1 and lonLat2. Returns the direction of the arc
-///          at lonLat1 (first) and lonLat2 (second). Angle is normalised to the
-///          range of atan2 (usually (-180, 180]). All input and output values
-///          are in units of degrees.
-/// @ref     https://en.wikipedia.org/wiki/Great-circle_navigation
-///
-inline std::pair<double, double> greatCircleCourse(const Point2& lonLat1,
-                                                   const Point2& lonLat2) {
-
-  const auto lambda1 = lonLat1[0] * Constants::degreesToRadians();
-  const auto lambda2 = lonLat2[0] * Constants::degreesToRadians();
-  const auto phi1 = lonLat1[1] * Constants::degreesToRadians();
-  const auto phi2 = lonLat2[1] * Constants::degreesToRadians();
-
-  const auto sinLambda12 = std::sin(lambda2 - lambda1);
-  const auto cosLambda12 = std::cos(lambda2 - lambda1);
-  const auto sinPhi1 = std::sin(phi1);
-  const auto sinPhi2 = std::sin(phi2);
-  const auto cosPhi1 = std::cos(phi1);
-  const auto cosPhi2 = std::cos(phi2);
-
-  const auto alpha1 =
-      std::atan2(cosPhi2 * sinLambda12,
-                 cosPhi1 * sinPhi2 - sinPhi1 * cosPhi2 * cosLambda12);
-
-  const auto alpha2 =
-      std::atan2(cosPhi1 * sinLambda12,
-                 -cosPhi2 * sinPhi1 + sinPhi2 * cosPhi1 * cosLambda12);
-
-  return std::make_pair(alpha1 * Constants::radiansToDegrees(),
-                        alpha2 * Constants::radiansToDegrees());
-};
-
 //------------------------------------------------------------------------------------------------------
 
 }  // namespace util

diff --git a/src/tests/util/test_unitsphere.cc b/src/tests/util/test_unitsphere.cc
@@ -6,7 +6,7 @@
  */
 
 #include "atlas/util/Point.h"
-#include "atlas/util/UnitSphere.h"
+#include "atlas/util/Geometry.h"
 
 #include "tests/AtlasTestEnvironment.h"
 
@@ -17,6 +17,8 @@ using namespace atlas::util;
 
 CASE("great-circle course") {
 
+  geometry::UnitSphere g;
+
   const auto point1 = PointLonLat(-71.6, -33.0);  // Valparaiso
   const auto point2 = PointLonLat(121.8, 31.4);   // Shanghai
   const auto point3 = PointLonLat(0., 89.);
@@ -27,8 +29,8 @@ CASE("great-circle course") {
   const auto targetCourse3 = 0.;
   const auto targetCourse4 = 180.;
 
-  const auto[ course1, course2 ] = greatCircleCourse(point1, point2);
-  const auto[ course3, course4 ] = greatCircleCourse(point3, point4);
+  const auto[ course1, course2 ] = g.greatCircleCourse(point1, point2);
+  const auto[ course3, course4 ] = g.greatCircleCourse(point3, point4);
 
   EXPECT_APPROX_EQ(course1, targetCourse1, 0.01);
   EXPECT_APPROX_EQ(course2, targetCourse2, 0.01);