diff --git a/src/core/observables/CylindricalLBProfileObservable.hpp b/src/core/observables/CylindricalLBProfileObservable.hpp index 5edd4b7ca7c..e2c27d1cf33 100644 --- a/src/core/observables/CylindricalLBProfileObservable.hpp +++ b/src/core/observables/CylindricalLBProfileObservable.hpp @@ -47,15 +47,14 @@ class CylindricalLBProfileObservable : public CylindricalProfileObservable { limits[0], limits[1], limits[2], n_bins[0], n_bins[1], n_bins[2], sampling_density); for (auto &p : sampling_positions) { - double theta; - Utils::Vector3d rotation_axis; auto p_cart = Utils::transform_coordinate_cylinder_to_cartesian( p, Utils::Vector3d{{0.0, 0.0, 1.0}}); // We have to rotate the coordinates since the utils function assumes // z-axis symmetry. - std::tie(theta, rotation_axis) = - Utils::rotation_params(Utils::Vector3d{{0.0, 0.0, 1.0}}, axis); - p_cart = Utils::vec_rotate(rotation_axis, theta, p_cart); + constexpr auto z_axis = Utils::Vector3d{{0.0, 0.0, 1.0}}; + auto const theta = angle_between(z_axis, axis); + auto const rot_axis = Utils::vector_product(z_axis, axis).normalize(); + p_cart = Utils::vec_rotate(rot_axis, theta, p_cart); p = p_cart + center; } } diff --git a/src/utils/include/utils/math/coordinate_transformation.hpp b/src/utils/include/utils/math/coordinate_transformation.hpp index 5316016e2de..09d17ff97c7 100644 --- a/src/utils/include/utils/math/coordinate_transformation.hpp +++ b/src/utils/include/utils/math/coordinate_transformation.hpp @@ -19,57 +19,168 @@ #ifndef UTILS_COORDINATE_TRANSFORMATION_HPP #define UTILS_COORDINATE_TRANSFORMATION_HPP +/** + * @file + * Convert coordinates from the Cartesian system to the cylindrical system. + * The transformation functions are provided with three overloads: + * - one function for the trivial Cartesian <-> cylindrical transformation + * - one function to transform from/to a cylindrical system with custom axis + * (extra @p axis argument, keep in mind the angle phi is under-defined) + * - one function to transform from/to an oriented cylindrical system with + * custom axis (extra @p orientation argument, the angle phi is well-defined) + */ + #include "utils/Vector.hpp" #include "utils/constants.hpp" +#include "utils/math/interval.hpp" #include "utils/math/vec_rotate.hpp" +#include + namespace Utils { -/** \brief Transform the given 3D position to cylinder coordinates with - * longitudinal axis aligned with axis parameter. +/** + * @brief Coordinate transformation from Cartesian to cylindrical coordinates. + * The origins and z-axis of the coordinate systems co-incide. + * The @f$ \phi = 0 @f$ direction corresponds to the x-axis in the + * original coordinate system. + * @param pos %Vector to transform */ inline Vector3d -transform_coordinate_cartesian_to_cylinder(const Vector3d &pos, - const Vector3d &axis) { +transform_coordinate_cartesian_to_cylinder(Vector3d const &pos) { + auto const r = std::sqrt(pos[0] * pos[0] + pos[1] * pos[1]); + auto const phi = std::atan2(pos[1], pos[0]); + return {r, phi, pos[2]}; +} + +/** + * @brief Coordinate transformation from Cartesian to cylindrical coordinates + * with change of basis. The origins of the coordinate systems co-incide. + * + * If the parameter @p axis is not equal to [0, 0, 1], the value + * of the angle @f$ \phi @f$ in cylindrical coordinates is under-defined. + * To fully define it, it is necessary to provide an orientation vector + * in Cartesian coordinates that will be used as the reference point + * (i.e. such that @f$ \phi = 0 @f$). + * + * If you don't need performance, you can call the overloaded function + * @ref transform_coordinate_cartesian_to_cylinder(Vector3d const &, Vector3d const &, Vector3d const &) directly. Otherwise, call this + * function once on the orientation vector and cache its @f$ \phi @f$ angle; + * this angle can be later passed as argument @p phi0 for a 50-80% speed-up. + * + * @param pos %Vector to transform + * @param axis Longitudinal axis of the cylindrical coordinates + * @param phi0 Angle offset to get @f$ \phi = 0 @f$ (obtained by calling + * this function on the orientation vector instead of @p pos) + */ +inline Vector3d transform_coordinate_cartesian_to_cylinder(Vector3d const &pos, + Vector3d const &axis, + double phi0 = 0.) { static auto const z_axis = Vector3d{{0, 0, 1}}; - double theta; - Vector3d rotation_axis; - std::tie(theta, rotation_axis) = rotation_params(axis, z_axis); - auto const rotated_pos = vec_rotate(rotation_axis, theta, pos); - auto const r = std::sqrt(rotated_pos[0] * rotated_pos[0] + - rotated_pos[1] * rotated_pos[1]); - auto const phi = std::atan2(rotated_pos[1], rotated_pos[0]); - return Vector3d{r, phi, rotated_pos[2]}; + auto const angle = angle_between(axis, z_axis); + auto const rotation_axis = Utils::vector_product(axis, z_axis).normalize(); + auto const pos_rotated = vec_rotate(rotation_axis, angle, pos); + auto pos_cylinder = transform_coordinate_cartesian_to_cylinder(pos_rotated); + if (phi0 != 0.) { + pos_cylinder[1] = interval(pos_cylinder[1] - phi0, -pi(), pi()); + } + return pos_cylinder; +} + +/** + * @brief Coordinate transformation from Cartesian to cylindrical coordinates + * with change of basis. The origins of the coordinate systems co-incide. + * @param pos %Vector to transform + * @param axis Longitudinal axis of the cylindrical coordinates + * @param orientation Reference point (in untransformed coordinates) for + * which @f$ \phi = 0 @f$ + */ +inline Vector3d transform_coordinate_cartesian_to_cylinder( + Vector3d const &pos, Vector3d const &axis, Vector3d const &orientation) { + auto const phi0 = + transform_coordinate_cartesian_to_cylinder(orientation, axis)[1]; + return transform_coordinate_cartesian_to_cylinder(pos, axis, phi0); } /** - * @brief Coordinate transformation from cylinder to cartesian coordinates. + * @brief Coordinate transformation from cylindrical to Cartesian coordinates. + * The origins and z-axis of the coordinate systems co-incide. + * The @f$ \phi = 0 @f$ direction corresponds to the x-axis in the + * transformed coordinate system. + * @param pos %Vector to transform */ inline Vector3d -transform_coordinate_cylinder_to_cartesian(Vector3d const &pos, - Vector3d const &axis) { - Vector3d const transformed{ - {pos[0] * std::cos(pos[1]), pos[0] * std::sin(pos[1]), pos[2]}}; +transform_coordinate_cylinder_to_cartesian(Vector3d const &pos) { + auto const &rho = pos[0]; + auto const &phi = pos[1]; + auto const &z = pos[2]; + return {rho * std::cos(phi), rho * std::sin(phi), z}; +} + +/** + * @brief Coordinate transformation from cylindrical to Cartesian coordinates + * with change of basis. The origins of the coordinate systems co-incide. + * + * If the parameter @p axis is not equal to [0, 0, 1], the value + * of the angle @f$ \phi @f$ in cylindrical coordinates is under-defined. + * To fully define it, it is necessary to provide an orientation vector + * in Cartesian coordinates that will be used as the reference point + * (i.e. such that @f$ \phi = 0 @f$). + * + * If you don't need performance, you can call the overloaded function + * @ref transform_coordinate_cylinder_to_cartesian(Vector3d const &, Vector3d const &, Vector3d const &) directly. Otherwise, call function + * @ref transform_coordinate_cartesian_to_cylinder(Vector3d const &, Vector3d const &, double) once on the orientation vector and cache + * its @f$ \phi @f$ angle; this angle can be later passed as argument + * @p phi0 for an 80-100% speed-up. + * + * @param pos %Vector to transform + * @param axis Longitudinal axis of the cylindrical coordinates + * @param phi0 Angle offset to get @f$ \phi = 0 @f$ + */ +inline Vector3d transform_coordinate_cylinder_to_cartesian(Vector3d const &pos, + Vector3d const &axis, + double phi0 = 0.) { static auto const z_axis = Vector3d{{0, 0, 1}}; - double theta; - Vector3d rotation_axis; - std::tie(theta, rotation_axis) = rotation_params(z_axis, axis); - auto const rotated_pos = vec_rotate(rotation_axis, theta, transformed); - return rotated_pos; + auto const angle = angle_between(z_axis, axis); + auto const rotation_axis = Utils::vector_product(z_axis, axis).normalize(); + auto const pos_cyl = Vector3d{{pos[0], pos[1] + phi0, pos[2]}}; + auto const pos_cart = transform_coordinate_cylinder_to_cartesian(pos_cyl); + auto const pos_rot = vec_rotate(rotation_axis, angle, pos_cart); + return pos_rot; } -/** \brief Transform the given 3D vector to cylinder coordinates with - * symmetry axis aligned with axis parameter. +/** + * @brief Coordinate transformation from cylindrical to Cartesian coordinates + * with change of basis. The origins of the coordinate systems co-incide. + * @param pos %Vector to transform + * @param axis Longitudinal axis of the cylindrical coordinates + * @param orientation Reference point (in Cartesian coordinates) for + * which @f$ \phi = 0 @f$ + */ +inline Vector3d transform_coordinate_cylinder_to_cartesian( + Vector3d const &pos, Vector3d const &axis, Vector3d const &orientation) { + auto const phi0 = + transform_coordinate_cartesian_to_cylinder(orientation, axis)[1]; + return transform_coordinate_cylinder_to_cartesian(pos, axis, phi0); +} + +/** + * @brief Vector transformation from Cartesian to cylindrical coordinates. + * @param vec %Vector to transform + * @param axis Longitudinal axis of the cylindrical coordinates + * @param pos Origin of the vector */ inline Vector3d transform_vector_cartesian_to_cylinder(Vector3d const &vec, Vector3d const &axis, Vector3d const &pos) { static auto const z_axis = Vector3d{{0, 0, 1}}; - double theta; - Vector3d rotation_axis; - std::tie(theta, rotation_axis) = rotation_params(axis, z_axis); - auto const rotated_pos = vec_rotate(rotation_axis, theta, pos); - auto const rotated_vec = vec_rotate(rotation_axis, theta, vec); + auto const angle = angle_between(axis, z_axis); + auto const rotation_axis = Utils::vector_product(axis, z_axis).normalize(); + auto const rotated_pos = vec_rotate(rotation_axis, angle, pos); + auto const rotated_vec = vec_rotate(rotation_axis, angle, vec); // v_r = (x * v_x + y * v_y) / sqrt(x^2 + y^2) auto const v_r = (rotated_pos[0] * rotated_vec[0] + rotated_pos[1] * rotated_vec[1]) / diff --git a/src/utils/include/utils/math/interval.hpp b/src/utils/include/utils/math/interval.hpp new file mode 100644 index 00000000000..b2ddf33ee58 --- /dev/null +++ b/src/utils/include/utils/math/interval.hpp @@ -0,0 +1,36 @@ +/* + * Copyright (C) 2020 The ESPResSo project + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ +#ifndef ESPRESSO_INTERVAL_HPP +#define ESPRESSO_INTERVAL_HPP + +#include + +namespace Utils { +/** + * @brief Wrap around the value of @p val in the interval [low, high]. + */ +template T interval(T val, T low, T high) { + if (val == low or val == high) + return val; + auto const span = high - low; + return std::fmod(span + std::fmod(val - low, span), span) + low; +} +} // namespace Utils + +#endif // ESPRESSO_INTERVAL_HPP diff --git a/src/utils/include/utils/math/vec_rotate.hpp b/src/utils/include/utils/math/vec_rotate.hpp index 4dbc624f168..559f297766f 100644 --- a/src/utils/include/utils/math/vec_rotate.hpp +++ b/src/utils/include/utils/math/vec_rotate.hpp @@ -23,7 +23,6 @@ #include "utils/math/sqr.hpp" #include -#include namespace Utils { /** @@ -50,17 +49,10 @@ inline Vector3d vec_rotate(const Vector3d &axis, double alpha, } /** - * @brief Determine rotation angle and axis for rotating vec onto target_vec. - * @param vec Vector to be rotated - * @param target_vec Target vector - * @return rotation angle and rotation axis + * @brief Determine the angle between two vectors. */ -inline std::tuple -rotation_params(Vector3d const &vec, Vector3d const &target_vec) { - auto const theta = - std::acos(vec * target_vec / (vec.norm() * target_vec.norm())); - auto const rotation_axis = Utils::vector_product(vec, target_vec).normalize(); - return std::make_tuple(theta, rotation_axis); +inline double angle_between(Vector3d const &v1, Vector3d const &v2) { + return std::acos(v1 * v2 / (v1.norm() * v2.norm())); } } // namespace Utils diff --git a/src/utils/tests/CMakeLists.txt b/src/utils/tests/CMakeLists.txt index 1bb9474a482..87caee85934 100644 --- a/src/utils/tests/CMakeLists.txt +++ b/src/utils/tests/CMakeLists.txt @@ -1,6 +1,7 @@ include(unit_test) unit_test(NAME abs_test SRC abs_test.cpp DEPENDS EspressoUtils) +unit_test(NAME interval_test SRC interval_test.cpp DEPENDS EspressoUtils) unit_test(NAME Vector_test SRC Vector_test.cpp DEPENDS EspressoUtils) unit_test(NAME Factory_test SRC Factory_test.cpp DEPENDS EspressoUtils) unit_test(NAME NumeratedContainer_test SRC NumeratedContainer_test.cpp DEPENDS diff --git a/src/utils/tests/coordinate_transformation.cpp b/src/utils/tests/coordinate_transformation.cpp index 714a11700e1..34327232123 100644 --- a/src/utils/tests/coordinate_transformation.cpp +++ b/src/utils/tests/coordinate_transformation.cpp @@ -25,10 +25,22 @@ #include #include +#include using Utils::Vector3d; BOOST_AUTO_TEST_CASE(cartesian_to_cylinder_test) { + constexpr auto eps = 1e-14; + auto const pos = Vector3d{{1.0, 3.3, 2.0}}; + auto const cyl = transform_coordinate_cartesian_to_cylinder(pos); + BOOST_CHECK_SMALL(abs(cyl[0] - std::sqrt(pos[0] * pos[0] + pos[1] * pos[1])), + eps); + BOOST_CHECK_SMALL(abs(cyl[1] - std::atan2(pos[1], pos[0])), eps); + BOOST_CHECK_SMALL(abs(cyl[2] - pos[2]), eps); +} + +BOOST_AUTO_TEST_CASE(cartesian_to_cylinder_with_axis_test) { + constexpr auto eps = 1e-14; Vector3d const cart_coord{{1.0, 3.3, 2.0}}; auto const transformed_x = transform_coordinate_cartesian_to_cylinder( cart_coord, Vector3d{{1, 0, 0}}); @@ -51,13 +63,96 @@ BOOST_AUTO_TEST_CASE(cartesian_to_cylinder_test) { std::atan2(cart_coord[1], cart_coord[0]), cart_coord[2]}}; for (int i = 0; i < 3; ++i) { - BOOST_CHECK(transformed_x[i] == expected_x[i]); - BOOST_CHECK(transformed_y[i] == expected_y[i]); - BOOST_CHECK(transformed_z[i] == expected_z[i]); + BOOST_CHECK_SMALL(abs(transformed_x[i] - expected_x[i]), eps); + BOOST_CHECK_SMALL(abs(transformed_y[i] - expected_y[i]), eps); + BOOST_CHECK_SMALL(abs(transformed_z[i] - expected_z[i]), eps); + } +} + +BOOST_AUTO_TEST_CASE(cartesian_to_cylinder_with_axis_with_phi_test) { + constexpr auto eps = 1e-14; + // tilted orthogonal basis + auto const x = + (Vector3d{{1, 0, 0}} - (1. / 3.) * Vector3d{{1, 1, 1}}).normalize(); + auto const y = (Vector3d{{0, 1, -1}}).normalize(); + auto const z = (Vector3d{{1, 1, 1}}).normalize(); + // check simple transformation without orientation (phi is random) + { + auto const x_cyl = transform_coordinate_cartesian_to_cylinder(x, z); + auto const y_cyl = transform_coordinate_cartesian_to_cylinder(y, z); + auto const z_cyl = transform_coordinate_cartesian_to_cylinder(z, z); + auto const x_ref = Vector3d{{1.0, x_cyl[1], 0.0}}; + auto const y_ref = Vector3d{{1.0, y_cyl[1], 0.0}}; + auto const z_ref = Vector3d{{0.0, z_cyl[1], 1.0}}; + for (int i = 0; i < 3; ++i) { + BOOST_CHECK_SMALL(abs(x_cyl[i] - x_ref[i]), eps); + BOOST_CHECK_SMALL(abs(y_cyl[i] - y_ref[i]), eps); + BOOST_CHECK_SMALL(abs(z_cyl[i] - z_ref[i]), eps); + } + } + // check transformation with orientation (phi is only random for r=0) + { + auto const x_cyl = transform_coordinate_cartesian_to_cylinder(x, z, y); + auto const y_cyl = transform_coordinate_cartesian_to_cylinder(y, z, y); + auto const z_cyl = transform_coordinate_cartesian_to_cylinder(z, z, y); + auto const x_ref = Vector3d{{1.0, -Utils::pi() / 2.0, 0.0}}; + auto const y_ref = Vector3d{{1.0, 0.0, 0.0}}; + auto const z_ref = Vector3d{{0.0, z_cyl[1], 1.0}}; + for (int i = 0; i < 3; ++i) { + BOOST_CHECK_SMALL(abs(x_cyl[i] - x_ref[i]), eps); + BOOST_CHECK_SMALL(abs(y_cyl[i] - y_ref[i]), eps); + BOOST_CHECK_SMALL(abs(z_cyl[i] - z_ref[i]), eps); + } + } + // check transformation with orientation for another angle + { + auto const u = vec_rotate(z, Utils::pi() / 3.0, x); + auto const v = vec_rotate(z, Utils::pi() / 3.0, y); + auto const u_cyl = transform_coordinate_cartesian_to_cylinder(u, z, y); + auto const v_cyl = transform_coordinate_cartesian_to_cylinder(v, z, y); + auto const u_ref = Vector3d{{1.0, Utils::pi() * (1. / 3. - 1. / 2.), 0.0}}; + auto const v_ref = Vector3d{{1.0, Utils::pi() / 3.0, 0.0}}; + for (int i = 0; i < 3; ++i) { + BOOST_CHECK_SMALL(abs(u_cyl[i] - u_ref[i]), eps); + BOOST_CHECK_SMALL(abs(v_cyl[i] - v_ref[i]), eps); + } + } + // check transformation of random vectors + { + std::subtract_with_carry_engine rng(2); + auto const r_uniform = [&rng]() { + return static_cast(rng() - rng.min()) / (rng.max() - rng.min()); + }; + for (int trial = 0; trial < 100; ++trial) { + Vector3d const v1{r_uniform(), r_uniform(), r_uniform()}; + Vector3d const v2{r_uniform(), r_uniform(), r_uniform()}; + auto const a = Utils::vector_product(v1, v2) / v1.norm() / v2.norm(); + auto const v1_v1 = transform_coordinate_cartesian_to_cylinder(v1, a, v1); + auto const v2_v1 = transform_coordinate_cartesian_to_cylinder(v2, a, v1); + auto const v1_v2 = transform_coordinate_cartesian_to_cylinder(v1, a, v2); + Vector3d const v1_v1_ref{v1.norm(), 0.0, 0.0}; + Vector3d const v2_v1_ref{v2.norm(), Utils::angle_between(v1, v2), 0.0}; + Vector3d const v1_v2_ref{v1.norm(), -Utils::angle_between(v1, v2), 0.0}; + for (int i = 0; i < 3; ++i) { + BOOST_CHECK_SMALL(abs(v1_v1[i] - v1_v1_ref[i]), eps); + BOOST_CHECK_SMALL(abs(v2_v1[i] - v2_v1_ref[i]), eps); + BOOST_CHECK_SMALL(abs(v1_v2[i] - v1_v2_ref[i]), eps); + } + } } } BOOST_AUTO_TEST_CASE(cylinder_to_cartesian_test) { + constexpr auto eps = 1e-14; + auto const cyl = Vector3d{{1.0, Utils::pi() / 4, 2.0}}; + auto const pos = transform_coordinate_cylinder_to_cartesian(cyl); + BOOST_CHECK_SMALL(abs(pos[0] - std::sqrt(2) / 2), eps); + BOOST_CHECK_SMALL(abs(pos[1] - std::sqrt(2) / 2), eps); + BOOST_CHECK_SMALL(abs(pos[2] - cyl[2]), eps); +} + +BOOST_AUTO_TEST_CASE(cylinder_to_cartesian_with_axis_test) { + constexpr auto eps = 1e-14; Vector3d const cylinder_coord{{1.2, 3.123, 42.0}}; auto const transformed_x = transform_coordinate_cylinder_to_cartesian( cylinder_coord, Vector3d{{1, 0, 0}}); @@ -80,8 +175,58 @@ BOOST_AUTO_TEST_CASE(cylinder_to_cartesian_test) { {cylinder_coord[0] * std::cos(cylinder_coord[1]), cylinder_coord[0] * std::sin(cylinder_coord[1]), cylinder_coord[2]}}; for (int i = 0; i < 3; ++i) { - BOOST_CHECK(transformed_x[i] == expected_x[i]); - BOOST_CHECK(transformed_y[i] == expected_y[i]); - BOOST_CHECK(transformed_z[i] == expected_z[i]); + BOOST_CHECK_SMALL(abs(transformed_x[i] - expected_x[i]), eps); + BOOST_CHECK_SMALL(abs(transformed_y[i] - expected_y[i]), eps); + BOOST_CHECK_SMALL(abs(transformed_z[i] - expected_z[i]), eps); + } +} + +BOOST_AUTO_TEST_CASE(cylinder_to_cartesian_with_axis_with_phi_test) { + constexpr auto eps = 1e-14; + // tilted orthogonal basis + auto const x = + (Vector3d{{1, 0, 0}} - (1. / 3.) * Vector3d{{1, 1, 1}}).normalize(); + auto const y = (Vector3d{{0, 1, -1}}).normalize(); + auto const z = (Vector3d{{1, 1, 1}}).normalize(); + // check simple transformation without orientation + { + auto const x_cyl = transform_coordinate_cartesian_to_cylinder(x, z); + auto const y_cyl = transform_coordinate_cartesian_to_cylinder(y, z); + auto const z_cyl = transform_coordinate_cartesian_to_cylinder(z, z); + auto const x_cart = transform_coordinate_cylinder_to_cartesian(x_cyl, z); + auto const y_cart = transform_coordinate_cylinder_to_cartesian(y_cyl, z); + auto const z_cart = transform_coordinate_cylinder_to_cartesian(z_cyl, z); + for (int i = 0; i < 3; ++i) { + BOOST_CHECK_SMALL(abs(x_cart[i] - x[i]), eps); + BOOST_CHECK_SMALL(abs(y_cart[i] - y[i]), eps); + BOOST_CHECK_SMALL(abs(z_cart[i] - z[i]), eps); + } + } + // check transformation with orientation + { + auto const x_cyl = transform_coordinate_cartesian_to_cylinder(x, z, y); + auto const y_cyl = transform_coordinate_cartesian_to_cylinder(y, z, y); + auto const z_cyl = transform_coordinate_cartesian_to_cylinder(z, z, y); + auto const x_cart = transform_coordinate_cylinder_to_cartesian(x_cyl, z, y); + auto const y_cart = transform_coordinate_cylinder_to_cartesian(y_cyl, z, y); + auto const z_cart = transform_coordinate_cylinder_to_cartesian(z_cyl, z, y); + for (int i = 0; i < 3; ++i) { + BOOST_CHECK_SMALL(abs(x_cart[i] - x[i]), eps); + BOOST_CHECK_SMALL(abs(y_cart[i] - y[i]), eps); + BOOST_CHECK_SMALL(abs(z_cart[i] - z[i]), eps); + } + } + // check transformation with orientation for another angle + { + auto const u = vec_rotate(z, Utils::pi() / 3.0, x); + auto const v = vec_rotate(z, Utils::pi() / 3.0, y); + auto const u_cyl = transform_coordinate_cartesian_to_cylinder(u, z, y); + auto const v_cyl = transform_coordinate_cartesian_to_cylinder(v, z, y); + auto const u_cart = transform_coordinate_cylinder_to_cartesian(u_cyl, z, y); + auto const v_cart = transform_coordinate_cylinder_to_cartesian(v_cyl, z, y); + for (int i = 0; i < 3; ++i) { + BOOST_CHECK_SMALL(abs(u_cart[i] - u[i]), eps); + BOOST_CHECK_SMALL(abs(v_cart[i] - v[i]), eps); + } } } diff --git a/src/utils/tests/interval_test.cpp b/src/utils/tests/interval_test.cpp new file mode 100644 index 00000000000..8f619d3cd10 --- /dev/null +++ b/src/utils/tests/interval_test.cpp @@ -0,0 +1,38 @@ +/* + * Copyright (C) 2020 The ESPResSo project + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#define BOOST_TEST_MODULE interval test +#define BOOST_TEST_DYN_LINK +#include + +#include +#include + +#include + +BOOST_AUTO_TEST_CASE(interval_test) { + constexpr auto eps = 1e-14; + constexpr auto pi = Utils::pi(); + BOOST_CHECK_EQUAL(Utils::interval(-pi, -pi, pi), -pi); + BOOST_CHECK_EQUAL(Utils::interval(pi, -pi, pi), pi); + for (int i = -24; i < 24; ++i) { + double const theta = i * pi / 6; + double const wrapped = Utils::interval(theta, -pi, pi); + double const ref_val = std::atan2(std::sin(theta), std::cos(theta)); + BOOST_CHECK_SMALL(std::abs(std::fmod(wrapped - ref_val, 2 * pi)), eps); + } +} diff --git a/src/utils/tests/vec_rotate_test.cpp b/src/utils/tests/vec_rotate_test.cpp index 18b9b14ef45..f1d5727a2cd 100644 --- a/src/utils/tests/vec_rotate_test.cpp +++ b/src/utils/tests/vec_rotate_test.cpp @@ -22,11 +22,9 @@ #include #include #include -using Utils::vec_rotate; #include #include -#include BOOST_AUTO_TEST_CASE(rotation) { using std::cos; @@ -43,20 +41,16 @@ BOOST_AUTO_TEST_CASE(rotation) { auto const expected = cos(t) * v + sin(t) * vector_product(k, v) + (1. - cos(t)) * (k * v) * k; - auto const is = vec_rotate(k, t, v); + auto const is = Utils::vec_rotate(k, t, v); auto const rel_diff = (expected - is).norm() / expected.norm(); BOOST_CHECK(rel_diff < std::numeric_limits::epsilon()); } -BOOST_AUTO_TEST_CASE(rotation_params) { - Utils::Vector3d v1 = {1.0, 0.0, 0.0}; - Utils::Vector3d v2 = {1.0, 1.0, 0.0}; +BOOST_AUTO_TEST_CASE(angle_between) { + Utils::Vector3d const v1 = {1.0, 0.0, 0.0}; + Utils::Vector3d const v2 = {1.0, 1.0, 0.0}; - double angle; - Utils::Vector3d rotation_axis; - std::tie(angle, rotation_axis) = Utils::rotation_params(v1, v2); + auto const angle = Utils::angle_between(v1, v2); BOOST_CHECK_CLOSE(angle, Utils::pi() / 4.0, 1e-7); - BOOST_CHECK_SMALL((rotation_axis * v1), 1e-7); - BOOST_CHECK_SMALL((rotation_axis * v2), 1e-7); }