Utilities/detail/periodic.hpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include <cmath>
0012 #include <numbers>
0013
0014 namespace Acts::detail {
0015
0016 /// Wrap a periodic value back into the nominal range.
0017 template <typename T>
0018 inline T wrap_periodic(T value, T start, T range) {
0019   using std::floor;
0020   // only wrap if really necessary
0021   T diff = value - start;
0022   return ((0 <= diff) && (diff < range))
0023              ? value
0024              : (value - range * floor(diff / range));
0025 }
0026
0027 /// Compute the minimal `lhs - rhs` using the periodicity.
0028 ///
0029 /// Imagine you have two values within the nominal range: `l` is close to the
0030 /// lower edge and `u` is close to the upper edge. The naive difference between
0031 /// the two is almost as large as the range itself. If we move `l` to its
0032 /// equivalent value outside the nominal range, i.e. just above the upper edge,
0033 /// the effective absolute difference becomes smaller.
0034 ///
0035 /// @note The sign of the returned value can be different from `lhs - rhs`
0036 template <typename T>
0037 inline T difference_periodic(T lhs, T rhs, T range) {
0038   using std::fmod;
0039   T delta = fmod(lhs - rhs, range);
0040   // check if |delta| is larger than half the range. if that is the case, we
0041   // can move either rhs/lhs by one range/period to get a smaller |delta|.
0042   if ((2 * delta) < -range) {
0043     delta += range;
0044   } else if (range <= (2 * delta)) {
0045     delta -= range;
0046   }
0047   return delta;
0048 }
0049
0050 /// Calculate the equivalent angle in the [0, 2*pi) range.
0051 template <typename T>
0052 inline T radian_pos(T x) {
0053   return wrap_periodic<T>(x, T{0}, T{2 * std::numbers::pi});
0054 }
0055
0056 /// Calculate the equivalent angle in the [-pi, pi) range.
0057 template <typename T>
0058 inline T radian_sym(T x) {
0059   return wrap_periodic<T>(x, -std::numbers::pi_v<T>, T{2 * std::numbers::pi});
0060 }
0061
0062 /// Ensure both phi and theta direction angles are within the allowed range.
0063 ///
0064 /// @param[in] phi Transverse direction angle
0065 /// @param[in] theta Longitudinal direction angle
0066 /// @return pair<phi,theta> containing the updated angles
0067 ///
0068 /// The phi angle is truly cyclic, i.e. all values outside the nominal range
0069 /// [-pi,pi) have a corresponding value inside nominal range, independent from
0070 /// the theta angle. The theta angle is more complicated. Imagine that the two
0071 /// angles describe a position on the unit sphere. If theta moves outside its
0072 /// nominal range [0,pi], we are moving over one of the two poles of the unit
0073 /// sphere along the great circle defined by phi. The angles still describe a
0074 /// valid position on the unit sphere, but to describe it with angles within
0075 /// their nominal range, both phi and theta need to be updated; when moving over
0076 /// the poles, phi needs to be flipped by 180degree to allow theta to remain
0077 /// within its nominal range.
0078 template <typename T>
0079 inline std::pair<T, T> normalizePhiTheta(T phi, T theta) {
0080   // wrap to [0,2pi). while the nominal range of theta is [0,pi], it is
0081   // periodic, i.e. describes identical positions, in the full [0,2pi) range.
0082   // moving it first to the periodic range simplifies further steps as the
0083   // possible range of theta becomes fixed.
0084   theta = radian_pos(theta);
0085   if (std::numbers::pi_v<T> < theta) {
0086     // theta is in the second half of the great circle and outside its nominal
0087     // range. need to change both phi and theta to be within range.
0088     phi += std::numbers::pi_v<T>;
0089     theta = T{2 * std::numbers::pi} - theta;
0090   }
0091   return {radian_sym(phi), theta};
0092 }
0093
0094 }  // namespace Acts::detail