orange/surf/ConeAligned.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/ConeAligned.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Types.hh"
0010 #include "corecel/cont/Array.hh"
0011 #include "corecel/cont/Span.hh"
0012 #include "corecel/math/ArrayUtils.hh"
0013 #include "orange/OrangeTypes.hh"
0014 #include "orange/SenseUtils.hh"
0015
0016 #include "detail/QuadraticSolver.hh"
0017
0018 namespace celeritas
0019 {
0020 //---------------------------------------------------------------------------//
0021 /*!
0022  * Axis-aligned cone (infinite and double-sheeted).
0023  *
0024  * For a cone parallel to the x axis:
0025  * \f[
0026     (y - y_0)^2 + (z - z_0)^2 - t^2 (x - x_0)^2 = 0
0027    \f]
0028
0029  * where \em t is the tangent of the opening angle (\f$r/h\f$ for a finite cone
0030  * with radius \em r and height \em h). In the cone below, the tangent of the
0031  * inner angle is \f$ t = \tan(\theta) = r/h \f$. Negative values of \em t are
0032  * equivalent and valid, representing the four points on a slice of the
0033  * double-sheet cone.
0034  * \verbatim
0035   \        |   r    /
0036    o- - - -+- - - -o
0037     '--_   :   _--'
0038  h      --.:.--
0039  ---       O
0040        _--':'--_
0041     .--    :    --.
0042    o-------+-------o
0043   /        |        \
0044    \endverbatim
0045  *
0046  */
0047 template<Axis T>
0048 class ConeAligned
0049 {
0050   public:
0051     //@{
0052     //! \name Type aliases
0053     using Intersections = Array<real_type, 2>;
0054     using StorageSpan = Span<real_type const, 4>;
0055     //@}
0056
0057   private:
0058     static constexpr Axis U{T == Axis::x ? Axis::y : Axis::x};
0059     static constexpr Axis V{T == Axis::z ? Axis::y : Axis::z};
0060
0061   public:
0062     //// CLASS ATTRIBUTES ////
0063
0064     // Surface type identifier
0065     static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type();
0066
0067     //! Safety is intersection along surface normal
0068     static CELER_CONSTEXPR_FUNCTION bool simple_safety() { return false; }
0069
0070     //!@{
0071     //! Axes
0072     static CELER_CONSTEXPR_FUNCTION Axis t_axis() { return T; }
0073     static CELER_CONSTEXPR_FUNCTION Axis u_axis() { return U; }
0074     static CELER_CONSTEXPR_FUNCTION Axis v_axis() { return V; }
0075     //!@}
0076
0077   public:
0078     //// CONSTRUCTORS ////
0079
0080     // Construct with square of tangent for simplification
0081     static ConeAligned from_tangent_sq(Real3 const& origin, real_type tsq);
0082
0083     // Construct from origin and tangent of the angle of its opening
0084     ConeAligned(Real3 const& origin, real_type tangent);
0085
0086     // Construct from raw data
0087     template<class R>
0088     explicit inline CELER_FUNCTION ConeAligned(Span<R, StorageSpan::extent>);
0089
0090     //// ACCESSORS ////
0091
0092     //! Get the origin position along the normal axis
0093     CELER_FUNCTION Real3 const& origin() const { return origin_; }
0094
0095     //! Get the square of the tangent of the opening angle
0096     CELER_FUNCTION real_type tangent_sq() const { return tsq_; }
0097
0098     //! Get a view to the data for type-deleted storage
0099     CELER_FUNCTION StorageSpan data() const { return {&origin_[0], 4}; }
0100
0101     //// CALCULATION ////
0102
0103     // Determine the sense of the position relative to this surface
0104     inline CELER_FUNCTION SignedSense calc_sense(Real3 const& pos) const;
0105
0106     // Calculate all possible straight-line intersections with this surface
0107     inline CELER_FUNCTION Intersections calc_intersections(
0108         Real3 const& pos, Real3 const& dir, SurfaceState on_surface) const;
0109
0110     // Calculate outward normal at a position
0111     inline CELER_FUNCTION Real3 calc_normal(Real3 const& pos) const;
0112
0113   private:
0114     // Location of the vanishing point
0115     Real3 origin_;
0116
0117     // Quadric value
0118     real_type tsq_;
0119
0120     //! Private default constructor for manual construction
0121     ConeAligned() = default;
0122 };
0123
0124 //---------------------------------------------------------------------------//
0125 // TYPE ALIASES
0126 //---------------------------------------------------------------------------//
0127
0128 using ConeX = ConeAligned<Axis::x>;
0129 using ConeY = ConeAligned<Axis::y>;
0130 using ConeZ = ConeAligned<Axis::z>;
0131
0132 //---------------------------------------------------------------------------//
0133 // INLINE DEFINITIONS
0134 //---------------------------------------------------------------------------//
0135 /*!
0136  * Surface type identifier.
0137  */
0138 template<Axis T>
0139 CELER_CONSTEXPR_FUNCTION SurfaceType ConeAligned<T>::surface_type()
0140 {
0141     return T == Axis::x   ? SurfaceType::kx
0142            : T == Axis::y ? SurfaceType::ky
0143            : T == Axis::z ? SurfaceType::kz
0144                           : SurfaceType::size_;
0145 }
0146
0147 //---------------------------------------------------------------------------//
0148 /*!
0149  * Construct from raw data.
0150  */
0151 template<Axis T>
0152 template<class R>
0153 CELER_FUNCTION ConeAligned<T>::ConeAligned(Span<R, StorageSpan::extent> data)
0154     : origin_{data[0], data[1], data[2]}, tsq_{data[3]}
0155 {
0156 }
0157
0158 //---------------------------------------------------------------------------//
0159 /*!
0160  * Determine the sense of the position relative to this surface.
0161  */
0162 template<Axis T>
0163 CELER_FUNCTION SignedSense ConeAligned<T>::calc_sense(Real3 const& pos) const
0164 {
0165     real_type const x = pos[to_int(T)] - origin_[to_int(T)];
0166     real_type const y = pos[to_int(U)] - origin_[to_int(U)];
0167     real_type const z = pos[to_int(V)] - origin_[to_int(V)];
0168
0169     return real_to_sense((-tsq_ * ipow<2>(x)) + ipow<2>(y) + ipow<2>(z));
0170 }
0171
0172 //---------------------------------------------------------------------------//
0173 /*!
0174  * Calculate all possible straight-line intersections with this surface.
0175  *
0176  * \f[
0177     (y - yc)^2 + (z - zc)^2 - t^2 * (x - xc)^2 = 0
0178    \f]
0179  */
0180 template<Axis T>
0181 CELER_FUNCTION auto
0182 ConeAligned<T>::calc_intersections(Real3 const& pos,
0183                                    Real3 const& dir,
0184                                    SurfaceState on_surface) const
0185     -> Intersections
0186 {
0187     // Expand translated positions into 'xyz' coordinate system
0188     real_type const x = pos[to_int(T)] - origin_[to_int(T)];
0189     real_type const y = pos[to_int(U)] - origin_[to_int(U)];
0190     real_type const z = pos[to_int(V)] - origin_[to_int(V)];
0191
0192     real_type const u = dir[to_int(T)];
0193     real_type const v = dir[to_int(U)];
0194     real_type const w = dir[to_int(V)];
0195
0196     // Scaled direction
0197     real_type a = (-tsq_ * ipow<2>(u)) + ipow<2>(v) + ipow<2>(w);
0198     real_type half_b = (-tsq_ * x * u) + (y * v) + (z * w);
0199     real_type c = (-tsq_ * ipow<2>(x)) + ipow<2>(y) + ipow<2>(z);
0200
0201     return detail::QuadraticSolver::solve_general(a, half_b, c, on_surface);
0202 }
0203
0204 //---------------------------------------------------------------------------//
0205 /*!
0206  * Calculate outward normal at a position.
0207  */
0208 template<Axis T>
0209 CELER_FUNCTION Real3 ConeAligned<T>::calc_normal(Real3 const& pos) const
0210 {
0211     Real3 norm;
0212     for (auto i = to_int(Axis::x); i < to_int(Axis::size_); ++i)
0213     {
0214         norm[i] = pos[i] - origin_[i];
0215     }
0216     norm[to_int(T)] *= -tsq_;
0217
0218     return make_unit_vector(norm);
0219 }
0220
0221 //---------------------------------------------------------------------------//
0222 }  // namespace celeritas