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 
13 #include <AMReX_Random.H>
14 #include <AMReX_REAL.H>
15 
16 #include <cmath>
17 
18 
19 namespace impactx::distribution
20 {
21  struct Kurth4D
22  {
33  Kurth4D(amrex::ParticleReal const sigx, amrex::ParticleReal const sigy,
34  amrex::ParticleReal const sigt,amrex::ParticleReal const sigpx,
35  amrex::ParticleReal const sigpy,amrex::ParticleReal const sigpt,
36  amrex::ParticleReal const muxpx=0.0, amrex::ParticleReal const muypy=0.0,
37  amrex::ParticleReal const mutpt=0.0
38  )
39  : m_sigmaX(sigx),m_sigmaY(sigy),m_sigmaT(sigt),m_sigmaPx(sigpx),m_sigmaPy(sigpy),
40  m_sigmaPt(sigpt),m_muxpx(muxpx),m_muypy(muypy),m_mutpt(mutpt)
41  {
42  }
43 
56  amrex::ParticleReal & x,
57  amrex::ParticleReal & y,
58  amrex::ParticleReal & t,
59  amrex::ParticleReal & px,
60  amrex::ParticleReal & py,
61  amrex::ParticleReal & pt,
62  amrex::RandomEngine const& engine) const {
63 
64  using namespace amrex::literals;
65 
66  amrex::ParticleReal v,phi,r,u1,u2,ln1;
67  amrex::ParticleReal alpha,u,Lz,pmax,pr,pphi;
68  amrex::ParticleReal root,a1,a2;
69 
70  constexpr amrex::ParticleReal pi = 3.14159265358979_prt;
71 
72  // Sample and transform to define (x,y):
73  v = amrex::Random(engine);
74  phi = amrex::Random(engine);
75  phi = 2_prt*pi*phi;
76  r = sqrt(v);
77  x = r*cos(phi);
78  y = r*sin(phi);
79 
80  // Random samples used to define Lz:
81  u = amrex::Random(engine);
82  Lz = r*(2.0_prt*u-1.0_prt);
83 
84  // Random samples used to define pr:
85  alpha = amrex::Random(engine);
86  alpha = pi*alpha;
87  pmax = 1.0_prt - pow((Lz/r),2) - pow(r,2) + pow(Lz,2);
88  pmax = sqrt(pmax);
89  pr = pmax*cos(alpha);
90  pphi = Lz/r;
91 
92  // Transformations used to obtain (px,py):
93  px = pr*cos(phi)-pphi*sin(phi);
94  py = pr*sin(phi)+pphi*cos(phi);
95 
96  // Sample and transform to define (t,pt):
97  t = amrex::Random(engine);
98  t = 2.0_prt*(t-0.5_prt);
99  u1 = amrex::Random(engine);
100  u2 = amrex::Random(engine);
101  ln1 = sqrt(-2_prt*log(u1));
102  pt = ln1*cos(2_prt*pi*u2);
103 
104  // Scale to produce the identity covariance matrix:
105  amrex::ParticleReal const c = sqrt(3.0_prt);
106  x = 2_prt*x;
107  y = 2_prt*y;
108  t = c*t;
109  px = 2_prt*px;
110  py = 2_prt*py;
111  // pt = pt;
112 
113  // Transform to produce the desired second moments/correlations:
114  root = sqrt(1.0_prt-m_muxpx*m_muxpx);
115  a1 = m_sigmaX*x/root;
116  a2 = m_sigmaPx*(-m_muxpx*x/root+px);
117  x = a1;
118  px = a2;
119  root = sqrt(1.0_prt-m_muypy*m_muypy);
120  a1 = m_sigmaY*y/root;
121  a2 = m_sigmaPy*(-m_muypy*y/root+py);
122  y = a1;
123  py = a2;
124  root = sqrt(1.0_prt-m_mutpt*m_mutpt);
125  a1 = m_sigmaT*t/root;
126  a2 = m_sigmaPt*(-m_mutpt*t/root+pt);
127  t = a1;
128  pt = a2;
129  }
130  private:
131  amrex::ParticleReal m_sigmaX,m_sigmaY,m_sigmaT;
132  amrex::ParticleReal m_sigmaPx,m_sigmaPy,m_sigmaPt;
133  amrex::ParticleReal m_muxpx,m_muypy,m_mutpt;
134  };
135 
136 } // namespace impactx::distribution
137 
138 #endif // IMPACTX_DISTRIBUTION_KURTH4D
#define AMREX_FORCE_INLINE
#define AMREX_GPU_HOST_DEVICE
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:26
c
Definition: Kurth4D.H:22
Kurth4D(amrex::ParticleReal const sigx, amrex::ParticleReal const sigy, amrex::ParticleReal const sigt, amrex::ParticleReal const sigpx, amrex::ParticleReal const sigpy, amrex::ParticleReal const sigpt, amrex::ParticleReal const muxpx=0.0, amrex::ParticleReal const muypy=0.0, amrex::ParticleReal const mutpt=0.0)
Definition: Kurth4D.H:33
amrex::ParticleReal m_sigmaX
Definition: Kurth4D.H:131
amrex::ParticleReal m_sigmaY
Definition: Kurth4D.H:131
amrex::ParticleReal m_sigmaPx
related RMS sizes (length)
Definition: Kurth4D.H:132
amrex::ParticleReal m_sigmaPy
Definition: Kurth4D.H:132
amrex::ParticleReal m_sigmaPt
Definition: Kurth4D.H:132
amrex::ParticleReal m_muxpx
RMS momentum.
Definition: Kurth4D.H:133
amrex::ParticleReal m_sigmaT
Definition: Kurth4D.H:131
amrex::ParticleReal m_mutpt
Definition: Kurth4D.H:133
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator()(amrex::ParticleReal &x, amrex::ParticleReal &y, amrex::ParticleReal &t, amrex::ParticleReal &px, amrex::ParticleReal &py, amrex::ParticleReal &pt, amrex::RandomEngine const &engine) const
Definition: Kurth4D.H:55
amrex::ParticleReal m_muypy
Definition: Kurth4D.H:133