Skip to content

Commit

Permalink
Added great-circle course special case where lonlat1 == lonlat2.
Browse files Browse the repository at this point in the history
  • Loading branch information
@odlomax
odlomax committed Mar 4, 2024
1 parent ad0c800 commit 6a3f0c9
Showing 1 changed file with 64 additions and 58 deletions.
122 changes: 64 additions & 58 deletions src/atlas/util/Geometry.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,38 +3,36 @@
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
* In applying this licence, ECMWF does not waive the privileges and immGeometryies
* granted to it by virtue of its status as an intergovernmental organisation
* nor does it submit to any jurisdiction.
* In applying this licence, ECMWF does not waive the privileges and
* immGeometryies granted to it by virtue of its status as an intergovernmental
* organisation nor does it submit to any jurisdiction.
*/

#include "atlas/util/Geometry.h"

#include <cmath>

#include "eckit/geometry/Point2.h"
#include "eckit/geometry/Point3.h"

#include "atlas/library/config.h"
#include "atlas/runtime/Exception.h"
#include "atlas/util/Constants.h"
#include "eckit/geometry/Point2.h"
#include "eckit/geometry/Point3.h"

namespace atlas {

namespace geometry {
namespace detail {
void GeometrySphere::lonlat2xyz(const Point2& lonlat, Point3& xyz) const {
#if ATLAS_ECKIT_VERSION_AT_LEAST(1, 24, 0)
Sphere::convertSphericalToCartesian(radius_, lonlat, xyz, 0., true);
Sphere::convertSphericalToCartesian(radius_, lonlat, xyz, 0., true);
#else
Sphere::convertSphericalToCartesian(radius_, lonlat, xyz);
Sphere::convertSphericalToCartesian(radius_, lonlat, xyz);
#endif
}
void GeometrySphere::xyz2lonlat(const Point3& xyz, Point2& lonlat) const {
Sphere::convertCartesianToSpherical(radius_, xyz, lonlat);
Sphere::convertCartesianToSpherical(radius_, xyz, lonlat);
}


/// @brief Calculate great-cricle course between points
///
/// @details Calculates the direction (clockwise from north) of a great-circle
Expand All @@ -45,7 +43,10 @@ void GeometrySphere::xyz2lonlat(const Point3& xyz, Point2& lonlat) const {
/// @ref https://en.wikipedia.org/wiki/Great-circle_navigation
///
std::pair<double, double> greatCircleCourse(const Point2& lonLat1,
const Point2& lonLat2) {
const Point2& lonLat2) {
if (lonLat1 == lonLat2) {
return std::make_pair(0., 0.);
}

const auto lambda1 = lonLat1[0] * util::Constants::degreesToRadians();
const auto lambda2 = lonLat2[0] * util::Constants::degreesToRadians();
Expand All @@ -71,71 +72,76 @@ std::pair<double, double> greatCircleCourse(const Point2& lonLat1,
alpha2 * util::Constants::radiansToDegrees());
};


} // namespace detail
} // namespace geometry
} // namespace detail
} // namespace geometry

extern "C" {
// ------------------------------------------------------------------
// C wrapper interfaces to C++ routines
Geometry::Implementation* atlas__Geometry__new_name(const char* name) {
Geometry::Implementation* geometry;
{
Geometry handle{std::string{name}};
geometry = handle.get();
geometry->attach();
}
geometry->detach();
return geometry;
Geometry::Implementation* geometry;
{
Geometry handle{std::string{name}};
geometry = handle.get();
geometry->attach();
}
geometry->detach();
return geometry;
}
geometry::detail::GeometryBase* atlas__Geometry__new_radius(const double radius) {
Geometry::Implementation* geometry;
{
Geometry handle{radius};
geometry = handle.get();
geometry->attach();
}
geometry->detach();
return geometry;
geometry::detail::GeometryBase* atlas__Geometry__new_radius(
const double radius) {
Geometry::Implementation* geometry;
{
Geometry handle{radius};
geometry = handle.get();
geometry->attach();
}
geometry->detach();
return geometry;
}
void atlas__Geometry__delete(Geometry::Implementation* This) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
delete This;
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
delete This;
}
void atlas__Geometry__xyz2lonlat(Geometry::Implementation* This, const double x, const double y, const double z,
double& lon, double& lat) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
PointLonLat lonlat;
This->xyz2lonlat(PointXYZ{x, y, z}, lonlat);
lon = lonlat.lon();
lat = lonlat.lat();
void atlas__Geometry__xyz2lonlat(Geometry::Implementation* This, const double x,
const double y, const double z, double& lon,
double& lat) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
PointLonLat lonlat;
This->xyz2lonlat(PointXYZ{x, y, z}, lonlat);
lon = lonlat.lon();
lat = lonlat.lat();
}
void atlas__Geometry__lonlat2xyz(Geometry::Implementation* This, const double lon, const double lat, double& x,
void atlas__Geometry__lonlat2xyz(Geometry::Implementation* This,
const double lon, const double lat, double& x,
double& y, double& z) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
PointXYZ xyz;
This->lonlat2xyz(PointLonLat{lon, lat}, xyz);
x = xyz.x();
y = xyz.y();
z = xyz.z();
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
PointXYZ xyz;
This->lonlat2xyz(PointLonLat{lon, lat}, xyz);
x = xyz.x();
y = xyz.y();
z = xyz.z();
}
double atlas__Geometry__distance_lonlat(Geometry::Implementation* This, const double lon1, const double lat1,
double atlas__Geometry__distance_lonlat(Geometry::Implementation* This,
const double lon1, const double lat1,
const double lon2, const double lat2) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
return This->distance(PointLonLat{lon1, lat1}, PointLonLat{lon2, lat2});
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
return This->distance(PointLonLat{lon1, lat1}, PointLonLat{lon2, lat2});
}
double atlas__Geometry__distance_xyz(Geometry::Implementation* This, const double x1, const double y1, const double z1,
const double x2, const double y2, const double z2) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
return This->distance(PointXYZ{x1, y1, z1}, PointXYZ{x2, y2, z2});
double atlas__Geometry__distance_xyz(Geometry::Implementation* This,
const double x1, const double y1,
const double z1, const double x2,
const double y2, const double z2) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
return This->distance(PointXYZ{x1, y1, z1}, PointXYZ{x2, y2, z2});
}
double atlas__Geometry__radius(Geometry::Implementation* This) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
return This->radius();
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
return This->radius();
}
double atlas__Geometry__area(Geometry::Implementation* This) {
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
return This->area();
ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry");
return This->area();
}
}
// ------------------------------------------------------------------
Expand Down

0 comments on commit 6a3f0c9

Please sign in to comment.