surf/detail/QuadricConeConverter.hh

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 orange/surf/detail/QuadricConeConverter.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <optional>
0010
0011 #include "corecel/math/SoftEqual.hh"
0012
0013 #include "../ConeAligned.hh"
0014 #include "../SimpleQuadric.hh"
0015
0016 namespace celeritas
0017 {
0018 namespace detail
0019 {
0020 //---------------------------------------------------------------------------//
0021 /*!
0022  * Try to convert a simple quadric to a cone.
0023  *
0024  * The simple quadric must have already undergone normalization (one
0025  * second-order term less than zero, two positive).
0026  */
0027 class QuadricConeConverter
0028 {
0029   public:
0030     // Construct with tolerance
0031     inline QuadricConeConverter(real_type tol);
0032
0033     // Try converting to a cone with this orientation
0034     template<Axis T>
0035     std::optional<ConeAligned<T>>
0036     operator()(AxisTag<T>, SimpleQuadric const& sq) const;
0037
0038   private:
0039     SoftEqual<> soft_equal_;
0040 };
0041
0042 //---------------------------------------------------------------------------//
0043 // INLINE DEFINITIONS
0044 //---------------------------------------------------------------------------//
0045 /*!
0046  * Construct with tolerance.
0047  */
0048 QuadricConeConverter::QuadricConeConverter(real_type tol) : soft_equal_{tol} {}
0049
0050 //---------------------------------------------------------------------------//
0051 /*!
0052  * Try converting to a cone with this orientation.
0053  *
0054  * Cone:
0055  * \verbatim
0056     -t^2 (x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = 0
0057  * \endverbatim
0058  * Expanded:
0059  * \verbatim
0060            -t^2 x^2    + y^2    + z^2
0061     + 2 t^2 x_0 x - 2y_0 y - 2z_0 z
0062     + (- t^2 x_0^2 + y_0^2 + z_0^2) = 0
0063  * \endverbatim
0064
0065  * SQ:
0066  * \verbatim
0067    a x^2 + b y^2 + (b + \epsilon) z^2 + e x + f y + g z  + h = 0
0068  * \endverbatim
0069  * Normalized (\c epsilon = 0)
0070  * \verbatim
0071    a/b x^2 + y^2 + z^2 + e/b x + f/b y + g/b z  + h/b = 0
0072  * \endverbatim
0073  * Match terms:
0074  * \verbatim
0075         -t^2  = a/b
0076     2 t^2 x_0 = e/b --> x_0 = e / (2 * a)
0077        -2 y_0 = f/b --> y_0 = f / (-2 * b)
0078        -2 z_0 = g/b --> z_0 = g / (-2 * b)
0079    - t^2 x_0^2 + y_0^2 + z_0^2 = h/b    [if not, it's a hyperboloid!]
0080  * \endverbatim
0081  */
0082 template<Axis T>
0083 std::optional<ConeAligned<T>>
0084 QuadricConeConverter::operator()(AxisTag<T>, SimpleQuadric const& sq) const
0085 {
0086     // Other coordinate system
0087     constexpr auto U = ConeAligned<T>::u_axis();
0088     constexpr auto V = ConeAligned<T>::v_axis();
0089
0090     auto second = sq.second();
0091     if (!(second[to_int(T)] < 0))
0092     {
0093         // Not the negative component
0094         return {};
0095     }
0096     if (!soft_equal_(second[to_int(U)], second[to_int(V)]))
0097     {
0098         // Not a circular cone
0099         return {};
0100     }
0101     CELER_ASSERT(second[to_int(U)] > 0 && second[to_int(V)] > 0);
0102
0103     // Normalize so U, V second-order coefficients are 1
0104     real_type const norm = (second[to_int(U)] + second[to_int(V)]) / 2;
0105     real_type const tsq = -second[to_int(T)] / norm;
0106
0107     // Calculate origin from first-order coefficients
0108     Real3 origin = make_array(sq.first());
0109     origin[to_int(T)] /= -2 * second[to_int(T)];
0110     origin[to_int(U)] /= -2 * norm;
0111     origin[to_int(V)] /= -2 * norm;
0112
0113     real_type const expected_h_b = -tsq * ipow<2>(origin[to_int(T)])
0114                                    + ipow<2>(origin[to_int(U)])
0115                                    + ipow<2>(origin[to_int(V)]);
0116     if (!soft_equal_(expected_h_b, sq.zeroth() / norm))
0117     {
0118         // Leftover constant: it's a hyperboloid!
0119         return {};
0120     }
0121
0122     // Clear potential signed zeros before returning
0123     origin += real_type{0};
0124     return ConeAligned<T>::from_tangent_sq(origin, tsq);
0125 }
0126
0127 //---------------------------------------------------------------------------//
0128 }  // namespace detail
0129 }  // namespace celeritas