Utilities/detail/periodic.hpp

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