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