celeritas/field/ZHelixStepper.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2020-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 celeritas/field/ZHelixStepper.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <type_traits>
0011
0012 #include "corecel/Types.hh"
0013 #include "corecel/math/Algorithms.hh"
0014
0015 #include "Types.hh"
0016
0017 namespace celeritas
0018 {
0019 class UniformZField;
0020 //---------------------------------------------------------------------------//
0021 /*!
0022  * Analytically step along a helical path for a uniform Z magnetic field.
0023  *
0024  * The analytical solution for the motion of a charged particle in a uniform
0025  * magnetic field along the z-direction.
0026  */
0027 template<class EquationT>
0028 class ZHelixStepper
0029 {
0030     static_assert(
0031         std::is_same<std::remove_cv_t<std::remove_reference_t<
0032                          typename std::remove_reference_t<EquationT>::Field_t>>,
0033                      UniformZField>::value,
0034         "ZHelix stepper only works with UniformZField");
0035
0036   public:
0037     //!@{
0038     //! \name Type aliases
0039     using result_type = FieldStepperResult;
0040     //!@}
0041
0042   public:
0043     //! Construct with the equation of motion
0044     explicit CELER_FUNCTION ZHelixStepper(EquationT&& eq)
0045         : calc_rhs_(::celeritas::forward<EquationT>(eq))
0046     {
0047     }
0048
0049     // Adaptive step size control
0050     CELER_FUNCTION auto
0051     operator()(real_type step, OdeState const& beg_state) const -> result_type;
0052
0053   private:
0054     //// DATA ////
0055
0056     // Evaluate the equation of the motion
0057     EquationT calc_rhs_;
0058
0059     //// HELPER TYPES ////
0060     enum class Helicity : bool
0061     {
0062         positive,
0063         negative
0064     };
0065
0066     //// HELPER FUNCTIONS ////
0067
0068     // Analytical solution for a given step along a helix trajectory
0069     CELER_FUNCTION OdeState move(real_type step,
0070                                  real_type radius,
0071                                  Helicity helicity,
0072                                  OdeState const& beg_state,
0073                                  OdeState const& rhs) const;
0074
0075     //// COMMON PROPERTIES ////
0076
0077     static CELER_CONSTEXPR_FUNCTION real_type tolerance()
0078     {
0079         if constexpr (std::is_same_v<real_type, double>)
0080             return 1e-10;
0081         else if constexpr (std::is_same_v<real_type, float>)
0082             return 1e-5f;
0083     }
0084 };
0085
0086 //---------------------------------------------------------------------------//
0087 // DEDUCTION GUIDES
0088 //---------------------------------------------------------------------------//
0089 template<class EquationT>
0090 CELER_FUNCTION ZHelixStepper(EquationT&&) -> ZHelixStepper<EquationT>;
0091
0092 //---------------------------------------------------------------------------//
0093 // INLINE DEFINITIONS
0094 //---------------------------------------------------------------------------//
0095 /*!
0096  * An explicit helix stepper with analytical solutions at the end and the
0097  * middle point for a given step.
0098  *
0099  * Assuming that the magnetic field is uniform and chosen to be parallel to
0100  * the z-axis, \f$B = (0, 0, B_z)\f$, without loss of generality, the motion of
0101  * a charged particle is described by a helix trajectory. For this algorithm,
0102  * the radius of the helix, \f$R = m gamma v/(qB)\f$ and the helicity, defined
0103  * as \f$ -sign(q B_z)\f$ are evaluated through the right hand side of the ODE
0104  * equation where q is the charge of the particle.
0105  */
0106 template<class E>
0107 CELER_FUNCTION auto
0108 ZHelixStepper<E>::operator()(real_type step,
0109                              OdeState const& beg_state) const -> result_type
0110 {
0111     result_type result;
0112
0113     // Evaluate the right hand side of the equation
0114     OdeState rhs = calc_rhs_(beg_state);
0115
0116     // Calculate the radius of the helix
0117     real_type radius = std::sqrt(dot_product(beg_state.mom, beg_state.mom)
0118                                  - ipow<2>(beg_state.mom[2]))
0119                        / norm(rhs.mom);
0120
0121     // Set the helicity: 1(-1) for negative(positive) charge with Bz > 0
0122     Helicity helicity = Helicity(rhs.mom[0] / rhs.pos[1] > 0);
0123
0124     // State after the half step
0125     result.mid_state
0126         = this->move(real_type(0.5) * step, radius, helicity, beg_state, rhs);
0127
0128     // State after the full step
0129     result.end_state = this->move(step, radius, helicity, beg_state, rhs);
0130
0131     // Solution are exact, but assign a tolerance for numerical treatments
0132     result.err_state.pos.fill(ZHelixStepper::tolerance());
0133     result.err_state.mom.fill(ZHelixStepper::tolerance());
0134
0135     return result;
0136 }
0137
0138 //---------------------------------------------------------------------------//
0139 /*!
0140  * Integration for a given step length on a helix.
0141  *
0142  * Equations of a charged particle motion in a uniform magnetic field,
0143  * \f$B(0, 0, B_z)\f$ along the curved trajectory \f$ds = v dt\f$ are
0144  * \f[
0145  *  d^2 x/ds^2 =  q/p (dy/ds) B_z
0146  *  d^2 y/ds^2 = -q/p (dx/ds) B_z
0147  *  d^2 z/ds^2 =  0
0148  * \f]
0149  * where \em q and \em p are the charge and the absolute momentum of the
0150  * particle, respectively. Since the motion in the perpendicular plane with
0151  * respected to the to the magnetic field is circular with a constant
0152  * \f$p_{\perp}\f$, the final ODE state of the perpendicular motion on the
0153  * circle for a given step length \em s is
0154  * \f[
0155  *  (x, y) = M(\phi) (x_0, y_0)^T
0156  *  (px, py) = M(\phi) (px_0, py_0)^T
0157  * \f]
0158  * where \f$\phi = s/R\f$ is the azimuth angle of the particle position between
0159  * the start and the end position and \f$M(\phi)\f$ is the rotational matrix.
0160  * The solution for the parallel direction along the field is trivial.
0161  */
0162 template<class E>
0163 CELER_FUNCTION OdeState ZHelixStepper<E>::move(real_type step,
0164                                                real_type radius,
0165                                                Helicity helicity,
0166                                                OdeState const& beg_state,
0167                                                OdeState const& rhs) const
0168 {
0169     // Solution for position and momentum after moving delta_phi on the helix
0170     real_type del_phi = (helicity == Helicity::positive) ? step / radius
0171                                                          : -step / radius;
0172     real_type sin_phi = std::sin(del_phi);
0173     real_type cos_phi = std::cos(del_phi);
0174
0175     OdeState end_state;
0176     end_state.pos = {(beg_state.pos[0] * cos_phi - beg_state.pos[1] * sin_phi),
0177                      (beg_state.pos[0] * sin_phi + beg_state.pos[1] * cos_phi),
0178                      beg_state.pos[2] + del_phi * radius * rhs.pos[2]};
0179
0180     end_state.mom = {rhs.pos[0] * cos_phi - rhs.pos[1] * sin_phi,
0181                      rhs.pos[0] * sin_phi + rhs.pos[1] * cos_phi,
0182                      rhs.pos[2]};
0183
0184     real_type momentum = norm(beg_state.mom);
0185     for (int i = 0; i < 3; ++i)
0186     {
0187         end_state.mom[i] *= momentum;
0188     }
0189
0190     return end_state;
0191 }
0192
0193 //---------------------------------------------------------------------------//
0194 }  // namespace celeritas