Cut gaussian beam (#868)

* add option to cut gaussian

* cut gaussian distrubution at injection

* add test for cut gaussian beam

* increase tolerance

* doc: mention that cutting the gaussian results in a lower total charge

* more explicit

* beware the const-able

1 files changed, 21 insertions, 13 deletions

diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index 7010da425..1b5d81f38 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -228,10 +228,12 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame(Real& x, Real& y, Real
 }
 
 void
-PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m,
-                                           Real x_rms, Real y_rms, Real z_rms,
-                                           Real q_tot, long npart,
-                                           int do_symmetrize) {
+PhysicalParticleContainer::AddGaussianBeam (
+    const Real x_m, const Real y_m, const Real z_m,
+    const Real x_rms, const Real y_rms, const Real z_rms,
+    const Real x_cut, const Real y_cut, const Real z_cut,
+    const Real q_tot, long npart,
+    const int do_symmetrize) {
 
     std::mt19937_64 mt(0451);
     std::normal_distribution<double> distx(x_m, x_rms);
@@ -256,17 +258,20 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m,
         }
         for (long i = 0; i < npart; ++i) {
 #if (defined WARPX_DIM_3D) || (WARPX_DIM_RZ)
-            Real weight = q_tot/npart/charge;
-            Real x = distx(mt);
-            Real y = disty(mt);
-            Real z = distz(mt);
+            const Real weight = q_tot/(npart*charge);
+            const Real x = distx(mt);
+            const Real y = disty(mt);
+            const Real z = distz(mt);
 #elif (defined WARPX_DIM_XZ)
-            Real weight = q_tot/npart/charge/y_rms;
-            Real x = distx(mt);
-            Real y = 0.;
-            Real z = distz(mt);
+            const Real weight = q_tot/(npart*charge*y_rms);
+            const Real x = distx(mt);
+            constexpr Real y = 0.;
+            const Real z = distz(mt);
 #endif
-            if (plasma_injector->insideBounds(x, y, z)) {
+            if (plasma_injector->insideBounds(x, y, z)  &&
+                std::abs( x - x_m ) < x_cut * x_rms     &&
+                std::abs( y - y_m ) < y_cut * y_rms     &&
+                std::abs( z - z_m ) < z_cut * z_rms   ) {
                 XDim3 u = plasma_injector->getMomentum(x, y, z);
                 u.x *= PhysConst::c;
                 u.y *= PhysConst::c;
@@ -388,6 +393,9 @@ PhysicalParticleContainer::AddParticles (int lev)
                         plasma_injector->x_rms,
                         plasma_injector->y_rms,
                         plasma_injector->z_rms,
+                        plasma_injector->x_cut,
+                        plasma_injector->y_cut,
+                        plasma_injector->z_cut,
                         plasma_injector->q_tot,
                         plasma_injector->npart,
                         plasma_injector->do_symmetrize);