diff --git a/src/atlas/projection/detail/ProjectionImpl.h b/src/atlas/projection/detail/ProjectionImpl.h index 655134b11..810dd99f1 100644 --- a/src/atlas/projection/detail/ProjectionImpl.h +++ b/src/atlas/projection/detail/ProjectionImpl.h @@ -10,9 +10,12 @@ #pragma once +#include +#include #include #include +#include "atlas/runtime/Exception.h" #include "atlas/util/Factory.h" #include "atlas/util/NormaliseLongitude.h" #include "atlas/util/Object.h" @@ -41,97 +44,97 @@ class ProjectionImpl : public util::Object { using Spec = atlas::util::Config; class Jacobian : public std::array, 2> { public: + + using std::array, 2>::array; + static Jacobian identity() { - Jacobian id; - id[0] = {1.0, 0.0}; - id[1] = {0.0, 1.0}; - return id; + return Jacobian{1., 0., 0., 1.}; } - Jacobian inverse() const { - const Jacobian& jac = *this; - Jacobian inv; - double det = jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0]; - inv[0][0] = +jac[1][1] / det; - inv[0][1] = -jac[0][1] / det; - inv[1][0] = -jac[1][0] / det; - inv[1][1] = +jac[0][0] / det; - return inv; - }; + Jacobian(double j00, double j01, double j10, double j11) : + array{std::array{j00, j01}, + std::array{j10, j11}} {} - Jacobian transpose() const { - Jacobian tra = *this; - std::swap(tra[0][1], tra[1][0]); - return tra; - }; + Jacobian(std::initializer_list> list) : + Jacobian{*(list.begin()->begin()), + *(list.begin()->begin() + 1), + *((list.begin() + 1)->begin()), + *((list.begin() + 1)->begin() + 1)} {} - Jacobian operator-(const Jacobian& jac2) const { - const Jacobian& jac1 = *this; - Jacobian jac; - jac[0][0] = jac1[0][0] - jac2[0][0]; - jac[0][1] = jac1[0][1] - jac2[0][1]; - jac[1][0] = jac1[1][0] - jac2[1][0]; - jac[1][1] = jac1[1][1] - jac2[1][1]; - return jac; + Jacobian operator-(const Jacobian& jac) const { + return Jacobian{(*this)[0][0] - jac[0][0], (*this)[0][1] - jac[0][1], + (*this)[1][0] - jac[1][0], (*this)[1][1] - jac[1][1]}; } - Jacobian operator+(const Jacobian& jac2) const { - const Jacobian& jac1 = *this; - Jacobian jac; - jac[0][0] = jac1[0][0] + jac2[0][0]; - jac[0][1] = jac1[0][1] + jac2[0][1]; - jac[1][0] = jac1[1][0] + jac2[1][0]; - jac[1][1] = jac1[1][1] + jac2[1][1]; - return jac; + Jacobian operator+(const Jacobian& jac) const { + return Jacobian{(*this)[0][0] + jac[0][0], (*this)[0][1] + jac[0][1], + (*this)[1][0] + jac[1][0], (*this)[1][1] + jac[1][1]}; + } + + Jacobian operator*(double a) const { + return Jacobian{(*this)[0][0] * a, (*this)[0][1] * a, + (*this)[1][0] * a, (*this)[1][1] * a}; + } + + Jacobian operator*(const Jacobian& jac) const { + return Jacobian{(*this)[0][0] * jac[0][0] + (*this)[0][1] * jac[1][0], + (*this)[0][0] * jac[0][1] + (*this)[0][1] * jac[1][1], + (*this)[1][0] * jac[0][0] + (*this)[1][1] * jac[1][0], + (*this)[1][0] * jac[0][1] + (*this)[1][1] * jac[1][1]}; + } + + Point2 operator*(const Point2& x) const { + return Point2{(*this)[0][0] * x[0] + (*this)[0][1] * x[1], + (*this)[1][0] * x[0] + (*this)[1][1] * x[1]}; } double norm() const { - const Jacobian& jac = *this; - return sqrt(jac[0][0] * jac[0][0] + jac[0][1] * jac[0][1] + jac[1][0] * jac[1][0] + jac[1][1] * jac[1][1]); + return std::sqrt((*this)[0][0] * (*this)[0][0] + (*this)[0][1] * (*this)[0][1] + + (*this)[1][0] * (*this)[1][0] + (*this)[1][1] * (*this)[1][1]); + } + + double determinant() const { + return (*this)[0][0] * (*this)[1][1] - (*this)[0][1] * (*this)[1][0]; + } + + Jacobian inverse() const { + return Jacobian{(*this)[1][1], -(*this)[0][1], -(*this)[1][0], (*this)[0][0]} * (1. / determinant()); } - Jacobian operator*(const Jacobian& jac2) const { - const Jacobian& jac1 = *this; - Jacobian jac; - jac[0][0] = jac1[0][0] * jac2[0][0] + jac1[0][1] * jac2[1][0]; - jac[0][1] = jac1[0][0] * jac2[0][1] + jac1[0][1] * jac2[1][1]; - jac[1][0] = jac1[1][0] * jac2[0][0] + jac1[1][1] * jac2[1][0]; - jac[1][1] = jac1[1][0] * jac2[0][1] + jac1[1][1] * jac2[1][1]; - return jac; + Jacobian transpose() const { + return Jacobian{(*this)[1][1], (*this)[0][1], (*this)[1][0], (*this)[0][0]}; } double dx_dlon() const { - const Jacobian& jac = *this; - return jac[JDX][JDLON]; + return (*this)[JDX][JDLON]; } double dy_dlon() const { - const Jacobian& jac = *this; - return jac[JDY][JDLON]; + return (*this)[JDY][JDLON]; } double dx_dlat() const { - const Jacobian& jac = *this; - return jac[JDX][JDLAT]; + return (*this)[JDX][JDLAT]; } double dy_dlat() const { - const Jacobian& jac = *this; - return jac[JDY][JDLAT]; + return (*this)[JDY][JDLAT]; } double dlon_dx() const { - const Jacobian& jac = *this; - return jac[JDLON][JDX]; + return (*this)[JDLON][JDX]; } double dlon_dy() const { - const Jacobian& jac = *this; - return jac[JDLON][JDY]; + return (*this)[JDLON][JDY]; } double dlat_dx() const { - const Jacobian& jac = *this; - return jac[JDLAT][JDX]; + return (*this)[JDLAT][JDX]; } double dlat_dy() const { - const Jacobian& jac = *this; - return jac[JDLAT][JDY]; + return (*this)[JDLAT][JDY]; + } + + friend std::ostream& operator<<(std::ostream& os, const Jacobian& jac) { + os << jac[0][0] << " " << jac[0][1] << "\n" + << jac[1][0] << " " << jac[1][1]; + return os; } private: @@ -147,6 +150,7 @@ class ProjectionImpl : public util::Object { }; }; + public: static const ProjectionImpl* create(const eckit::Parametrisation& p); static const ProjectionImpl* create(const std::string& type, const eckit::Parametrisation& p); diff --git a/src/tests/projection/test_jacobian.cc b/src/tests/projection/test_jacobian.cc index 9206998ec..3b36fd619 100644 --- a/src/tests/projection/test_jacobian.cc +++ b/src/tests/projection/test_jacobian.cc @@ -163,6 +163,67 @@ CASE("test_lonlat") { doTaylorTest(StructuredGrid("N16"), 1e-9); } +CASE("test_constructors") { + + Projection::Jacobian a = {}; + Projection::Jacobian b = {1., 2., 3., 4.}; + Projection::Jacobian c = {{5., 6.}, + {7., 8.}}; + + Log::info() << "a =" << std::endl; + Log::info() << a << std::endl; + Log::info() << "b =" << std::endl; + Log::info() << b << std::endl; + Log::info() << "c =" << std::endl; + Log::info() << c << std::endl; + + EXPECT_EQ(a[0][0], 0.); + EXPECT_EQ(a[0][1], 0.); + EXPECT_EQ(a[1][0], 0.); + EXPECT_EQ(a[1][1], 0.); + + EXPECT_EQ(b[0][0], 1.); + EXPECT_EQ(b[0][1], 2.); + EXPECT_EQ(b[1][0], 3.); + EXPECT_EQ(b[1][1], 4.); + + EXPECT_EQ(c[0][0], 5.); + EXPECT_EQ(c[0][1], 6.); + EXPECT_EQ(c[1][0], 7.); + EXPECT_EQ(c[1][1], 8.); + +} + +CASE("test_multiplication") { + + const Projection::Jacobian A = {{0., -1.}, + {1., 0.}}; + const Point2 x = {2., 1.}; + + const auto Ax = A * x; + const auto AinvAx = A.inverse() * Ax; + const auto Atimes2 = A * 2.; + + Log::info() << "A =" << std::endl; + Log::info() << A << std::endl; + Log::info() << "2A =" << std::endl; + Log::info() << Atimes2 << std::endl; + Log::info() << "x =" << std::endl; + Log::info() << x << std::endl; + Log::info() << "Ax =" << std::endl; + Log::info() << Ax << std::endl; + Log::info() << "A^-1 Ax =" << std::endl; + Log::info() << AinvAx << std::endl; + + EXPECT_EQ(Atimes2[0][0], 0.); + EXPECT_EQ(Atimes2[0][1], -2.); + EXPECT_EQ(Atimes2[1][0], 2.); + EXPECT_EQ(Atimes2[1][1], 0.); + EXPECT_EQ(Ax, Point2(-1., 2.)); + EXPECT_EQ(AinvAx, x); + +} + //-----------------------------------------------------------------------------