ImpactX
Semigaussian.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_SEMIGAUSSIAN
11 #define IMPACTX_DISTRIBUTION_SEMIGAUSSIAN
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 Semigaussian
26  {
38  amrex::ParticleReal lambdax,
39  amrex::ParticleReal lambday,
40  amrex::ParticleReal lambdat,
41  amrex::ParticleReal lambdapx,
42  amrex::ParticleReal lambdapy,
43  amrex::ParticleReal lambdapt,
44  amrex::ParticleReal muxpx=0.0,
45  amrex::ParticleReal muypy=0.0,
46  amrex::ParticleReal mutpt=0.0
47  )
48  : m_lambdaX(lambdax), m_lambdaY(lambday), m_lambdaT(lambdat), m_lambdaPx(lambdapx), m_lambdaPy(lambdapy),
49  m_lambdaPt(lambdapt), 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 ln1,u1,u2;
98  amrex::ParticleReal root,a1,a2;
99  amrex::ParticleReal phi,v,r;
100 
101  // Generate a 3D uniform distribution (x,y,t) within a cylinder:
102 
103  phi = amrex::Random(engine);
104  phi = 2_prt*pi*phi;
105  v = amrex::Random(engine);
106  r = sqrt(v);
107  x = r*cos(phi);
108  y = r*sin(phi);
109  t = amrex::Random(engine);
110  t = 2_prt*(t-0.5_prt);
111 
112  // Scale to produce the identity covariance matrix:
113  amrex::ParticleReal const c = sqrt(3.0_prt);
114  x = 2_prt*x;
115  y = 2_prt*y;
116  t = c*t;
117 
118  // Generate three normal random variables (px,py,pt) using Box-Muller:
119 
120  u1 = amrex::Random(engine);
121  u2 = amrex::Random(engine);
122  ln1 = sqrt(-2_prt*log(u1));
123  px = ln1*cos(2_prt*pi*u2);
124  py = ln1*sin(2_prt*pi*u2);
125 
126  u1 = amrex::Random(engine);
127  u2 = amrex::Random(engine);
128  ln1 = sqrt(-2_prt*log(u1));
129  pt = ln1*cos(2_prt*pi*u2);
130 
131  // Transform to produce the desired second moments/correlations:
132  root = sqrt(1.0_prt-m_muxpx*m_muxpx);
133  a1 = m_lambdaX * x / root;
134  a2 = m_lambdaPx * (-m_muxpx * x / root + px);
135  x = a1;
136  px = a2;
137  root = sqrt(1.0_prt-m_muypy*m_muypy);
138  a1 = m_lambdaY * y / root;
139  a2 = m_lambdaPy * (-m_muypy * y / root + py);
140  y = a1;
141  py = a2;
142  root = sqrt(1.0_prt-m_mutpt*m_mutpt);
143  a1 = m_lambdaT * t / root;
144  a2 = m_lambdaPt * (-m_mutpt * t / root + pt);
145  t = a1;
146  pt = a2;
147  }
148 
149  private:
150  amrex::ParticleReal m_lambdaX, m_lambdaY, m_lambdaT;
151  amrex::ParticleReal m_lambdaPx, m_lambdaPy, m_lambdaPt;
152  amrex::ParticleReal m_muxpx, m_muypy, m_mutpt;
153  };
154 
155 } // namespace impactx::distribution
156 
157 #endif // IMPACTX_DISTRIBUTION_SEMIGAUSSIAN
#define AMREX_FORCE_INLINE
#define AMREX_GPU_HOST_DEVICE
static constexpr amrex::Real pi
constexpr std::enable_if_t< std::is_floating_point_v< T >, 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 > sqrt(const GpuComplex< T > &a_z) noexcept
Definition: All.H:27
@ t
fixed t as the independent variable
c
Definition: ReferenceParticle.H:30
Definition: Semigaussian.H:26
amrex::ParticleReal m_lambdaPx
related position axis intercepts (length) of the phase space ellipse
Definition: Semigaussian.H:151
amrex::ParticleReal m_mutpt
Definition: Semigaussian.H:152
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: Semigaussian.H:84
void finalize()
Definition: Semigaussian.H:69
Semigaussian(amrex::ParticleReal lambdax, amrex::ParticleReal lambday, amrex::ParticleReal lambdat, amrex::ParticleReal lambdapx, amrex::ParticleReal lambdapy, amrex::ParticleReal lambdapt, amrex::ParticleReal muxpx=0.0, amrex::ParticleReal muypy=0.0, amrex::ParticleReal mutpt=0.0)
Definition: Semigaussian.H:37
amrex::ParticleReal m_lambdaT
Definition: Semigaussian.H:150
amrex::ParticleReal m_muypy
Definition: Semigaussian.H:152
amrex::ParticleReal m_lambdaPt
Definition: Semigaussian.H:151
amrex::ParticleReal m_lambdaPy
Definition: Semigaussian.H:151
amrex::ParticleReal m_lambdaY
Definition: Semigaussian.H:150
amrex::ParticleReal m_lambdaX
Definition: Semigaussian.H:150
void initialize([[maybe_unused]] amrex::ParticleReal bunch_charge, [[maybe_unused]] RefPart const &ref)
Definition: Semigaussian.H:60
amrex::ParticleReal m_muxpx
related momentum axis intercepts of the phase space ellipse
Definition: Semigaussian.H:152