Skip to content

Commit

Permalink
Update SphericalVector to work with StructuredColumns as source fun…
Browse files Browse the repository at this point in the history 
…ctionspace. (#175)

* Realised StructuredColumns weight corrections cancel due to Euler identity.

* Added great-circle course tests for points on the pole.

* Added great-circle course special case where lonlat1 == lonlat2.
  • Loading branch information
@odlomax
odlomax authored Mar 6, 2024
1 parent 4e2a055 commit c3333e9
Show file tree
Hide file tree
Showing 4 changed files with 162 additions and 165 deletions.
9 changes: 2 additions & 7 deletions src/atlas/interpolation/method/sphericalvector/SphericalVector.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,13 +78,8 @@ void SphericalVector::do_setup(const FunctionSpace& source,
option::variables(2));
auto targetLonLats = target_.createField<double>(option::name("lonlat") |
option::variables(2));
sourceLonLats.array().copy(source_.lonlat());
targetLonLats.array().copy(target_.lonlat());
sourceLonLats.haloExchange();
targetLonLats.haloExchange();

const auto sourceLonLatsView = array::make_view<double, 2>(sourceLonLats);
const auto targetLonLatsView = array::make_view<double, 2>(targetLonLats);
const auto sourceLonLatsView = array::make_view<double, 2>(source_.lonlat());
const auto targetLonLatsView = array::make_view<double, 2>(target_.lonlat());

const auto unitSphere = geometry::UnitSphere{};

Expand Down
116 changes: 60 additions & 56 deletions src/atlas/util/Geometry.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
*
* 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
* In applying this licence, ECMWF does not waive the privileges and immunities
* granted to it by virtue of its status as an intergovernmental organisation
* nor does it submit to any jurisdiction.
*/
Expand All @@ -12,9 +12,6 @@

#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"
Expand All @@ -25,16 +22,15 @@ 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 +41,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 +70,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
Loading

0 comments on commit c3333e9

Please sign in to comment.