ThinDipole.H

Go to the documentation of this file.

/* Copyright 2022-2023 The Regents of the University of California, through Lawrence
 *           Berkeley National Laboratory (subject to receipt of any required
 *           approvals from the U.S. Dept. of Energy). All rights reserved.
 *
 * This file is part of ImpactX.
 *
 * Authors: Chad Mitchell, Axel Huebl
 * License: BSD-3-Clause-LBNL
 */
#ifndef IMPACTX_THINDIPOLE_H
#define IMPACTX_THINDIPOLE_H
 
#include "particles/ImpactXParticleContainer.H"
#include "mixin/alignment.H"
#include "mixin/beamoptic.H"
#include "mixin/thin.H"
#include "mixin/lineartransport.H"
#include "mixin/named.H"
#include "mixin/nofinalize.H"
 
#include <AMReX_Extension.H>
#include <AMReX_Math.H>
#include <AMReX_REAL.H>
#include <AMReX_SIMD.H>
 
#include <cmath>
#include <stdexcept>
 
 
namespace impactx::elements
{
    struct ThinDipole
    : public mixin::Named,
      public mixin::BeamOptic<ThinDipole>,
      public mixin::LinearTransport<ThinDipole>,
      public mixin::Thin,
      public mixin::Alignment,
      public mixin::NoFinalize,
      public amrex::simd::Vectorized<amrex::simd::native_simd_size_particlereal>
    {
        static constexpr auto type = "ThinDipole";
        using PType = ImpactXParticleContainer::ParticleType;

        ThinDipole (
            amrex::ParticleReal theta,
            amrex::ParticleReal rc,
            amrex::ParticleReal dx = 0,
            amrex::ParticleReal dy = 0,
            amrex::ParticleReal rotation_degree = 0,
            std::optional<std::string> name = std::nullopt
        )
        : Named(std::move(name)),
          Alignment(dx, dy, rotation_degree),
          m_theta(theta * degree2rad),
          m_rc(rc)
        {
        }

        void reverse () { m_theta = -m_theta; }

        using BeamOptic::operator();

        void compute_constants (RefPart const & refpart)
        {
            using namespace amrex::literals; // for _rt and _prt
 
            Alignment::compute_constants(refpart);
 
            // access reference particle to find relativistic beta
            m_beta = refpart.beta();
 
            // compute the effective (equivalent) arc length and curvature
            m_ds = m_theta * m_rc;
            m_kx = 1.0_prt / m_rc;
        }

        template<typename T_Real=amrex::ParticleReal, typename T_IdCpu=uint64_t>
        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        void operator() (
            T_Real const & AMREX_RESTRICT x,
            [[maybe_unused]] T_Real const & AMREX_RESTRICT y,
            T_Real & AMREX_RESTRICT t,
            T_Real & AMREX_RESTRICT px,
            [[maybe_unused]] T_Real const & AMREX_RESTRICT py,
            T_Real const & AMREX_RESTRICT pt,
            [[maybe_unused]] T_IdCpu const & AMREX_RESTRICT idcpu,
            [[maybe_unused]] RefPart const & AMREX_RESTRICT refpart
        ) const
        {
            using namespace amrex::literals; // for _rt and _prt
            using amrex::Math::powi;
            using namespace std; // for cmath(float)
 
            // access position data
            T_Real const xout = x;
            // T_Real const yout = y;
            T_Real const tout = t;
 
            // initialize output values of momenta
            T_Real pxout = px;
            // T_Real pyout = py;
            // T_Real ptout = pt;
 
            // compute the function expressing dp/p in terms of pt (labeled f in Ripken etc.)
            T_Real const f = -1.0_prt + sqrt(1.0_prt - 2.0_prt * pt / m_beta + powi<2>(pt));
            T_Real const fprime = (1.0_prt - m_beta*pt) / (m_beta * (1.0_prt + f));
 
            // advance position and momentum
            // x = xout;
            pxout = px - powi<2>(m_kx) * m_ds * xout + m_kx * m_ds * f; //eq (3.2b)
 
            // y = yout;
            // pyout = py;
 
            t = tout + m_kx * xout * m_ds * fprime; //eq (3.2e)
            // ptout = pt;
 
            // assign updated momenta
            px = pxout;
            // py = pyout;
            // pt = ptout;
        }

        using Thin::operator();

        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        void operator() (RefPart & AMREX_RESTRICT refpart) const
        {
            // assign input reference particle values
            amrex::ParticleReal const px = refpart.px;
            amrex::ParticleReal const pz = refpart.pz;
 
            // calculate expensive terms once
            auto const [sin_theta, cos_theta] = amrex::Math::sincos(m_theta);
 
            // advance position and momentum (thin dipole)
            refpart.px = px*cos_theta - pz*sin_theta;
            refpart.pz = pz*cos_theta + px*sin_theta;
        }

        using LinearTransport::operator();

        AMREX_GPU_HOST AMREX_FORCE_INLINE
        Map6x6
        transport_map (RefPart const & AMREX_RESTRICT refpart) const
        {
            using namespace amrex::literals; // for _rt and _prt
            using amrex::Math::powi;
 
            amrex::ParticleReal const beta = refpart.beta();
            amrex::ParticleReal const ds = m_theta * m_rc;     // effective arc length
            amrex::ParticleReal const kx = 1.0_prt / m_rc;     // curvature
 
            Map6x6 R = Map6x6::Identity();
            R(2,1) = -powi<2>(kx) * ds;      // geometric focusing
            R(2,6) = -kx * ds / beta;        // chromatic dispersion kick
            R(5,1) = +kx * ds / beta;        // path-length coupling to x
 
            // apply the transverse rotation (roll) alignment error
            return rotate_aligned_map(R);
        }
 
        amrex::ParticleReal m_theta;  
        amrex::ParticleReal m_rc;     
 
      private:
        // constants that are independent of the individually tracked particle,
        // see: compute_constants() to refresh
        amrex::ParticleReal m_beta;  
        amrex::ParticleReal m_ds;    
        amrex::ParticleReal m_kx;    
    };
 
} // namespace impactx
 
IMPACTX_PUSH_EXTERN_TEMPLATE(impactx::elements::ThinDipole)
 
#endif // IMPACTX_THINDIPOLE_H