ImpactX
Kurth4D.H
Go to the documentation of this file.
1 /* Copyright 2022-2023 The Regents of the University of California, through Lawrence
2  * Berkeley National Laboratory (subject to receipt of any required
3  * approvals from the U.S. Dept. of Energy). All rights reserved.
4  *
5  * This file is part of ImpactX.
6  *
7  * Authors: Chad Mitchell, Axel Huebl
8  * License: BSD-3-Clause-LBNL
9  */
10 #ifndef IMPACTX_DISTRIBUTION_KURTH4D
11 #define IMPACTX_DISTRIBUTION_KURTH4D
12 
14 
15 #include <ablastr/constant.H>
16 
17 #include <AMReX_Random.H>
18 #include <AMReX_REAL.H>
19 
20 #include <cmath>
21 
22 
23 namespace impactx::distribution
24 {
25  struct Kurth4D
26  {
38  amrex::ParticleReal sigx,
39  amrex::ParticleReal sigy,
40  amrex::ParticleReal sigt,
41  amrex::ParticleReal sigpx,
42  amrex::ParticleReal sigpy,
43  amrex::ParticleReal sigpt,
44  amrex::ParticleReal muxpx=0.0,
45  amrex::ParticleReal muypy=0.0,
46  amrex::ParticleReal mutpt=0.0
47  )
48  : m_sigmaX(sigx), m_sigmaY(sigy), m_sigmaT(sigt), m_sigmaPx(sigpx), m_sigmaPy(sigpy),
49  m_sigmaPt(sigpt), m_muxpx(muxpx), m_muypy(muypy), m_mutpt(mutpt)
50  {
51  }
52 
60  void initialize ([[maybe_unused]] amrex::ParticleReal bunch_charge, [[maybe_unused]] RefPart const & ref)
61  {
62  }
63 
68  void
70  {
71  }
72 
85  amrex::ParticleReal & AMREX_RESTRICT x,
86  amrex::ParticleReal & AMREX_RESTRICT y,
87  amrex::ParticleReal & AMREX_RESTRICT t,
88  amrex::ParticleReal & AMREX_RESTRICT px,
89  amrex::ParticleReal & AMREX_RESTRICT py,
90  amrex::ParticleReal & AMREX_RESTRICT pt,
91  amrex::RandomEngine const & engine
92  ) const
93  {
94  using namespace amrex::literals;
96 
97  amrex::ParticleReal v,phi,r,u1,u2,ln1;
98  amrex::ParticleReal alpha,u,Lz,pmax,pr,pphi;
99  amrex::ParticleReal root,a1,a2;
100 
101  // Sample and transform to define (x,y):
102  v = amrex::Random(engine);
103  phi = amrex::Random(engine);
104  phi = 2_prt*pi*phi;
105  r = sqrt(v);
106  x = r*cos(phi);
107  y = r*sin(phi);
108 
109  // Random samples used to define Lz:
110  u = amrex::Random(engine);
111  Lz = r*(2.0_prt*u-1.0_prt);
112 
113  // Random samples used to define pr:
114  alpha = amrex::Random(engine);
115  alpha = pi*alpha;
116  pmax = 1.0_prt - pow((Lz/r),2) - pow(r,2) + pow(Lz,2);
117  pmax = sqrt(pmax);
118  pr = pmax*cos(alpha);
119  pphi = Lz/r;
120 
121  // Transformations used to obtain (px,py):
122  px = pr*cos(phi)-pphi*sin(phi);
123  py = pr*sin(phi)+pphi*cos(phi);
124 
125  // Sample and transform to define (t,pt):
126  t = amrex::Random(engine);
127  t = 2.0_prt*(t-0.5_prt);
128  u1 = amrex::Random(engine);
129  u2 = amrex::Random(engine);
130  ln1 = sqrt(-2_prt*log(u1));
131  pt = ln1*cos(2_prt*pi*u2);
132 
133  // Scale to produce the identity covariance matrix:
134  amrex::ParticleReal const c = sqrt(3.0_prt);
135  x = 2_prt*x;
136  y = 2_prt*y;
137  t = c*t;
138  px = 2_prt*px;
139  py = 2_prt*py;
140  // pt = pt;
141 
142  // Transform to produce the desired second moments/correlations:
143  root = sqrt(1.0_prt-m_muxpx*m_muxpx);
144  a1 = m_sigmaX*x/root;
145  a2 = m_sigmaPx*(-m_muxpx*x/root+px);
146  x = a1;
147  px = a2;
148  root = sqrt(1.0_prt-m_muypy*m_muypy);
149  a1 = m_sigmaY*y/root;
150  a2 = m_sigmaPy*(-m_muypy*y/root+py);
151  y = a1;
152  py = a2;
153  root = sqrt(1.0_prt-m_mutpt*m_mutpt);
154  a1 = m_sigmaT*t/root;
155  a2 = m_sigmaPt*(-m_mutpt*t/root+pt);
156  t = a1;
157  pt = a2;
158  }
159 
160  private:
161  amrex::ParticleReal m_sigmaX, m_sigmaY, m_sigmaT;
162  amrex::ParticleReal m_sigmaPx, m_sigmaPy, m_sigmaPt;
163  amrex::ParticleReal m_muxpx, m_muypy, m_mutpt;
164  };
165 
166 } // namespace impactx::distribution
167 
168 #endif // IMPACTX_DISTRIBUTION_KURTH4D
#define AMREX_FORCE_INLINE
#define AMREX_GPU_HOST_DEVICE
static constexpr amrex::Real pi
constexpr std::enable_if_t< std::is_floating_point< T >::value, T > pi()
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > log(const GpuComplex< T > &a_z) noexcept
Real Random()
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > pow(const GpuComplex< T > &a_z, const T &a_y) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > sqrt(const GpuComplex< T > &a_z) noexcept
Definition: All.H:27
@ t
fixed t as the independent variable
c
Definition: ReferenceParticle.H:30
Definition: Kurth4D.H:26
amrex::ParticleReal m_sigmaX
Definition: Kurth4D.H:161
void initialize([[maybe_unused]] amrex::ParticleReal bunch_charge, [[maybe_unused]] RefPart const &ref)
Definition: Kurth4D.H:60
amrex::ParticleReal m_sigmaY
Definition: Kurth4D.H:161
amrex::ParticleReal m_sigmaPx
related RMS sizes (length)
Definition: Kurth4D.H:162
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator()(amrex::ParticleReal &AMREX_RESTRICT x, amrex::ParticleReal &AMREX_RESTRICT y, amrex::ParticleReal &AMREX_RESTRICT t, amrex::ParticleReal &AMREX_RESTRICT px, amrex::ParticleReal &AMREX_RESTRICT py, amrex::ParticleReal &AMREX_RESTRICT pt, amrex::RandomEngine const &engine) const
Definition: Kurth4D.H:84
amrex::ParticleReal m_sigmaPy
Definition: Kurth4D.H:162
amrex::ParticleReal m_sigmaPt
Definition: Kurth4D.H:162
Kurth4D(amrex::ParticleReal sigx, amrex::ParticleReal sigy, amrex::ParticleReal sigt, amrex::ParticleReal sigpx, amrex::ParticleReal sigpy, amrex::ParticleReal sigpt, amrex::ParticleReal muxpx=0.0, amrex::ParticleReal muypy=0.0, amrex::ParticleReal mutpt=0.0)
Definition: Kurth4D.H:37
amrex::ParticleReal m_muxpx
RMS momentum.
Definition: Kurth4D.H:163
amrex::ParticleReal m_sigmaT
Definition: Kurth4D.H:161
amrex::ParticleReal m_mutpt
Definition: Kurth4D.H:163
amrex::ParticleReal m_muypy
Definition: Kurth4D.H:163
void finalize()
Definition: Kurth4D.H:69