VecGeom/volumes/EllipticUtilities.h

0001 // This file is part of VecGeom and is distributed under the
0002 // conditions in the file LICENSE.txt in the top directory.
0003 // For the full list of authors see CONTRIBUTORS.txt and `git log`.
0004
0005 /// Collection of auxiliary utilities for elliptical shapes
0006 /// @file volumes/EllipticUtilities.h
0007 /// @author Original implementation by Evgueni Tcherniaev (evc)
0008
0009 #ifndef VECGEOM_ELLIPTICUTILITIES_H_
0010 #define VECGEOM_ELLIPTICUTILITIES_H_
0011
0012 #include "VecGeom/base/Global.h"
0013 #include "VecGeom/base/RNG.h"
0014 #include "VecGeom/base/Vector2D.h"
0015 #include "VecGeom/base/Vector3D.h"
0016
0017 namespace vecgeom {
0018
0019 inline namespace VECGEOM_IMPL_NAMESPACE {
0020
0021 namespace EllipticUtilities {
0022
0023 /// Computes the Complete Elliptical Integral of the Second Kind
0024 /// @param e eccentricity
0025 //
0026 VECCORE_ATT_HOST_DEVICE
0027 VECGEOM_FORCE_INLINE
0028 Precision comp_ellint_2(Precision e)
0029 {
0030   // The algorithm is based upon Carlson B.C., "Computation of real
0031   // or complex elliptic integrals", Numerical Algorithms,
0032   // Volume 10, Issue 1, 1995 (see equations 2.36 - 2.39)
0033   //
0034   const Precision eps = 1. / 134217728.; // 1 / 2^27
0035
0036   Precision a = 1.;
0037   Precision b = vecCore::math::Sqrt((1. - e) * (1. + e));
0038   if (b == 1.) return kHalfPi;
0039   if (b == 0.) return 1.;
0040
0041   Precision x = 1.;
0042   Precision y = b;
0043   Precision S = 0.;
0044   Precision M = 1.;
0045   while (x - y > eps * y) {
0046     Precision tmp = (x + y) * 0.5;
0047
0048     y = vecCore::math::Sqrt(x * y);
0049     x = tmp;
0050     M += M;
0051     S += M * (x - y) * (x - y);
0052   }
0053   return 0.5 * kHalfPi * ((a + b) * (a + b) - S) / (x + y);
0054 }
0055
0056 /// Computes perimeter of ellipse (X/A)^2 + (Y/B)^2 = 1
0057 /// @param A X semi-axis
0058 /// @param B Y semi-axis
0059 //
0060 VECCORE_ATT_HOST_DEVICE
0061 VECGEOM_FORCE_INLINE
0062 Precision EllipsePerimeter(Precision A, Precision B)
0063 {
0064   Precision dx  = vecCore::math::Abs(A);
0065   Precision dy  = vecCore::math::Abs(B);
0066   Precision a   = vecCore::math::Max(dx, dy);
0067   Precision b   = vecCore::math::Min(dx, dy);
0068   Precision b_a = b / a;
0069   Precision e   = vecCore::math::Sqrt((1. - b_a) * (1. + b_a));
0070   return 4. * a * comp_ellint_2(e);
0071 }
0072
0073 /// Computes the lateral surface area of elliptic cone
0074 /// @param pA X semi-axis of ellipse at base (Z = 0)
0075 /// @param pB Y semi-axis of ellipse at base (Z = 0)
0076 /// @param pH height
0077 //
0078 VECCORE_ATT_HOST_DEVICE
0079 VECGEOM_FORCE_INLINE
0080 Precision EllipticalConeLateralArea(Precision pA, Precision pB, Precision pH)
0081 {
0082   Precision x  = vecCore::math::Abs(pA);
0083   Precision y  = vecCore::math::Abs(pB);
0084   Precision h  = vecCore::math::Abs(pH);
0085   Precision a  = vecCore::math::Max(x, y);
0086   Precision b  = vecCore::math::Min(x, y);
0087   Precision aa = a * a;
0088   Precision bb = b * b;
0089   Precision hh = h * h;
0090   Precision e  = vecCore::math::Sqrt(((a - b) * (a + b) * hh) / ((hh + bb) * aa));
0091   return 2. * a * vecCore::math::Sqrt(hh + bb) * comp_ellint_2(e);
0092 }
0093
0094 /// Picks random point inside ellipse (x/a)^2 + (y/b)^2 = 1
0095 /// @param a X semi-axis
0096 /// @param b Y semi-axis
0097 //
0098 VECCORE_ATT_HOST_DEVICE
0099 VECGEOM_FORCE_INLINE
0100 Vector2D<Precision> RandomPointInEllipse(Precision a, Precision b)
0101 {
0102   Precision aa = (a * a == 0.) ? 0. : 1. / (a * a);
0103   Precision bb = (b * b == 0.) ? 0. : 1. / (b * b);
0104   for (int i = 0; i < 1000; ++i) {
0105     Precision x = a * (2. * RNG::Instance().uniform() - 1.);
0106     Precision y = b * (2. * RNG::Instance().uniform() - 1.);
0107     if (x * x * aa + y * y * bb <= 1.) return Vector2D<Precision>(x, y);
0108   }
0109   return Vector2D<Precision>(0., 0.);
0110 }
0111
0112 /// Picks random point on ellipse (x/a)^2 + (y/b)^2 = 1
0113 /// @param a X semi-axis
0114 /// @param b Y semi-axis
0115 //
0116 VECCORE_ATT_HOST_DEVICE
0117 VECGEOM_FORCE_INLINE
0118 Vector2D<Precision> RandomPointOnEllipse(Precision a, Precision b)
0119 {
0120   Precision A      = vecCore::math::Abs(a);
0121   Precision B      = vecCore::math::Abs(b);
0122   Precision mu_max = vecCore::math::Max(A, B);
0123
0124   Precision x, y;
0125   for (int i = 0; i < 1000; ++i) {
0126     Precision phi = kTwoPi * RNG::Instance().uniform();
0127     x             = vecCore::math::Cos(phi);
0128     y             = vecCore::math::Sin(phi);
0129     Precision Bx  = B * x;
0130     Precision Ay  = A * y;
0131     Precision mu  = vecCore::math::Sqrt(Bx * Bx + Ay * Ay);
0132     if (mu_max * RNG::Instance().uniform() <= mu) break;
0133   }
0134   return Vector2D<Precision>(A * x, B * y);
0135 }
0136
0137 } // namespace EllipticUtilities
0138 } // namespace VECGEOM_IMPL_NAMESPACE
0139 } // namespace vecgeom
0140
0141 #endif /* VECGEOM_ELLIPTICUTILITIES_H_ */