ImpactX
Kurth6D.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_KURTH6D
11 #define IMPACTX_DISTRIBUTION_KURTH6D
12 
13 #include <AMReX_Random.H>
14 #include <AMReX_REAL.H>
15 
16 #include <cmath>
17 
18 
19 namespace impactx
20 {
21 namespace distribution
22 {
23  struct Kurth6D
24  {
37  Kurth6D(amrex::ParticleReal const sigx, amrex::ParticleReal const sigy,
38  amrex::ParticleReal const sigt,amrex::ParticleReal const sigpx,
39  amrex::ParticleReal const sigpy,amrex::ParticleReal const sigpt,
40  amrex::ParticleReal const muxpx=0.0, amrex::ParticleReal const muypy=0.0,
41  amrex::ParticleReal const mutpt=0.0
42  )
43  : m_sigmaX(sigx),m_sigmaY(sigy),m_sigmaT(sigt),m_sigmaPx(sigpx),m_sigmaPy(sigpy),
44  m_sigmaPt(sigpt),m_muxpx(muxpx),m_muypy(muypy),m_mutpt(mutpt)
45  {
46  }
47 
60  amrex::ParticleReal & x,
61  amrex::ParticleReal & y,
62  amrex::ParticleReal & t,
63  amrex::ParticleReal & px,
64  amrex::ParticleReal & py,
65  amrex::ParticleReal & pt,
66  amrex::RandomEngine const& engine) const {
67 
68  using namespace amrex::literals;
69 
70  amrex::ParticleReal v,costheta,sintheta,phi,r;
71  amrex::ParticleReal L,alpha,pmax,pr,beta,p1,p2;
72  amrex::ParticleReal root,a1,a2;
73 
74  constexpr amrex::ParticleReal pi = 3.14159265358979_prt;
75 
76  // Random samples used to define (x,y,z):
77  v = amrex::Random(engine);
78  costheta = amrex::Random(engine);
79  costheta = 2_prt*(costheta-0.5_prt);
80  sintheta = sqrt(1_prt-pow(costheta,2));
81  phi = amrex::Random(engine);
82  phi = 2_prt*pi*phi;
83 
84  // Transformations for (x,y,t):
85  r = pow(v,1_prt/3_prt);
86  x = r*sintheta*cos(phi);
87  y = r*sintheta*sin(phi);
88  t = r*costheta;
89 
90  // Random samples used to define L:
91  L = amrex::Random(engine);
92  L = r*sqrt(L);
93 
94  // Random samples used to define pr:
95  alpha = amrex::Random(engine);
96  alpha = pi*alpha;
97  pmax = 1_prt - pow(L/r,2) - pow(r,2) + pow(L,2);
98  pmax = sqrt(pmax);
99  pr = pmax*cos(alpha);
100 
101  // Random samples used to define ptangent:
102  beta = amrex::Random(engine);
103  beta = 2_prt*pi*beta;
104  p1 = L/r*cos(beta); // This is phi component
105  p2 = L/r*sin(beta); // This is theta component
106 
107  // Transformation from spherical to Cartesian coord.:
108  px = pr*sintheta*cos(phi) + p2*costheta*cos(phi) - p1*sin(phi);
109  py = pr*sintheta*sin(phi) + p2*costheta*sin(phi) + p1*cos(phi);
110  pt = pr*costheta - p2*sintheta;
111 
112  // Scale to produce the identity covariance matrix:
113  amrex::ParticleReal c = sqrt(5.0_prt);
114  x = c*x;
115  y = c*y;
116  t = c*t;
117  px = c*px;
118  py = c*py;
119  pt = c*pt;
120 
121  // Transform to produce the desired second moments/correlations:
122  root = sqrt(1.0_prt-m_muxpx*m_muxpx);
123  a1 = m_sigmaX*x/root;
124  a2 = m_sigmaPx*(-m_muxpx*x/root+px);
125  x = a1;
126  px = a2;
127  root = sqrt(1.0_prt-m_muypy*m_muypy);
128  a1 = m_sigmaY*y/root;
129  a2 = m_sigmaPy*(-m_muypy*y/root+py);
130  y = a1;
131  py = a2;
132  root = sqrt(1.0_prt-m_mutpt*m_mutpt);
133  a1 = m_sigmaT*t/root;
134  a2 = m_sigmaPt*(-m_mutpt*t/root+pt);
135  t = a1;
136  pt = a2;
137  }
138  private:
139  amrex::ParticleReal m_sigmaX,m_sigmaY,m_sigmaT;
140  amrex::ParticleReal m_sigmaPx,m_sigmaPy,m_sigmaPt;
141  amrex::ParticleReal m_muxpx,m_muypy,m_mutpt;
142  };
143 
144 } // namespace distribution
145 } // namespace impactx
146 
147 #endif // IMPACTX_DISTRIBUTION_KURTH6D
amrex::ParticleReal m_muypy
Definition: Kurth6D.H:141
constexpr std::enable_if_t< std::is_floating_point< T >::value, T > pi()
amrex::ParticleReal m_muxpx
RMS momentum.
Definition: Kurth6D.H:141
Definition: ImpactX.cpp:31
amrex::ParticleReal m_sigmaT
Definition: Kurth6D.H:139
Real Random()
c
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: Kurth6D.H:59
amrex::ParticleReal m_sigmaY
Definition: Kurth6D.H:139
#define AMREX_FORCE_INLINE
amrex::ParticleReal m_sigmaPt
Definition: Kurth6D.H:140
amrex::ParticleReal m_sigmaX
Definition: Kurth6D.H:139
amrex::ParticleReal m_mutpt
Definition: Kurth6D.H:141
#define AMREX_GPU_HOST_DEVICE
Definition: Kurth6D.H:23
Kurth6D(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: Kurth6D.H:37
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > pow(const GpuComplex< T > &a_z, const T &a_y) noexcept
amrex::ParticleReal m_sigmaPx
related RMS sizes (length)
Definition: Kurth6D.H:140
amrex::ParticleReal m_sigmaPy
Definition: Kurth6D.H:140
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > sqrt(const GpuComplex< T > &a_z) noexcept