EIC code displayed by LXR |
|
|
|
|||
File indexing completed on 2025-09-18 09:09:41
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 Euclidian (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 // Calculate the Euclidian (2) distance between two points 0054 template<class T, size_type N> 0055 [[nodiscard]] inline CELER_FUNCTION T distance(Array<T, N> const& x, 0056 Array<T, N> const& y); 0057 0058 //---------------------------------------------------------------------------// 0059 // Calculate a cartesian unit vector from spherical coordinates 0060 template<class T> 0061 [[nodiscard]] inline CELER_FUNCTION Array<T, 3> 0062 from_spherical(T costheta, T phi); 0063 0064 //---------------------------------------------------------------------------// 0065 // Rotate the direction 'dir' according to the reference rotation axis 'rot' 0066 template<class T> 0067 [[nodiscard]] inline CELER_FUNCTION Array<T, 3> 0068 rotate(Array<T, 3> const& dir, Array<T, 3> const& rot); 0069 0070 //---------------------------------------------------------------------------// 0071 // INLINE DEFINITIONS 0072 //---------------------------------------------------------------------------// 0073 /*! 0074 * Increment a vector by another vector multiplied by a scalar. 0075 * 0076 * \f[ 0077 * y \gets \alpha x + y 0078 * \f] 0079 * 0080 * Note that this uses \c celeritas::fma which supports types other than 0081 * floating point. 0082 */ 0083 template<class T, size_type N> 0084 CELER_FUNCTION void axpy(T a, Array<T, N> const& x, Array<T, N>* y) 0085 { 0086 CELER_EXPECT(y); 0087 for (size_type i = 0; i != N; ++i) 0088 { 0089 (*y)[i] = fma(a, x[i], (*y)[i]); 0090 } 0091 } 0092 0093 //---------------------------------------------------------------------------// 0094 /*! 0095 * Dot product of two vectors. 0096 * 0097 * Note that this uses \c celeritas::fma which supports types other than 0098 * floating point. 0099 */ 0100 template<class T, size_type N> 0101 CELER_FUNCTION T dot_product(Array<T, N> const& x, Array<T, N> const& y) 0102 { 0103 T result{}; 0104 for (size_type i = 0; i != N; ++i) 0105 { 0106 result = fma(x[i], y[i], result); 0107 } 0108 return result; 0109 } 0110 0111 //---------------------------------------------------------------------------// 0112 /*! 0113 * Cross product of two space vectors. 0114 */ 0115 template<class T> 0116 CELER_FUNCTION Array<T, 3> 0117 cross_product(Array<T, 3> const& x, Array<T, 3> const& y) 0118 { 0119 return {x[1] * y[2] - x[2] * y[1], 0120 x[2] * y[0] - x[0] * y[2], 0121 x[0] * y[1] - x[1] * y[0]}; 0122 } 0123 0124 //---------------------------------------------------------------------------// 0125 /*! 0126 * Calculate the Euclidian (2) norm of a vector. 0127 */ 0128 template<class T, size_type N> 0129 CELER_FUNCTION T norm(Array<T, N> const& v) 0130 { 0131 return std::sqrt(dot_product(v, v)); 0132 } 0133 0134 //---------------------------------------------------------------------------// 0135 /*! 0136 * Construct a unit vector. 0137 * 0138 * Unit vectors have an Euclidian norm magnitude of 1. 0139 */ 0140 template<class T, size_type N> 0141 CELER_FUNCTION Array<T, N> make_unit_vector(Array<T, N> const& v) 0142 { 0143 Array<T, N> result{v}; 0144 T const scale_factor = 1 / norm(result); 0145 for (auto& el : result) 0146 { 0147 el *= scale_factor; 0148 } 0149 return result; 0150 } 0151 0152 //---------------------------------------------------------------------------// 0153 /*! 0154 * Calculate the Euclidian (2) distance between two points. 0155 */ 0156 template<class T, size_type N> 0157 CELER_FUNCTION T distance(Array<T, N> const& x, Array<T, N> const& y) 0158 { 0159 T dist_sq = 0; 0160 for (size_type i = 0; i != N; ++i) 0161 { 0162 dist_sq += ipow<2>(y[i] - x[i]); 0163 } 0164 return std::sqrt(dist_sq); 0165 } 0166 0167 //---------------------------------------------------------------------------// 0168 /*! 0169 * Calculate a Cartesian vector from spherical coordinates. 0170 * 0171 * Theta is the angle between the \em z axis and the outgoing vector, and \c 0172 * phi is the angle between the \em x axis and the projection of the vector 0173 * onto the \em x-y plane. 0174 */ 0175 template<class T> 0176 inline CELER_FUNCTION Array<T, 3> from_spherical(T costheta, T phi) 0177 { 0178 CELER_EXPECT(costheta >= -1 && costheta <= 1); 0179 0180 T const sintheta = std::sqrt(1 - costheta * costheta); 0181 return {sintheta * std::cos(phi), sintheta * std::sin(phi), costheta}; 0182 } 0183 0184 //---------------------------------------------------------------------------// 0185 /*! 0186 * Rotate a direction about the given scattering direction. 0187 * 0188 * The equivalent to calling the Shift transport code's \code 0189 void cartesian_vector_transform( 0190 double costheta, 0191 double phi, 0192 Vector_View vector); 0193 * \endcode 0194 * is the call 0195 * \code 0196 vector = rotate(from_spherical(costheta, phi), vector); 0197 * \endcode 0198 * 0199 * This code effectively decomposes the given rotation vector \c rot into two 0200 * sequential transform matrices, one with an angle \em theta about the \em y 0201 * axis and one about \em phi rotating around the \em z axis. These two angles 0202 * are the spherical coordinate transform of the given \c rot cartesian 0203 * direction vector. 0204 * 0205 * There is some extra code in here to deal with loss of precision when the 0206 * incident direction is along the \em z axis. As \c rot approaches \em z, the 0207 * azimuthal angle \f$ \phi \f$ must be calculated carefully from both the 0208 * \em x and \em y components of the vector, not independently. 0209 * If \c rot actually equals \em z 0210 * then the azimuthal angle is completely indeterminate so we arbitrarily 0211 * choose \f$ \phi = 0 \f$. 0212 * 0213 * This function is often used for calculating exiting scattering angles. In 0214 * that case, \c dir is the exiting angle from the scattering calculation, and 0215 * \c rot is the original direction of the particle. The direction vectors are 0216 * defined as 0217 * \f[ 0218 \vec{\Omega} 0219 = \sin\theta\cos\phi\vec{i} 0220 + \sin\theta\sin\phi\vec{j} 0221 + \cos\theta\vec{k} \,. 0222 * \f] 0223 */ 0224 template<class T> 0225 inline CELER_FUNCTION Array<T, 3> 0226 rotate(Array<T, 3> const& dir, Array<T, 3> const& rot) 0227 { 0228 CELER_EXPECT(is_soft_unit_vector(dir)); 0229 CELER_EXPECT(is_soft_unit_vector(rot)); 0230 0231 // Direction enumeration 0232 enum 0233 { 0234 X = 0, 0235 Y = 1, 0236 Z = 2 0237 }; 0238 0239 // Transform direction vector into theta, phi so we can use it as a 0240 // rotation matrix 0241 T sintheta = std::sqrt(1 - ipow<2>(rot[Z])); 0242 T cosphi; 0243 T sinphi; 0244 0245 if (sintheta >= detail::RealVecTraits<T>::min_accurate_sintheta()) 0246 { 0247 // Typical case: far enough from z axis to assume the X and Y 0248 // components have a hypotenuse of 1 within epsilon tolerance 0249 T const inv_sintheta = 1 / (sintheta); 0250 cosphi = rot[X] * inv_sintheta; 0251 sinphi = rot[Y] * inv_sintheta; 0252 } 0253 else if (sintheta > 0) 0254 { 0255 // Avoid catastrophic roundoff error by normalizing x/y components 0256 cosphi = rot[X] / hypot(rot[X], rot[Y]); 0257 sinphi = std::sqrt(1 - ipow<2>(cosphi)); 0258 } 0259 else 0260 { 0261 // NaN or 0: choose an arbitrary azimuthal angle for the incident dir 0262 cosphi = 1; 0263 sinphi = 0; 0264 } 0265 0266 Array<T, 3> result 0267 = {(rot[Z] * dir[X] + sintheta * dir[Z]) * cosphi - sinphi * dir[Y], 0268 (rot[Z] * dir[X] + sintheta * dir[Z]) * sinphi + cosphi * dir[Y], 0269 -sintheta * dir[X] + rot[Z] * dir[Z]}; 0270 0271 // Always normalize to prevent roundoff error from propagating 0272 return make_unit_vector(result); 0273 } 0274 0275 //---------------------------------------------------------------------------// 0276 } // 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 |
|