EIC code displayed by LXR |
|
|
|
|||
File indexing completed on 2026-04-17 08:20:20
0001 //------------------------------- -*- C++ -*- -------------------------------// 0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details 0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT) 0004 //---------------------------------------------------------------------------// 0005 //! \file corecel/math/ArrayUtils.hh 0006 //! \brief Math functions using celeritas::Array 0007 //---------------------------------------------------------------------------// 0008 #pragma once 0009 0010 #include <cmath> 0011 0012 #include "corecel/Assert.hh" 0013 #include "corecel/Types.hh" 0014 #include "corecel/cont/Array.hh" 0015 0016 #include "Algorithms.hh" 0017 #include "ArraySoftUnit.hh" 0018 #include "SoftEqual.hh" 0019 0020 #include "detail/ArrayUtilsImpl.hh" 0021 0022 namespace celeritas 0023 { 0024 //---------------------------------------------------------------------------// 0025 // Perform y <- ax + y 0026 template<class T, size_type N> 0027 inline CELER_FUNCTION void axpy(T a, Array<T, N> const& x, Array<T, N>* y); 0028 0029 //---------------------------------------------------------------------------// 0030 // Calculate product of two vectors 0031 template<class T, size_type N> 0032 [[nodiscard]] inline CELER_FUNCTION T dot_product(Array<T, N> const& x, 0033 Array<T, N> const& y); 0034 0035 //---------------------------------------------------------------------------// 0036 // Calculate product of two vectors 0037 template<class T> 0038 [[nodiscard]] inline CELER_FUNCTION Array<T, 3> 0039 cross_product(Array<T, 3> const& x, Array<T, 3> const& y); 0040 0041 //---------------------------------------------------------------------------// 0042 // Calculate the Euclidean (2) norm of a vector 0043 template<class T, size_type N> 0044 [[nodiscard]] inline CELER_FUNCTION T norm(Array<T, N> const& vec); 0045 0046 //---------------------------------------------------------------------------// 0047 // Construct a vector with unit magnitude 0048 template<class T, size_type N> 0049 [[nodiscard]] inline CELER_FUNCTION Array<T, N> 0050 make_unit_vector(Array<T, N> const& v); 0051 0052 //---------------------------------------------------------------------------// 0053 // Return x - (x . y) * y for a unit vector y 0054 template<class T, size_type N> 0055 [[nodiscard]] inline CELER_FUNCTION Array<T, N> 0056 make_orthogonal(Array<T, N> const& x, Array<T, N> const& y); 0057 0058 //---------------------------------------------------------------------------// 0059 // Calculate the Euclidean (2) distance between two points 0060 template<class T, size_type N> 0061 [[nodiscard]] inline CELER_FUNCTION T distance(Array<T, N> const& x, 0062 Array<T, N> const& y); 0063 0064 //---------------------------------------------------------------------------// 0065 // Calculate a cartesian unit vector from spherical coordinates 0066 template<class T> 0067 [[nodiscard]] inline CELER_FUNCTION Array<T, 3> 0068 from_spherical(T costheta, T phi); 0069 0070 //---------------------------------------------------------------------------// 0071 // Rotate the direction 'dir' according to the reference rotation axis 'rot' 0072 template<class T> 0073 [[nodiscard]] inline CELER_FUNCTION Array<T, 3> 0074 rotate(Array<T, 3> const& dir, Array<T, 3> const& rot); 0075 0076 //---------------------------------------------------------------------------// 0077 // INLINE DEFINITIONS 0078 //---------------------------------------------------------------------------// 0079 /*! 0080 * Increment a vector by another vector multiplied by a scalar. 0081 * 0082 * \f[ 0083 * y \gets \alpha x + y 0084 * \f] 0085 * 0086 * Note that this uses \c celeritas::fma which supports types other than 0087 * floating point. 0088 */ 0089 template<class T, size_type N> 0090 CELER_FUNCTION void axpy(T a, Array<T, N> const& x, Array<T, N>* y) 0091 { 0092 CELER_EXPECT(y); 0093 for (size_type i = 0; i != N; ++i) 0094 { 0095 (*y)[i] = fma(a, x[i], (*y)[i]); 0096 } 0097 } 0098 0099 //---------------------------------------------------------------------------// 0100 /*! 0101 * Dot product of two vectors. 0102 * 0103 * Note that this uses \c celeritas::fma which supports types other than 0104 * floating point. 0105 */ 0106 template<class T, size_type N> 0107 CELER_FUNCTION T dot_product(Array<T, N> const& x, Array<T, N> const& y) 0108 { 0109 T result{}; 0110 for (size_type i = 0; i != N; ++i) 0111 { 0112 result = fma(x[i], y[i], result); 0113 } 0114 return result; 0115 } 0116 0117 //---------------------------------------------------------------------------// 0118 /*! 0119 * Cross product of two space vectors. 0120 */ 0121 template<class T> 0122 CELER_FUNCTION Array<T, 3> 0123 cross_product(Array<T, 3> const& x, Array<T, 3> const& y) 0124 { 0125 return {x[1] * y[2] - x[2] * y[1], 0126 x[2] * y[0] - x[0] * y[2], 0127 x[0] * y[1] - x[1] * y[0]}; 0128 } 0129 0130 //---------------------------------------------------------------------------// 0131 /*! 0132 * Calculate the Euclidean (2) norm of a vector. 0133 */ 0134 template<class T, size_type N> 0135 CELER_FUNCTION T norm(Array<T, N> const& v) 0136 { 0137 return std::sqrt(dot_product(v, v)); 0138 } 0139 0140 //---------------------------------------------------------------------------// 0141 /*! 0142 * Construct a unit vector. 0143 * 0144 * Unit vectors have an Euclidean norm magnitude of 1. 0145 */ 0146 template<class T, size_type N> 0147 CELER_FUNCTION Array<T, N> make_unit_vector(Array<T, N> const& v) 0148 { 0149 Array<T, N> result{v}; 0150 T const scale_factor = 1 / norm(result); 0151 for (auto& el : result) 0152 { 0153 el *= scale_factor; 0154 } 0155 return result; 0156 } 0157 0158 //---------------------------------------------------------------------------// 0159 /*! 0160 * Return the component of \em x that is orthogonal to the unit vector \em y. 0161 * 0162 * In this implementation, \em y must be normalized, and the result is not 0163 * normalized. 0164 * 0165 * \f[ 0166 \mathbf{x}' \gets \mathbf{x} - (\mathbf{x} \cdot \mathbf{y}) \mathbf{y} 0167 \, , \quad \|\mathbf{y}\| = 1 0168 \f] 0169 */ 0170 template<class T, size_type N> 0171 [[nodiscard]] inline CELER_FUNCTION Array<T, N> 0172 make_orthogonal(Array<T, N> const& x, Array<T, N> const& y) 0173 { 0174 CELER_EXPECT(is_soft_unit_vector(y)); 0175 Array<T, N> result{x}; 0176 axpy(-dot_product(x, y), y, &result); 0177 return result; 0178 } 0179 0180 //---------------------------------------------------------------------------// 0181 /*! 0182 * Calculate the Euclidean (2) distance between two points. 0183 */ 0184 template<class T, size_type N> 0185 CELER_FUNCTION T distance(Array<T, N> const& x, Array<T, N> const& y) 0186 { 0187 T dist_sq = 0; 0188 for (size_type i = 0; i != N; ++i) 0189 { 0190 dist_sq += ipow<2>(y[i] - x[i]); 0191 } 0192 return std::sqrt(dist_sq); 0193 } 0194 0195 //---------------------------------------------------------------------------// 0196 /*! 0197 * Calculate a Cartesian vector from spherical coordinates. 0198 * 0199 * Theta is the angle between the \em z axis and the outgoing vector, and \c 0200 * phi is the angle between the \em x axis and the projection of the vector 0201 * onto the \em x-y plane. 0202 */ 0203 template<class T> 0204 inline CELER_FUNCTION Array<T, 3> from_spherical(T costheta, T phi) 0205 { 0206 CELER_EXPECT(costheta >= -1 && costheta <= 1); 0207 0208 T const sintheta = std::sqrt(1 - costheta * costheta); 0209 return {sintheta * std::cos(phi), sintheta * std::sin(phi), costheta}; 0210 } 0211 0212 //---------------------------------------------------------------------------// 0213 /*! 0214 * Rotate a direction about the given scattering direction. 0215 * 0216 * The equivalent to calling the Shift transport code's \code 0217 void cartesian_vector_transform( 0218 double costheta, 0219 double phi, 0220 Vector_View vector); 0221 * \endcode 0222 * is the call 0223 * \code 0224 vector = rotate(from_spherical(costheta, phi), vector); 0225 * \endcode 0226 * 0227 * This code effectively decomposes the given rotation vector \c rot into two 0228 * sequential transform matrices, one with an angle \em theta about the \em y 0229 * axis and one about \em phi rotating around the \em z axis. These two angles 0230 * are the spherical coordinate transform of the given \c rot cartesian 0231 * direction vector. 0232 * 0233 * There is some extra code in here to deal with loss of precision when the 0234 * incident direction is along the \em z axis. As \c rot approaches \em z, the 0235 * azimuthal angle \f$ \phi \f$ must be calculated carefully from both the 0236 * \em x and \em y components of the vector, not independently. 0237 * If \c rot actually equals \em z 0238 * then the azimuthal angle is completely indeterminate so we arbitrarily 0239 * choose \f$ \phi = 0 \f$. 0240 * 0241 * This function is often used for calculating exiting scattering angles. In 0242 * that case, \c dir is the exiting angle from the scattering calculation, and 0243 * \c rot is the original direction of the particle. The direction vectors are 0244 * defined as 0245 * \f[ 0246 \vec{\Omega} 0247 = \sin\theta\cos\phi\vec{i} 0248 + \sin\theta\sin\phi\vec{j} 0249 + \cos\theta\vec{k} \,. 0250 * \f] 0251 */ 0252 template<class T> 0253 inline CELER_FUNCTION Array<T, 3> 0254 rotate(Array<T, 3> const& dir, Array<T, 3> const& rot) 0255 { 0256 CELER_EXPECT(is_soft_unit_vector(dir)); 0257 CELER_EXPECT(is_soft_unit_vector(rot)); 0258 0259 // Direction enumeration 0260 enum 0261 { 0262 X = 0, 0263 Y = 1, 0264 Z = 2 0265 }; 0266 0267 // Transform direction vector into theta, phi so we can use it as a 0268 // rotation matrix 0269 T sintheta = std::sqrt(1 - ipow<2>(rot[Z])); 0270 T cosphi; 0271 T sinphi; 0272 0273 if (sintheta >= detail::RealVecTraits<T>::min_accurate_sintheta()) 0274 { 0275 // Typical case: far enough from z axis to assume the X and Y 0276 // components have a hypotenuse of 1 within epsilon tolerance 0277 T const inv_sintheta = 1 / (sintheta); 0278 cosphi = rot[X] * inv_sintheta; 0279 sinphi = rot[Y] * inv_sintheta; 0280 } 0281 else if (sintheta > 0) 0282 { 0283 // Avoid catastrophic roundoff error by normalizing x/y components 0284 cosphi = rot[X] / hypot(rot[X], rot[Y]); 0285 sinphi = std::sqrt(1 - ipow<2>(cosphi)); 0286 } 0287 else 0288 { 0289 // NaN or 0: choose an arbitrary azimuthal angle for the incident dir 0290 cosphi = 1; 0291 sinphi = 0; 0292 } 0293 0294 Array<T, 3> result 0295 = {(rot[Z] * dir[X] + sintheta * dir[Z]) * cosphi - sinphi * dir[Y], 0296 (rot[Z] * dir[X] + sintheta * dir[Z]) * sinphi + cosphi * dir[Y], 0297 -sintheta * dir[X] + rot[Z] * dir[Z]}; 0298 0299 // Always normalize to prevent roundoff error from propagating 0300 return make_unit_vector(result); 0301 } 0302 0303 //---------------------------------------------------------------------------// 0304 } // namespace celeritas
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |
|