Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-05-09 08:52:18

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/orangeinp/IntersectRegion.hh
0006 //! \brief Contains IntersectRegionInterface and concrete daughters
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009 
0010 #include <vector>
0011 
0012 #include "corecel/Macros.hh"
0013 #include "corecel/cont/Array.hh"
0014 #include "corecel/cont/EnumArray.hh"
0015 #include "corecel/grid/GridTypes.hh"
0016 #include "corecel/math/Turn.hh"
0017 #include "orange/OrangeTypes.hh"
0018 
0019 namespace celeritas
0020 {
0021 struct JsonPimpl;
0022 
0023 namespace orangeinp
0024 {
0025 class IntersectSurfaceBuilder;
0026 
0027 //---------------------------------------------------------------------------//
0028 /*!
0029  * Interface class for building non-reentrant spatial regions.
0030  *
0031  * This is a building block for constructing more complex objects out of
0032  * smaller spatial regions. A \c shape object will have a single intersect
0033  * region, and a \c solid object region may have multiple adjacent intersect
0034  * regions.
0035  *
0036  * Convex regions should be as minimal as possible and rely on transformations
0037  * to change axes, displacement, etc. As a general rule, the exterior bounding
0038  * box of a intersect region should be <em>centered on the origin</em>, and
0039  * objects should be aligned along the \em z axis.
0040  *
0041  * When implementing this class, prefer to build simpler surfaces (planes)
0042  * before complex ones (cones) in case we implement short-circuiting logic,
0043  * since expressions are currently sorted.
0044  *
0045  * \note Additional methods such as volume calculation may be added here later.
0046  */
0047 class IntersectRegionInterface
0048 {
0049   public:
0050     //! Construct surfaces that are AND-ed into this region
0051     virtual void build(IntersectSurfaceBuilder&) const = 0;
0052 
0053     //! Write the region to a JSON object
0054     virtual void output(JsonPimpl*) const = 0;
0055 
0056     virtual ~IntersectRegionInterface() = default;
0057 
0058   protected:
0059     //!@{
0060     //! Allow construction and assignment only through daughter classes
0061     IntersectRegionInterface() = default;
0062     CELER_DEFAULT_COPY_MOVE(IntersectRegionInterface);
0063     //!@}
0064 };
0065 
0066 //---------------------------------------------------------------------------//
0067 /*!
0068  * A rectangular parallelepiped/cuboid centered on the origin.
0069  *
0070  * The box is constructed with half-widths.
0071  */
0072 class Box final : public IntersectRegionInterface
0073 {
0074   public:
0075     // Construct with half-widths
0076     explicit Box(Real3 const& halfwidths);
0077 
0078     // Build surfaces
0079     void build(IntersectSurfaceBuilder&) const final;
0080 
0081     // Output to JSON
0082     void output(JsonPimpl*) const final;
0083 
0084     //// ACCESSORS ////
0085 
0086     //! Half-width for each axis
0087     Real3 const& halfwidths() const { return hw_; }
0088 
0089   private:
0090     Real3 hw_;
0091 };
0092 
0093 //---------------------------------------------------------------------------//
0094 /*!
0095  * A closed truncated cone along the *z* axis centered on the origin.
0096  *
0097  * A quadric cone technically defines two opposing cones that touch at a single
0098  * vanishing point, but this cone is required to be truncated so that the
0099  * vanishing point is on our outside the cone.
0100  *
0101  * The midpoint along the \em z axis of the cone is the origin. A cone is \em
0102  * not allowed to have equal radii: for that, use a cylinder. However, it \em
0103  * may have a single radius of zero, which puts the vanishing point on one end
0104  * of the cone.
0105  *
0106  * This intersect region, along with the Cylinder, is a base component of the
0107  * G4Polycone (PCON).
0108  */
0109 class Cone final : public IntersectRegionInterface
0110 {
0111   public:
0112     // Construct with Z half-height and lo, hi radii
0113     Cone(Real2 const& radii, real_type halfheight);
0114 
0115     //// INTERFACE ////
0116 
0117     // Build surfaces
0118     void build(IntersectSurfaceBuilder&) const final;
0119 
0120     // Output to JSON
0121     void output(JsonPimpl*) const final;
0122 
0123     //// TEMPLATE INTERFACE ////
0124 
0125     // Whether this encloses another cone
0126     bool encloses(Cone const& other) const;
0127 
0128     //// ACCESSORS ////
0129 
0130     //! Lower and upper radii
0131     Real2 const& radii() const { return radii_; }
0132 
0133     //! Half-height along Z
0134     real_type halfheight() const { return hh_; }
0135 
0136   private:
0137     Real2 radii_;
0138     real_type hh_;
0139 };
0140 
0141 //---------------------------------------------------------------------------//
0142 /*!
0143  * A *z*-aligned cylinder centered on the origin.
0144  *
0145  * The cylinder is defined with a radius and half-height.
0146  */
0147 class Cylinder final : public IntersectRegionInterface
0148 {
0149   public:
0150     // Construct with radius
0151     Cylinder(real_type radius, real_type halfheight);
0152 
0153     // Build surfaces
0154     void build(IntersectSurfaceBuilder&) const final;
0155 
0156     // Output to JSON
0157     void output(JsonPimpl*) const final;
0158 
0159     //// TEMPLATE INTERFACE ////
0160 
0161     // Whether this encloses another cylinder
0162     bool encloses(Cylinder const& other) const;
0163 
0164     //// ACCESSORS ////
0165 
0166     //! Radius
0167     real_type radius() const { return radius_; }
0168 
0169     //! Half-height along Z
0170     real_type halfheight() const { return hh_; }
0171 
0172   private:
0173     real_type radius_;
0174     real_type hh_;
0175 };
0176 
0177 //---------------------------------------------------------------------------//
0178 /*!
0179  * An axis-aligned ellipsoid centered at the origin.
0180  *
0181  * The ellipsoid is constructed with the three radial lengths. For a length
0182  * scale \em L , the quadric it creates has second-order terms that are \f$
0183  * O(1) \f$ and a zeroth order term that's \f$ O(L^2) \f$. Translations on that
0184  * length scale will preserve the accuracy of the quadratic solution.
0185  *
0186  * There are many scalings of the quadric equation that produce unitary
0187  * second-order terms if the ellipsoid's radii are identical: \f[
0188   \frac{k}{r_x^2} x^2 +  \frac{k}{r_y^2} y^2 +  \frac{k}{r_z^2} z^2 = k
0189  * \f]
0190  * but we make the ad hoc decision to choose \f$ k = \min r_i \max r_i \f$
0191  * to avoid irrational normalization constants, which makes unit tests and
0192  * output easier to read.
0193  */
0194 class Ellipsoid final : public IntersectRegionInterface
0195 {
0196   public:
0197     // Construct with radius along each Cartesian axis
0198     explicit Ellipsoid(Real3 const& radii);
0199 
0200     // Build surfaces
0201     void build(IntersectSurfaceBuilder&) const final;
0202 
0203     // Output to JSON
0204     void output(JsonPimpl*) const final;
0205 
0206     //// TEMPLATE INTERFACE ////
0207 
0208     // Whether this encloses another ellipsoid
0209     bool encloses(Ellipsoid const& other) const;
0210 
0211     //// ACCESSORS ////
0212 
0213     //! Radius along each axis
0214     Real3 const& radii() const { return radii_; }
0215 
0216     //! Get the radius along a single axis
0217     real_type radius(Axis ax) const { return radii_[to_int(ax)]; }
0218 
0219   private:
0220     Real3 radii_;
0221 };
0222 
0223 //---------------------------------------------------------------------------//
0224 /*!
0225  * A *z*-aligned cylinder with an elliptical cross section.
0226  *
0227  * The elliptical cylinder is defined with a two radii and a half-height,
0228  * such that the centroid of the bounding box is origin. The quadric
0229  * coefficient of the cylindrical component, \f[
0230   x^2 / r_x^2 + y^2 / r_y^2 = 1 \,
0231   \f]
0232  * are scaled based on the radii of the cylinder, reducing to
0233  * \f$ x^2 + y^2 = R^2 \f$ when the radii are equal.
0234  */
0235 class EllipticalCylinder final : public IntersectRegionInterface
0236 {
0237   public:
0238     // Construct with x- and y-radii and half-height in z
0239     EllipticalCylinder(Real2 const& radii, real_type halfheight);
0240 
0241     // Build surfaces
0242     void build(IntersectSurfaceBuilder&) const final;
0243 
0244     // Output to JSON
0245     void output(JsonPimpl*) const final;
0246 
0247     //// TEMPLATE INTERFACE ////
0248 
0249     // Whether this encloses another ellipsoid
0250     bool encloses(EllipticalCylinder const& other) const;
0251 
0252     //// ACCESSORS ////
0253 
0254     //! Radius along each axis
0255     Real2 const& radii() const { return radii_; }
0256 
0257     //! Half-height along Z
0258     real_type halfheight() const { return hh_; }
0259 
0260     // Get the radius along the x/y axis
0261     real_type radius(Axis ax) const;
0262 
0263   private:
0264     Real2 radii_;
0265     real_type hh_;
0266 };
0267 
0268 //---------------------------------------------------------------------------//
0269 /*!
0270  * A finite *z*-aligned cone with an elliptical cross section.
0271  *
0272  * The elliptical cone is defined in an analogous fashion to the regular
0273  * (i.e., circular) cone. A half-height (hh) defines the z-extents, such
0274  * that the centroid of the outer bounding box is the origin. The lower radii
0275  * are the x- and y-radii at the plane z = -hh. The upper radii are the x- and
0276  * y-radii at the plane z = hh. There are several restrictions on these radii:
0277  *
0278  * 1) Either the lower or upper radii may be (0, 0); this is the only permitted
0279  *    way for the elliptical cone to include the vertex.
0280  * 2) The aspect ratio of the elliptical cross sections is constant. Thus, the
0281  *    aspect ratio at z = -hh must equal the aspect ratio at z = hh.
0282  * 3) Degenerate elliptical cones with lower_radii == upper_radii (i.e.,
0283  *    elliptical cylinders) are not permitted.
0284  * 4) Degenerate elliptical cones where lower or upper radii are equal to
0285  *    (0, x) or (x, 0), where x is non-zero, are not permitted.
0286  *
0287  * The elliptical surface can be expressed as:
0288  *
0289  * \f[
0290    (x/r_x)^2 + (y/r_y)^2 = (v-z)^2,
0291  * \f]
0292  *
0293  * where v is the location of the vertex. The r_x, r_y, and v can be calculated
0294  * from the lower and upper radii as given by \c G4EllipticalCone:
0295  * \verbatim
0296    r_x = (lower_radii[X] - upper_radii[X])/(2 hh),
0297    r_y = (lower_radii[Y] - upper_radii[Y])/(2 hh),
0298      v = hh (lower_radii[X] + upper_radii[X])/(lower_radii[X] -
0299  upper_radii[X]).
0300    \endverbatim
0301  */
0302 class EllipticalCone final : public IntersectRegionInterface
0303 {
0304   public:
0305     // Construct with x- and y-radii and half-height in z
0306     EllipticalCone(Real2 const& lower_radii,
0307                    Real2 const& upper_radii,
0308                    real_type halfheight);
0309 
0310     // Build surfaces
0311     void build(IntersectSurfaceBuilder&) const final;
0312 
0313     // Output to JSON
0314     void output(JsonPimpl*) const final;
0315 
0316     //// TEMPLATE INTERFACE ////
0317 
0318     // Whether this encloses another elliptical cone
0319     bool encloses(EllipticalCone const& other) const;
0320 
0321     //// ACCESSORS ////
0322 
0323     //! Radii along the x- and y-axes at z=-hh
0324     Real2 const& lower_radii() const { return lower_radii_; }
0325 
0326     //! Radii along the x- and y-axes at z=-hh
0327     Real2 const& upper_radii() const { return upper_radii_; }
0328 
0329     //! Half-height along Z
0330     real_type halfheight() const { return hh_; }
0331 
0332     // Get the bottom/top radius along the x/y axis
0333     real_type radius(Bound b, Axis ax) const;
0334 
0335   private:
0336     Real2 lower_radii_;
0337     Real2 upper_radii_;
0338     real_type hh_;
0339 };
0340 
0341 //---------------------------------------------------------------------------//
0342 /*!
0343  * Region formed by extruding + scaling a convex polygon along a line segment.
0344  *
0345  * The convex polygon is supplied as a set of points on the XY plane in
0346  * counterclockwise order. The line segment and scaling factors are specified
0347  * by providing a line segment point and scaling factor for the top and bottom
0348  * polygon faces of the region. The line segment point of the top face must
0349  * have a z value greater than that of the bottom face. Along the line segment,
0350  * the size of the polygon is linearly scaled in accordance with scaling
0351  * factors.
0352  *
0353  * As is done in Geant4, construction is done by first applying scaling factors
0354  * to the upper and lower polygons via scalar multiplication with each polygon
0355  * point, then the points on the line are used to offset the upper and lower
0356  * polygons.
0357  */
0358 class ExtrudedPolygon final : public IntersectRegionInterface
0359 {
0360   public:
0361     //!@{
0362     //! \name Type aliases
0363     using VecReal2 = std::vector<Real2>;
0364 
0365     //! Specifies the top or bottom face of the ExtrudedPolygon
0366     struct PolygonFace
0367     {
0368         //! Start or end point of the line segment the polygon is extruded
0369         //! along
0370         Real3 line_segment_point{};
0371 
0372         //! The fractional amount this face is scaled
0373         real_type scaling_factor{};
0374     };
0375     //!@}
0376 
0377   public:
0378     // Construct from a convex polygon and bottom/top faces
0379     ExtrudedPolygon(VecReal2 const& polygon,
0380                     PolygonFace const& bot_face,
0381                     PolygonFace const& top_face);
0382 
0383     // Build surfaces
0384     void build(IntersectSurfaceBuilder&) const final;
0385 
0386     // Output to JSON
0387     void output(JsonPimpl*) const final;
0388 
0389     //// ACCESSORS ////
0390 
0391     //! Polygon points (2D)
0392     VecReal2 polygon() const { return polygon_; }
0393 
0394     //! Bottom point of the line segment
0395     Real3 bot_line_segment_point() const { return line_segment_[Bound::lo]; }
0396 
0397     //! Top point of the line segment
0398     Real3 top_line_segment_point() const { return line_segment_[Bound::hi]; }
0399 
0400     //! Bottom scaling factor
0401     real_type bot_scaling_factor() const
0402     {
0403         return scaling_factors_[Bound::lo];
0404     }
0405 
0406     //! Top scaling factor
0407     real_type top_scaling_factor() const
0408     {
0409         return scaling_factors_[Bound::hi];
0410     }
0411 
0412   private:
0413     //// TYPES ////
0414     using Range = celeritas::Array<real_type, 2>;
0415 
0416     //// DATA ////
0417 
0418     VecReal2 polygon_;
0419 
0420     EnumArray<Bound, Real3> line_segment_;
0421     EnumArray<Bound, real_type> scaling_factors_;
0422 
0423     Range x_range_;
0424     Range y_range_;
0425 
0426     //// HELPER FUNCTIONS ////
0427 
0428     // Calculate the min/max x or y values of the extruded region
0429     Range calc_range(VecReal2 const& polygon, size_type dim);
0430 };
0431 
0432 //---------------------------------------------------------------------------//
0433 /*!
0434  * A generalized polygon with parallel flat faces along the *z* axis.
0435  *
0436  * A GenPrism, like VecGeom's GenTrap, ROOT's Arb8, and Geant4's
0437  * G4GenericTrap, represents a generalized volume with polyhedral faces on two
0438  * parallel planes perpendicular to the \em z axis. Unlike those other codes,
0439  * the number of faces can be arbitrary in number.
0440  *
0441  * The faces have an orientation and ordering so that \em twisted faces can be
0442  * constructed by joining corresponding points using straight-line "vertical"
0443  * edges, directly matching the G4GenericTrap definition, but directly
0444  * generating a hyperbolic paraboloid for each twisted face.
0445  *
0446  * Trapezoids constructed from the helper functions will have sides that are
0447  * same ordering as a prism: the rightward face is first (normal is along the
0448  * \em +x axis), then the others follow counterclockwise.
0449  */
0450 class GenPrism final : public IntersectRegionInterface
0451 {
0452   public:
0453     //!@{
0454     //! \name Type aliases
0455     using VecReal2 = std::vector<Real2>;
0456     //!@}
0457 
0458     //! Regular trapezoidal top/bottom face
0459     struct TrapFace
0460     {
0461         //! Half the vertical distance between horizontal edges
0462         real_type hy{};
0463         //! Top horizontal edge half-length
0464         real_type hx_lo{};
0465         //! Bottom horizontal edge half-length
0466         real_type hx_hi{};
0467         //! Shear angle between horizontal line centers and Y axis
0468         Turn alpha;
0469     };
0470 
0471   public:
0472     // Helper function to construct a Trd shape from hz and two rectangles,
0473     // one for each z-face
0474     static GenPrism from_trd(real_type halfz, Real2 const& lo, Real2 const& hi);
0475 
0476     // Helper function to construct a general trap from its half-height and
0477     // the two trapezoids defining its lower and upper faces
0478     static GenPrism from_trap(real_type hz,
0479                               Turn theta,
0480                               Turn phi,
0481                               TrapFace const& lo,
0482                               TrapFace const& hi);
0483 
0484     // Construct from half Z height and 4 vertices for top and bottom planes
0485     GenPrism(real_type halfz, VecReal2 const& lo, VecReal2 const& hi);
0486 
0487     // Build surfaces
0488     void build(IntersectSurfaceBuilder&) const final;
0489 
0490     // Output to JSON
0491     void output(JsonPimpl*) const final;
0492 
0493     //// ACCESSORS ////
0494 
0495     //! Half-height along Z
0496     real_type halfheight() const { return hh_; }
0497 
0498     //! Polygon on -z face
0499     VecReal2 const& lower() const { return lo_; }
0500 
0501     //! Polygon on +z face
0502     VecReal2 const& upper() const { return hi_; }
0503 
0504     //! Number of sides (points on the Z face)
0505     size_type num_sides() const { return lo_.size(); }
0506 
0507     // Calculate the cosine of the twist angle for a given side
0508     real_type calc_twist_cosine(size_type size_idx) const;
0509 
0510   private:
0511     enum class Degenerate
0512     {
0513         none,
0514         lo,
0515         hi
0516     };
0517 
0518     real_type hh_;  //!< half-height
0519     VecReal2 lo_;  //!< corners of the -z face (CCW)
0520     VecReal2 hi_;  //!< corners of the +z face (CCW)
0521     Degenerate degen_{Degenerate::none};  //!< no plane on this z axis
0522     real_type length_scale_{};
0523 };
0524 
0525 //---------------------------------------------------------------------------//
0526 /*!
0527  * An axis-aligned infinite half-space to use for truncation operations.
0528  *
0529  * An "inside" sense means to include everything *below* the position on the
0530  * axis, and an "outside" sense means to include only what's *above* the
0531  * position.
0532  */
0533 class InfPlane : public IntersectRegionInterface
0534 {
0535   public:
0536     // Construct with sense, axis, and position
0537     InfPlane(Sense sense, Axis axis, real_type position);
0538 
0539     // Build surfaces
0540     void build(IntersectSurfaceBuilder&) const final;
0541 
0542     // Output to JSON
0543     void output(JsonPimpl*) const final;
0544 
0545     //// ACCESSORS ////
0546 
0547     //! Get the sense (inside or outside)
0548     Sense sense() const { return sense_; }
0549 
0550     //! Get the axis (x, y, or z)
0551     Axis axis() const { return axis_; }
0552 
0553     //! Get the position along the axis
0554     real_type position() const { return position_; }
0555 
0556   private:
0557     Sense sense_;
0558     Axis axis_;
0559     real_type position_;
0560 };
0561 
0562 //---------------------------------------------------------------------------//
0563 /*!
0564  * An open wedge shape from the *z* axis.
0565  *
0566  * The wedge is defined by an interior angle that \em must be less than or
0567  * equal to 180 degrees (half a turn) and \em must be more than zero. It can be
0568  * subtracted, or its negation can be subtracted. The start angle is mapped
0569  * onto \f$[0, 1)\f$ on construction.
0570  */
0571 class InfAziWedge final : public IntersectRegionInterface
0572 {
0573   public:
0574     // Construct from an angular range less than 180
0575     InfAziWedge(Turn start, Turn stop);
0576 
0577     // Build surfaces
0578     void build(IntersectSurfaceBuilder&) const final;
0579 
0580     // Output to JSON
0581     void output(JsonPimpl*) const final;
0582 
0583     //// ACCESSORS ////
0584 
0585     //! Starting angle
0586     Turn start() const { return start_; }
0587 
0588     //! stop angle
0589     Turn stop() const { return stop_; }
0590 
0591   private:
0592     Turn start_;
0593     Turn stop_;
0594 };
0595 
0596 //---------------------------------------------------------------------------//
0597 /*!
0598  * Select a polar (latitudinal) region.
0599  *
0600  * This uses an equatorial plane and up to two cones to slice a
0601  * polar-coordinate region from the origin.  A polar wedge always defines a
0602  * region in a single hemisphere: either \f$ z >= 0 \f$ or \f$ z <= 0 \f$,
0603  * corresponding to an stop range of [0, .25] turns or [0.25, 0.5] turns.
0604  */
0605 class InfPolarWedge final : public IntersectRegionInterface
0606 {
0607   public:
0608     // Construct from a starting angle and stop angle
0609     InfPolarWedge(Turn start, Turn stop);
0610 
0611     // Build surfaces
0612     void build(IntersectSurfaceBuilder&) const final;
0613 
0614     // Output to JSON
0615     void output(JsonPimpl*) const final;
0616 
0617     //// ACCESSORS ////
0618 
0619     //! Starting angle
0620     Turn start() const { return start_; }
0621 
0622     //! stop angle
0623     Turn stop() const { return stop_; }
0624 
0625   private:
0626     Turn start_;
0627     Turn stop_;
0628 
0629     static constexpr Turn north_pole{0};
0630     static constexpr Turn equator{0.25};
0631     static constexpr Turn south_pole{0.5};
0632 };
0633 
0634 //---------------------------------------------------------------------------//
0635 /*!
0636  * An involute "blade" centered on the origin.
0637  *
0638  * This is the intersection of two parallel involutes with a cylindrical shell.
0639  * The three radii, which must be in ascending order, are that of the involute,
0640  * the inner cylinder, and the outer cylinder.
0641  *
0642  * The "chirality" of the involute is viewed from the \em +z axis looking down:
0643  * whether it spirals to the right or left.
0644  */
0645 class Involute final : public IntersectRegionInterface
0646 {
0647   public:
0648     // Construct with radius
0649     explicit Involute(Real3 const& radii,
0650                       Real2 const& displacement,
0651                       Chirality chirality,
0652                       real_type halfheight);
0653 
0654     // Build surfaces
0655     void build(IntersectSurfaceBuilder&) const final;
0656 
0657     // Output to JSON
0658     void output(JsonPimpl*) const final;
0659 
0660     //// ACCESSORS ////
0661 
0662     //! Radii: Rdius of involute, minimum radius, maximum radius
0663     Real3 const& radii() const { return radii_; }
0664     //! Displacement angle
0665     Real2 const& displacement_angle() const { return a_; }
0666     //!  Angular bounds of involute
0667     Real2 const& t_bounds() const { return t_bounds_; }
0668     //! Chirality of involute: turning left or right
0669     Chirality chirality() const { return sign_; }
0670     //! Halfheight
0671     real_type halfheight() const { return hh_; }
0672 
0673   private:
0674     Real3 radii_;
0675     Real2 a_;
0676     Real2 t_bounds_;
0677     Chirality sign_;
0678     real_type hh_;
0679 };
0680 
0681 //---------------------------------------------------------------------------//
0682 /*!
0683  * A finite *z*-aligned parabolid.
0684  *
0685  * The paraboloid is defined in an analogous fashion to the cone. A half-height
0686  * (hh) defines the z-extents, such that the centroid of the outer bounding box
0687  * is the origin. The lower and upper radii correspond to the radii at
0688  * \f$ z = \pm h \f$. Either the lower or upper radii may be 0, i.e.,
0689  * the solid may include the vertex. Degenerate cases where the lower and upper
0690  * radii are equal are not permitted: a cylinder should be used instead.
0691  *
0692  * A paraboloid with these properties is expressed in SimpleQuadric form as:
0693  * \f[
0694     x^2 + y^2 + \frac{(R_{\mathrm{lo}}^2 - R_{\mathrm{hi}}^2)}{h} z
0695     - \frac{R_{\mathrm{lo}}^2 + R_{\mathrm{hi}}^2}{2} = 0,
0696  * \f]
0697  * where \f$R_{\mathrm{lo}}\f$ and \f$R_\mathrm{hi}\f$ correspond to the lower
0698  * and upper radii, respectively, and \f$h\f$ is the full height. Note that the
0699  * scaling is such that as \f$ R_{\mathrm{lo}} \to R_{\mathrm{hi}} \f$ this
0700  * approaches the cylindrical equation \f$ x^2 + y^2 = R^2 \f$.
0701  *
0702  */
0703 class Paraboloid final : public IntersectRegionInterface
0704 {
0705   public:
0706     // Construct with lower/upper radii and the half-height
0707     Paraboloid(real_type lower_radius,
0708                real_type upper_radius,
0709                real_type halfheight);
0710 
0711     // Build surfaces
0712     void build(IntersectSurfaceBuilder&) const final;
0713 
0714     // Output to JSON
0715     void output(JsonPimpl*) const final;
0716 
0717     //// TEMPLATE INTERFACE ////
0718 
0719     // Whether this encloses another paraboloid
0720     bool encloses(Paraboloid const& other) const;
0721 
0722     //// ACCESSORS ////
0723 
0724     //! Radius at z=-hh
0725     real_type lower_radius() const { return r_lo_; }
0726 
0727     //! Radius at z=hh
0728     real_type upper_radius() const { return r_hi_; }
0729 
0730     //! Half-height along Z
0731     real_type halfheight() const { return hh_; }
0732 
0733   private:
0734     real_type r_lo_;
0735     real_type r_hi_;
0736     real_type hh_;
0737 };
0738 
0739 //---------------------------------------------------------------------------//
0740 /*!
0741  * A general parallelepiped centered on the origin.
0742  *
0743  * A parallelepiped is a shape having 3 pairs of parallel faces out of
0744  * which one is parallel with the \em x-y plane (\em z faces). All faces are
0745  * parallelograms in the general case. The \em z faces have 2 edges parallel
0746  * with the \em x axis. Note that all angle parameters are expressed in terms
0747  * of fractions of a 360-degree turn.
0748  *
0749  * The shape has the center in the origin and it is defined by:
0750  *
0751  *   - \c halfedges: a 3-vector (dY, dY, dZ) with half-lengths of the
0752  *     projections of the edges on X, Y, Z. The lower Z face is positioned at
0753  *     `-dZ`, and the upper one at `+dZ`.
0754  *   - \c alpha angle between the segment defined by the centers of the
0755  *     X-parallel edges and Y axis. Validity range is `(-1/4, 1/4)`;
0756  *   - \c theta polar angle of the shape's main axis, e.g. the segment defined
0757  *     by the centers of the Z faces. Validity range is `[0, 1/4)`;
0758  *   - \c phi azimuthal angle of the shape's main axis (as explained above).
0759  *     Validity range is `[0, 1)`.
0760  */
0761 class Parallelepiped final : public IntersectRegionInterface
0762 {
0763   public:
0764     // Construct with half widths and 3 angles
0765     Parallelepiped(Real3 const& halfedges, Turn alpha, Turn theta, Turn phi);
0766 
0767     // Build surfaces
0768     void build(IntersectSurfaceBuilder&) const final;
0769 
0770     // Output to JSON
0771     void output(JsonPimpl*) const final;
0772 
0773     //// ACCESSORS ////
0774 
0775     //! Half-lengths of edge projections along each axis
0776     Real3 const& halfedges() const { return hpr_; }
0777     //! Angle between slanted *y* edges and the *y* axis (in turns)
0778     Turn alpha() const { return alpha_; }
0779     //! Polar angle of main axis (in turns)
0780     Turn theta() const { return theta_; }
0781     //! Azimuthal angle of main axis (in turns)
0782     Turn phi() const { return phi_; }
0783 
0784   private:
0785     // half-lengths
0786     Real3 hpr_;
0787     // angle between slanted y-edges and the y-axis
0788     Turn alpha_;
0789     // polar angle of main axis
0790     Turn theta_;
0791     // azimuthal angle of main axis
0792     Turn phi_;
0793 };
0794 
0795 //---------------------------------------------------------------------------//
0796 /*!
0797  * A regular, z-extruded polygon centered on the origin.
0798  *
0799  * This is the base component of a G4Polyhedra (PGON). The default rotation is
0800  * to put the first point at \f$ y = 0 \f$. Looking at an x-y
0801  * slice for a prism with apothem \em a, odd-numbered prisms have a plane at
0802  * \f$ x=-a \f$:
0803  * - \f$ n=3 \f$ is a triangle pointing right,
0804  * - \f$ n=4 \f$ is a diamond,
0805  * - \f$ n=5 \f$ is a tilted pentagon (flat on the left), and
0806  * - \f$ n=6 \f$ is a flat-top hexagon.
0807  *
0808  *
0809  * The "orientation" parameter is a scaled counterclockwise rotation on
0810  * \f$[0, 1)\f$, where zero preserves the orientation described above, and
0811  * unity would replicate the original shape but with the "p0" face being where
0812  * the "p1" originally was. With a value of 0.5:
0813  * - \f$ n=3 \f$ is a triangle pointing left
0814  * - \f$ n=4 \f$ is an upright square,
0815  * - \f$ n=5 \f$ is a tilted pentagon (flat on the right), and
0816  * - \f$ n=6 \f$ is a pointy-top hexagon.
0817  */
0818 class Prism final : public IntersectRegionInterface
0819 {
0820   public:
0821     // Construct with inner radius (apothem), half height, and orientation
0822     Prism(int num_sides,
0823           real_type apothem,
0824           real_type halfheight,
0825           real_type orientation);
0826 
0827     // Build surfaces
0828     void build(IntersectSurfaceBuilder&) const final;
0829 
0830     // Output to JSON
0831     void output(JsonPimpl*) const final;
0832 
0833     //// TEMPLATE INTERFACE ////
0834 
0835     // Whether this encloses another cylinder
0836     bool encloses(Prism const& other) const;
0837 
0838     //// ACCESSORS ////
0839 
0840     //! Number of sides
0841     int num_sides() const { return num_sides_; }
0842     //! Inner radius
0843     real_type apothem() const { return apothem_; }
0844     //! Half the Z height
0845     real_type halfheight() const { return hh_; }
0846     //! Rotation factor
0847     real_type orientation() const { return orientation_; }
0848 
0849   private:
0850     // Number of sides
0851     int num_sides_;
0852 
0853     // Distance from center to midpoint of its side
0854     real_type apothem_;
0855 
0856     // Half the z height
0857     real_type hh_;
0858 
0859     // Rotational offset: 0 has point at (r, 0), 1 is congruent with 0
0860     real_type orientation_;
0861 };
0862 
0863 //---------------------------------------------------------------------------//
0864 /*!
0865  * A sphere centered on the origin.
0866  *
0867  * \note Be aware there's also a sphere *surface* at orange/surf/Sphere.hh in a
0868  * different namespace.
0869  */
0870 class Sphere final : public IntersectRegionInterface
0871 {
0872   public:
0873     // Construct with radius
0874     explicit Sphere(real_type radius);
0875 
0876     // Build surfaces
0877     void build(IntersectSurfaceBuilder&) const final;
0878 
0879     // Output to JSON
0880     void output(JsonPimpl*) const final;
0881 
0882     //// TEMPLATE INTERFACE ////
0883 
0884     // Whether this encloses another sphere
0885     bool encloses(Sphere const& other) const;
0886 
0887     //// ACCESSORS ////
0888 
0889     //! Radius
0890     real_type radius() const { return radius_; }
0891 
0892   private:
0893     real_type radius_;
0894 };
0895 
0896 //---------------------------------------------------------------------------//
0897 }  // namespace orangeinp
0898 }  // namespace celeritas