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