ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
coordinate_transformation.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20#pragma once
21
22/**
23 * @file
24 * Convert coordinates from the Cartesian system to the cylindrical system.
25 * The transformation functions are provided with three overloads:
26 * - one function for the trivial Cartesian <-> cylindrical transformation
27 * - one function to transform from/to a cylindrical system with custom axis
28 * (extra @p axis argument, keep in mind the angle phi is under-defined)
29 * - one function to transform from/to an oriented cylindrical system with
30 * custom axis (extra @p orientation argument, the angle phi is well-defined)
31 */
32
33#include "utils/Vector.hpp"
35#include "utils/matrix.hpp"
36
37#include <cassert>
38#include <cmath>
39
40namespace Utils {
41
42/**
43 * @brief Basis change.
44 */
45inline Vector3d basis_change(Vector3d const &b1, Vector3d const &b2,
46 Vector3d const &b3, Vector3d const &v,
47 bool reverse = false) {
48 auto const e_x = b1.normalized();
49 auto const e_y = b2.normalized();
50 auto const e_z = b3.normalized();
51 auto const M = Matrix<double, 3, 3>{
52 {e_x[0], e_x[1], e_x[2]},
53 {e_y[0], e_y[1], e_y[2]},
54 {e_z[0], e_z[1],
55 e_z[2]}}.transposed();
56 if (reverse) {
57 return M * v;
58 }
59 return M.inversed() * v;
60}
61
62/**
63 * @brief Coordinate transformation from Cartesian to cylindrical coordinates.
64 * The origins and z-axis of the coordinate systems co-incide.
65 * The @f$ \phi = 0 @f$ direction corresponds to the x-axis in the
66 * original coordinate system.
67 * @param pos Vector to transform
68 */
69inline Vector3d
71 auto const r = std::sqrt(pos[0] * pos[0] + pos[1] * pos[1]);
72 auto const phi = std::atan2(pos[1], pos[0]);
73 return {r, phi, pos[2]};
74}
75
76/**
77 * @brief Coordinate transformation from Cartesian to cylindrical coordinates
78 * with change of basis. The origins of the coordinate systems co-incide.
79 *
80 * If the parameter @p axis is not equal to <tt>[0, 0, 1]</tt>, the value
81 * of the angle @f$ \phi @f$ in cylindrical coordinates is under-defined.
82 * To fully define it, it is necessary to provide an orientation vector
83 * in Cartesian coordinates that will be used as the reference point
84 * (i.e. such that @f$ \phi = 0 @f$), by default it is the x-axis.
85 *
86 * @param pos Vector to transform
87 * @param axis Longitudinal axis of the cylindrical coordinates
88 * @param orientation Reference point (in untransformed coordinates) for
89 * which @f$ \phi = 0 @f$
90 */
92 Vector3d const &pos, Vector3d const &axis, Vector3d const &orientation) {
93 // check that axis and orientation are orthogonal
94 assert(std::abs(axis * orientation) <
95 5 * std::numeric_limits<double>::epsilon());
96 auto const rotation_axis = vector_product(axis, orientation);
97 auto const pos_t = basis_change(orientation, rotation_axis, axis, pos);
99}
100
101/**
102 * @brief Coordinate transformation from cylindrical to Cartesian coordinates.
103 * The origins and z-axis of the coordinate systems co-incide.
104 * The @f$ \phi = 0 @f$ direction corresponds to the x-axis in the
105 * transformed coordinate system.
106 * @param pos Vector to transform
107 */
108inline Vector3d
110 auto const &rho = pos[0];
111 auto const &phi = pos[1];
112 auto const &z = pos[2];
113 return {rho * std::cos(phi), rho * std::sin(phi), z};
114}
115
116/**
117 * @brief Coordinate transformation from cylindrical to Cartesian coordinates
118 * with change of basis. The origins of the coordinate systems co-incide.
119 *
120 * If the parameter @p axis is not equal to <tt>[0, 0, 1]</tt>, the value
121 * of the angle @f$ \phi @f$ in cylindrical coordinates is under-defined.
122 * To fully define it, it is necessary to provide an orientation vector
123 * in Cartesian coordinates that will be used as the reference point
124 * (i.e. such that @f$ \phi = 0 @f$).
125 *
126 * @param pos Vector to transform
127 * @param axis Longitudinal axis of the cylindrical coordinates
128 * @param orientation Reference point (in Cartesian coordinates) for
129 * which @f$ \phi = 0 @f$
130 */
132 Vector3d const &pos, Vector3d const &axis, Vector3d const &orientation) {
133 // check that axis and orientation are orthogonal
134 assert(std::abs(axis * orientation) <
135 5 * std::numeric_limits<double>::epsilon());
136 auto const rotation_axis = vector_product(axis, orientation);
137 auto const pos_t = transform_coordinate_cylinder_to_cartesian(pos);
138 return basis_change(orientation, rotation_axis, axis, pos_t, true);
139}
140
141/**
142 * @brief Vector transformation from Cartesian to cylindrical coordinates.
143 * @param vec Vector to transform
144 * @param axis Longitudinal axis of the cylindrical coordinates
145 * @param pos Origin of the vector
146 */
148 Vector3d const &axis,
149 Vector3d const &pos) {
150 static auto const z_axis = Vector3d{{0, 0, 1}};
151 auto const angle = angle_between(axis, z_axis);
152 auto const rotation_axis = Utils::vector_product(axis, z_axis).normalize();
153 auto const rotated_pos = vec_rotate(rotation_axis, angle, pos);
154 auto const rotated_vec = vec_rotate(rotation_axis, angle, vec);
155 auto const r = std::sqrt(rotated_pos[0] * rotated_pos[0] +
156 rotated_pos[1] * rotated_pos[1]);
157 // v_r = (x * v_x + y * v_y) / sqrt(x^2 + y^2)
158 auto const v_r =
159 (rotated_pos[0] * rotated_vec[0] + rotated_pos[1] * rotated_vec[1]) / r;
160 // v_phi = (x * v_y - y * v_x ) / sqrt(x^2 + y^2)
161 auto const v_phi =
162 (rotated_pos[0] * rotated_vec[1] - rotated_pos[1] * rotated_vec[0]) / r;
163 return Vector3d{v_r, v_phi, rotated_vec[2]};
164}
165
166} // namespace Utils
Vector implementation and trait types for boost qvm interoperability.
Vector normalized() const
Definition Vector.hpp:157
Matrix implementation and trait types for boost qvm interoperability.
VectorXd< 3 > Vector3d
Definition Vector.hpp:164
Vector< T, 3 > vector_product(Vector< T, 3 > const &a, Vector< T, 3 > const &b)
Definition Vector.hpp:368
Vector3d vec_rotate(const Vector3d &axis, double angle, const Vector3d &vector)
Rotate a vector around an axis.
double angle_between(Vector3d const &v1, Vector3d const &v2)
Determine the angle between two vectors.
Vector3d transform_vector_cartesian_to_cylinder(Vector3d const &vec, Vector3d const &axis, Vector3d const &pos)
Vector transformation from Cartesian to cylindrical coordinates.
Vector3d transform_coordinate_cylinder_to_cartesian(Vector3d const &pos)
Coordinate transformation from cylindrical to Cartesian coordinates.
Vector3d transform_coordinate_cartesian_to_cylinder(Vector3d const &pos)
Coordinate transformation from Cartesian to cylindrical coordinates.
Vector3d basis_change(Vector3d const &b1, Vector3d const &b2, Vector3d const &b3, Vector3d const &v, bool reverse=false)
Basis change.
Matrix representation with static size.
Definition matrix.hpp:65
Matrix< T, Cols, Rows > transposed() const
Retrieve a transposed copy of the matrix.
Definition matrix.hpp:189