surf/detail/QuadraticSolver.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2021-2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file orange/surf/detail/QuadraticSolver.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <cmath>
0011
0012 #include "corecel/Types.hh"
0013 #include "corecel/cont/Array.hh"
0014 #include "corecel/math/Algorithms.hh"
0015 #include "orange/OrangeTypes.hh"
0016
0017 namespace celeritas
0018 {
0019 namespace detail
0020 {
0021 //---------------------------------------------------------------------------//
0022 /*!
0023  * Find positive, real, nonzero roots for quadratic functions.
0024  *
0025  * The quadratic equation \f[
0026    a x^2 + b^2 + c = 0
0027  * \f]
0028  * has two solutions mathematically, but we only want solutions where x is real
0029  * and positive.  Furthermore the equation is subject to catastrophic roundoff
0030  * due to floating point precision (see \c Tolerance::sqrt_quadratic and the
0031  * derivation in \c CylAligned ).
0032  *
0033  * \return An Intersections array where each item is a positive valid
0034  * intersection or the sentinel result \c no_intersection() .
0035  */
0036 class QuadraticSolver
0037 {
0038   public:
0039     //!@{
0040     //! \name Type aliases
0041     using Intersections = Array<real_type, 2>;
0042     //!@}
0043
0044     // Solve when possibly along a surface (zeroish a)
0045     static inline CELER_FUNCTION Intersections solve_general(
0046         real_type a, real_type half_b, real_type c, SurfaceState on_surface);
0047
0048     // Solve degenerate case when a ~ 0 but not on surface
0049     static inline CELER_FUNCTION Intersections
0050     solve_along_surface(real_type half_b, real_type c);
0051
0052   public:
0053     // Construct with nonzero a, and b/2
0054     inline CELER_FUNCTION QuadraticSolver(real_type a, real_type half_b);
0055
0056     // Solve fully general case
0057     inline CELER_FUNCTION Intersections operator()(real_type c) const;
0058
0059     // Solve degenerate case (known to be on surface)
0060     inline CELER_FUNCTION Intersections operator()() const;
0061
0062   private:
0063     //// DATA ////
0064     real_type a_inv_;
0065     real_type hba_;  // (b/2)/a
0066
0067     //! Fuzziness for "along surface" intercept
0068     static CELER_CONSTEXPR_FUNCTION real_type min_a()
0069     {
0070         return ipow<2>(Tolerance<>::sqrt_quadratic());
0071     }
0072 };
0073
0074 //---------------------------------------------------------------------------//
0075 // INLINE DEFINITIONS
0076 //---------------------------------------------------------------------------//
0077 /*!
0078  * Find all positive roots for general quadric surfaces.
0079  *
0080  * This is used for cones, simple quadrics, and general quadrics.
0081  */
0082 CELER_FUNCTION auto
0083 QuadraticSolver::solve_general(real_type a,
0084                                real_type half_b,
0085                                real_type c,
0086                                SurfaceState on_surface) -> Intersections
0087 {
0088     if (std::fabs(a) >= QuadraticSolver::min_a())
0089     {
0090         // Not along the surface
0091         QuadraticSolver solve(a, half_b);
0092         return on_surface == SurfaceState::on ? solve() : solve(c);
0093     }
0094     else if (on_surface == SurfaceState::off)
0095     {
0096         // Travelling parallel to the quadric's surface
0097         return QuadraticSolver::solve_along_surface(half_b, c);
0098     }
0099
0100     // On and along surface: no intersections
0101     return {no_intersection(), no_intersection()};
0102 }
0103
0104 //---------------------------------------------------------------------------//
0105 /*!
0106  * Solve degenerate case when parallel to but not on surface.
0107  *
0108  * This is the case when a is approximately zero.
0109  *
0110  * Note that as both a and b -> 0, |x| -> infinity. Thus, if b/a is
0111  * sufficiently small, no positive root is returned (because this case
0112  * corresponds to a ray crossing a surface at an extreme distance).
0113  */
0114 CELER_FUNCTION auto
0115 QuadraticSolver::solve_along_surface(real_type half_b,
0116                                      real_type c) -> Intersections
0117 {
0118     Intersections result;
0119     if (std::fabs(half_b) > QuadraticSolver::min_a())
0120     {
0121         // On and along surface: no intersections
0122         result[0] = -c / (2 * half_b);
0123         if (result[0] < 0)
0124         {
0125             result[0] = no_intersection();
0126         }
0127         result[1] = no_intersection();
0128     }
0129     else
0130     {
0131         result = {no_intersection(), no_intersection()};
0132     }
0133
0134     // On and along surface: no intersections
0135     return result;
0136 }
0137
0138 //---------------------------------------------------------------------------//
0139 /*!
0140  * Construct with non-degenerate a and  b/2.
0141  */
0142 CELER_FUNCTION QuadraticSolver::QuadraticSolver(real_type a, real_type half_b)
0143     : a_inv_(1 / a), hba_(half_b * a_inv_)
0144 {
0145     CELER_EXPECT(std::fabs(a) >= QuadraticSolver::min_a());
0146 }
0147
0148 //---------------------------------------------------------------------------//
0149 /*!
0150  * Find all positive roots of x^2 + (b/a)*x + (c/a) = 0.
0151  *
0152  * In this case, the quadratic formula can be written as: \f[
0153    x = -b/2 \pm \sqrt{(b/2)^2 - c}.
0154  * \f]
0155  *
0156  * Callers:
0157  * - General quadratic solve: not on nor along surface
0158  * - Sphere when not on surface
0159  * - Cylinder when not on surface
0160  */
0161 CELER_FUNCTION auto
0162 QuadraticSolver::operator()(real_type c) const -> Intersections
0163 {
0164     // Scale c by 1/a in accordance with scaling of b
0165     c *= a_inv_;
0166     real_type b2_4 = ipow<2>(hba_);  // (b/2)^2
0167
0168     Intersections result;
0169
0170     if (b2_4 > c)
0171     {
0172         // Two real roots, r1 and r2
0173         real_type t2 = std::sqrt(b2_4 - c);  // (b^2 - 4ac) / 4
0174         result[0] = -hba_ - t2;
0175         result[1] = -hba_ + t2;
0176
0177         if (result[1] <= 0)
0178         {
0179             // Both are nonpositive
0180             result[0] = no_intersection();
0181             result[1] = no_intersection();
0182         }
0183         else if (result[0] <= 0)
0184         {
0185             // Only first is nonpositive
0186             result[0] = no_intersection();
0187         }
0188     }
0189     else if (b2_4 == c)
0190     {
0191         // One real root, r1
0192         result[0] = -hba_;
0193         result[1] = no_intersection();
0194
0195         if (result[0] <= 0)
0196         {
0197             result[0] = no_intersection();
0198         }
0199     }
0200     else
0201     {
0202         // No real roots
0203         result = {no_intersection(), no_intersection()};
0204     }
0205
0206     return result;
0207 }
0208
0209 //---------------------------------------------------------------------------//
0210 /*!
0211  * Solve degenerate case (known to be "on" surface).
0212  *
0213  * Since x = 0 is a root, then c = 0 and x = -b is the other root. This will be
0214  * inaccurate if a particle is logically on a surface but not physically on it.
0215  */
0216 CELER_FUNCTION auto QuadraticSolver::operator()() const -> Intersections
0217 {
0218     Intersections result{-2 * hba_, no_intersection()};
0219
0220     if (result[0] <= 0)
0221     {
0222         result[0] = no_intersection();
0223     }
0224
0225     return result;
0226 }
0227
0228 //---------------------------------------------------------------------------//
0229 }  // namespace detail
0230 }  // namespace celeritas