ImpactX
Semigaussian.H
Go to the documentation of this file.
1 /* Copyright 2022 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_SEMIGAUSSIAN
11 #define IMPACTX_DISTRIBUTION_SEMIGAUSSIAN
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 Semigaussian
24  {
35  Semigaussian(amrex::ParticleReal const sigx, amrex::ParticleReal const sigy,
36  amrex::ParticleReal const sigt,amrex::ParticleReal const sigpx,
37  amrex::ParticleReal const sigpy,amrex::ParticleReal const sigpt,
38  amrex::ParticleReal const muxpx=0.0, amrex::ParticleReal const muypy=0.0,
39  amrex::ParticleReal const mutpt=0.0
40  )
41  : m_sigmaX(sigx),m_sigmaY(sigy),m_sigmaT(sigt),m_sigmaPx(sigpx),m_sigmaPy(sigpy),
42  m_sigmaPt(sigpt),m_muxpx(muxpx),m_muypy(muypy),m_mutpt(mutpt)
43  {
44  }
45 
58  amrex::ParticleReal & x,
59  amrex::ParticleReal & y,
60  amrex::ParticleReal & t,
61  amrex::ParticleReal & px,
62  amrex::ParticleReal & py,
63  amrex::ParticleReal & pt,
64  amrex::RandomEngine const& engine) const {
65 
66  using namespace amrex::literals;
67 
68  amrex::ParticleReal ln1,u1,u2;
69  amrex::ParticleReal root,a1,a2;
70  amrex::ParticleReal phi,v,r;
71 
72  constexpr amrex::ParticleReal pi = 3.14159265358979_prt;
73 
74  // Generate a 3D uniform distribution (x,y,t) within a cylinder:
75 
76  phi = amrex::Random(engine);
77  phi = 2_prt*pi*phi;
78  v = amrex::Random(engine);
79  r = sqrt(v);
80  x = r*cos(phi);
81  y = r*sin(phi);
82  t = amrex::Random(engine);
83  t = 2_prt*(t-0.5_prt);
84 
85  // Scale to produce the identity covariance matrix:
86  amrex::ParticleReal c = sqrt(3.0_prt);
87  x = 2_prt*x;
88  y = 2_prt*y;
89  t = c*t;
90 
91  // Generate three normal random variables (px,py,pt) using Box-Muller:
92 
93  u1 = amrex::Random(engine);
94  u2 = amrex::Random(engine);
95  ln1 = sqrt(-2_prt*log(u1));
96  px = ln1*cos(2_prt*pi*u2);
97  py = ln1*sin(2_prt*pi*u2);
98 
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  // Transform to produce the desired second moments/correlations:
105  root = sqrt(1.0_prt-m_muxpx*m_muxpx);
106  a1 = m_sigmaX*x/root;
107  a2 = m_sigmaPx*(-m_muxpx*x/root+px);
108  x = a1;
109  px = a2;
110  root = sqrt(1.0_prt-m_muypy*m_muypy);
111  a1 = m_sigmaY*y/root;
112  a2 = m_sigmaPy*(-m_muypy*y/root+py);
113  y = a1;
114  py = a2;
115  root = sqrt(1.0_prt-m_mutpt*m_mutpt);
116  a1 = m_sigmaT*t/root;
117  a2 = m_sigmaPt*(-m_mutpt*t/root+pt);
118  t = a1;
119  pt = a2;
120  }
121  private:
122  amrex::ParticleReal m_sigmaX,m_sigmaY,m_sigmaT;
123  amrex::ParticleReal m_sigmaPx,m_sigmaPy,m_sigmaPt;
124  amrex::ParticleReal m_muxpx,m_muypy,m_mutpt;
125  };
126 
127 } // namespace distribution
128 } // namespace impactx
129 
130 #endif // IMPACTX_DISTRIBUTION_SEMIGAUSSIAN
Semigaussian(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: Semigaussian.H:35
amrex::ParticleReal m_mutpt
Definition: Semigaussian.H:124
amrex::ParticleReal m_sigmaPy
Definition: Semigaussian.H:123
Real Random()
Definition: ImpactX.cpp:31
amrex::ParticleReal m_sigmaY
Definition: Semigaussian.H:122
c
#define AMREX_FORCE_INLINE
amrex::ParticleReal m_muypy
Definition: Semigaussian.H:124
#define AMREX_GPU_HOST_DEVICE
amrex::ParticleReal m_sigmaT
Definition: Semigaussian.H:122
amrex::ParticleReal m_muxpx
RMS momentum.
Definition: Semigaussian.H:124
amrex::ParticleReal m_sigmaPx
related RMS sizes (length)
Definition: Semigaussian.H:123
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > log(const GpuComplex< T > &a_z) noexcept
Definition: Semigaussian.H:23
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: Semigaussian.H:57
amrex::ParticleReal m_sigmaX
Definition: Semigaussian.H:122
amrex::ParticleReal m_sigmaPt
Definition: Semigaussian.H:123
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > sqrt(const GpuComplex< T > &a_z) noexcept