From 19486296becb22f1e848a60718d9a59d49431593 Mon Sep 17 00:00:00 2001
From: Chad Mitchell <ChadMitchell@lbl.gov>
Date: Fri, 11 Oct 2024 11:12:39 -0700
Subject: [PATCH 01/14] Add periodic masking option to Aperture.

---
 docs/source/usage/parameters.rst   |  2 ++
 docs/source/usage/python.rst       |  2 ++
 src/initialization/InitElement.cpp |  6 +++++-
 src/particles/elements/Aperture.H  | 21 ++++++++++++++++++---
 src/python/elements.cpp            | 16 +++++++++++++++-
 5 files changed, 42 insertions(+), 5 deletions(-)

diff --git a/docs/source/usage/parameters.rst b/docs/source/usage/parameters.rst
index 1b0f2644b..9d4c63620 100644
--- a/docs/source/usage/parameters.rst
+++ b/docs/source/usage/parameters.rst
@@ -401,6 +401,8 @@ Lattice Elements
 
             * ``<element_name>.xmax`` (``float``, in meters) maximum value of the horizontal coordinate
             * ``<element_name>.ymax`` (``float``, in meters) maximum value of the vertical coordinate
+            * ``<element_name>.repeat_x`` (``float``, in meters) horizontal period for repeated aperture masking (inactive by default)
+            * ``<element_name>.repeat_y`` (``float``, in meters) vertical period for repeated aperture masking (inactive by default)
             * ``<element_name>.shape`` (``string``) shape of the aperture boundary: ``rectangular`` (default) or ``elliptical``
             * ``<element_name>.dx`` (``float``, in meters) horizontal translation error
             * ``<element_name>.dy`` (``float``, in meters) vertical translation error
diff --git a/docs/source/usage/python.rst b/docs/source/usage/python.rst
index 73e9c0329..2c91d557f 100644
--- a/docs/source/usage/python.rst
+++ b/docs/source/usage/python.rst
@@ -936,6 +936,8 @@ This module provides elements for the accelerator lattice.
 
    :param xmax: maximum allowed value of the horizontal coordinate (meter)
    :param ymax: maximum allowed value of the vertical coordinate (meter)
+   :param repeat_x: horizontal period for repeated aperture masking (inactive by default) (meter)
+   :param repeat_y: vertical period for repeated aperture masking (inactive by default) (meter)
    :param shape: aperture boundary shape: ``"rectangular"`` (default) or ``"elliptical"``
    :param dx: horizontal translation error in m
    :param dy: vertical translation error in m
diff --git a/src/initialization/InitElement.cpp b/src/initialization/InitElement.cpp
index 71ec7e4aa..d666322b5 100644
--- a/src/initialization/InitElement.cpp
+++ b/src/initialization/InitElement.cpp
@@ -376,9 +376,13 @@ namespace detail
             auto a = detail::query_alignment(pp_element);
 
             amrex::Real xmax, ymax;
+            amrex::ParticleReal repeat_x = 0.0;
+            amrex::ParticleReal repeat_y = 0.0;
             std::string shape_str = "rectangular";
             pp_element.get("xmax", xmax);
             pp_element.get("ymax", ymax);
+            pp_element.queryAdd("repeat_x", repeat_x);
+            pp_element.queryAdd("repeat_y", repeat_y);
             pp_element.queryAdd("shape", shape_str);
             AMREX_ALWAYS_ASSERT_WITH_MESSAGE(shape_str == "rectangular" || shape_str == "elliptical",
                                              element_name + ".shape must be \"rectangular\" or \"elliptical\"");
@@ -386,7 +390,7 @@ namespace detail
                                         Aperture::Shape::rectangular :
                                         Aperture::Shape::elliptical;
 
-            m_lattice.emplace_back( Aperture(xmax, ymax, shape, a["dx"], a["dy"], a["rotation_degree"], element_name) );
+            m_lattice.emplace_back( Aperture(xmax, ymax, repeat_x, repeat_y, shape, a["dx"], a["dy"], a["rotation_degree"], element_name) );
         } else if (element_type == "beam_monitor")
         {
             std::string openpmd_name = element_name;
diff --git a/src/particles/elements/Aperture.H b/src/particles/elements/Aperture.H
index 9c286870e..8d185ca03 100644
--- a/src/particles/elements/Aperture.H
+++ b/src/particles/elements/Aperture.H
@@ -51,6 +51,8 @@ namespace impactx
          * @param shape aperture shape
          * @param xmax maximum value of horizontal coordinate (m)
          * @param ymax maximum value of vertical coordinate (m)
+         * @param repeat_x horizontal period for repeated masking, optional (m)
+         * @param repeat_y vertical period for repeated masking, optional (m)
          * @param dx horizontal translation error in m
          * @param dy vertical translation error in m
          * @param rotation_degree rotation error in the transverse plane [degrees]
@@ -59,6 +61,8 @@ namespace impactx
         Aperture (
             amrex::ParticleReal xmax,
             amrex::ParticleReal ymax,
+            amrex::ParticleReal repeat_x,
+            amrex::ParticleReal repeat_y,
             Shape shape,
             amrex::ParticleReal dx = 0,
             amrex::ParticleReal dy = 0,
@@ -67,7 +71,7 @@ namespace impactx
         )
         : Named(name),
           Alignment(dx, dy, rotation_degree),
-          m_shape(shape), m_xmax(xmax), m_ymax(ymax)
+          m_shape(shape), m_xmax(xmax), m_ymax(ymax), m_repeat_x(repeat_x), m_repeat_y(repeat_y)
         {
         }
 
@@ -104,8 +108,17 @@ namespace impactx
             shift_in(x, y, px, py);
 
             // scale horizontal and vertical coordinates
-            amrex::ParticleReal const u = x / m_xmax;
-            amrex::ParticleReal const v = y / m_ymax;
+            amrex::ParticleReal u;
+            amrex::ParticleReal v;
+            // first case: a single aperture
+            if (m_repeat_x == 0.0 || m_repeat_y == 0.0) {
+               u = x / m_xmax;
+               v = y / m_ymax;
+            // second case: a periodically repeated aperture
+            } else {
+               u = std::fmod(x,m_repeat_x) / m_xmax;
+               v = std::fmod(y,m_repeat_y) / m_ymax;
+            }
 
             // compare against the aperture boundary
             switch (m_shape)
@@ -133,6 +146,8 @@ namespace impactx
         Shape m_shape; //! aperture type (rectangular, elliptical)
         amrex::ParticleReal m_xmax; //! maximum horizontal coordinate
         amrex::ParticleReal m_ymax; //! maximum vertical coordinate
+        amrex::ParticleReal m_repeat_x; //! horizontal period for repeated masking
+        amrex::ParticleReal m_repeat_y; //! vertical period for repeated masking
 
     };
 
diff --git a/src/python/elements.cpp b/src/python/elements.cpp
index b6403e71b..4d162573a 100644
--- a/src/python/elements.cpp
+++ b/src/python/elements.cpp
@@ -286,6 +286,8 @@ void init_elements(py::module& m)
         .def(py::init([](
                  amrex::ParticleReal xmax,
                  amrex::ParticleReal ymax,
+                 amrex::ParticleReal repeat_x,
+                 amrex::ParticleReal repeat_y,
                  std::string const & shape,
                  amrex::ParticleReal dx,
                  amrex::ParticleReal dy,
@@ -299,10 +301,12 @@ void init_elements(py::module& m)
                  Aperture::Shape const s = shape == "rectangular" ?
                      Aperture::Shape::rectangular :
                      Aperture::Shape::elliptical;
-                 return new Aperture(xmax, ymax, s, dx, dy, rotation_degree, name);
+                 return new Aperture(xmax, ymax, repeat_x, repeat_y, s, dx, dy, rotation_degree, name);
              }),
              py::arg("xmax"),
              py::arg("ymax"),
+             py::arg("repeat_x") = 0,
+             py::arg("repeat_y") = 0,
              py::arg("shape") = "rectangular",
              py::arg("dx") = 0,
              py::arg("dy") = 0,
@@ -336,6 +340,16 @@ void init_elements(py::module& m)
             [](Aperture & ap, amrex::ParticleReal ymax) { ap.m_ymax = ymax; },
             "maximum vertical coordinate"
         )
+        .def_property("repeat_x",
+            [](Aperture & ap) { return ap.m_repeat_x; },
+            [](Aperture & ap, amrex::ParticleReal repeat_x) { ap.m_repeat_x = repeat_x; },
+            "horizontal period for repeated aperture masking"
+        )
+        .def_property("repeat_y",
+            [](Aperture & ap) { return ap.m_repeat_y; },
+            [](Aperture & ap, amrex::ParticleReal repeat_y) { ap.m_repeat_y = repeat_y; },
+            "vertical period for repeated aperture masking"
+        )
     ;
     register_beamoptics_push(py_Aperture);
 

From fea841c9e4247ac27af8acf2bb60d1a628b4e95a Mon Sep 17 00:00:00 2001
From: Chad Mitchell <ChadMitchell@lbl.gov>
Date: Fri, 11 Oct 2024 14:44:22 -0700
Subject: [PATCH 02/14] Add draft files for benchmark example.

---
 examples/aperture/README.rst                  |  52 +++++++++
 .../aperture/analysis_aperture_pepperpot.py   | 108 ++++++++++++++++++
 examples/aperture/input_aperture_pepperpot.in |  52 +++++++++
 examples/aperture/run_aperture_pepperpot.py   |  69 +++++++++++
 4 files changed, 281 insertions(+)
 create mode 100755 examples/aperture/analysis_aperture_pepperpot.py
 create mode 100644 examples/aperture/input_aperture_pepperpot.in
 create mode 100755 examples/aperture/run_aperture_pepperpot.py

diff --git a/examples/aperture/README.rst b/examples/aperture/README.rst
index 2592ebc38..bb21b4ced 100644
--- a/examples/aperture/README.rst
+++ b/examples/aperture/README.rst
@@ -50,3 +50,55 @@ We run the following script to analyze correctness:
    .. literalinclude:: analysis_aperture.py
       :language: python3
       :caption: You can copy this file from ``examples/aperture/analysis_aperture.py``.
+
+
+.. _examples-aperture-pepperpot:
+
+Aperture Collimation with Periodic Masking
+===========================================
+
+Proton beam undergoing collimation by a periodic array of rectangular apertures, such as those used in a pepperpot emittance measurement. 
+
+We use a 250 MeV proton beam with a horizontal rms beam size of 1.56 mm and a vertical rms beam size of 2.21 mm.
+
+After a short drift of 0.123 m, the beam intercepts a periodic array of apertures of period 1 mm (in the horizontal and vertical), each of which is 0.15 mm x 0.1 mm in size.
+
+In this test, the initial values of :math:`\sigma_x`, :math:`\sigma_y`, :math:`\sigma_t`, :math:`\epsilon_x`, :math:`\epsilon_y`, and :math:`\epsilon_t` must agree with nomin>
+The test fails if:
+
+* any of the final coordinates for the valid (not lost) particles lie outside the aperture boundary or
+* any of the lost particles are inside the aperture boundary or
+* if the sum of lost and kept particles is not equal to the initial particles
+
+
+Run
+---
+
+This example can be run as a Python script (``python3 run_aperture_pepperpot.py``) or with an app with an input file (``impactx input_aperture_pepperpot.in``).
+Each can also be prefixed with an `MPI executor <https://www.mpi-forum.org>`__, such as ``mpiexec -n 4 ...`` or ``srun -n 4 ...``, depending on the system.
+
+.. tab-set::
+
+   .. tab-item:: Python Script
+
+       .. literalinclude:: run_aperture_pepperpot.py
+          :language: python3
+          :caption: You can copy this file from ``examples/aperture/run_aperture_pepperpot.py``.
+
+   .. tab-item:: App Input File
+
+       .. literalinclude:: input_aperture_pepperpot.in
+          :language: ini
+          :caption: You can copy this file from ``examples/aperture/input_aperture_pepperpot.in``.
+
+Analyze
+-------
+
+We run the following script to analyze correctness:
+
+.. dropdown:: Script ``analysis_aperture_pepperpot.py``
+
+   .. literalinclude:: analysis_aperture_pepperpot.py
+      :language: python3
+      :caption: You can copy this file from ``examples/aperture/analysis_aperture_pepperpot.py``.
+
diff --git a/examples/aperture/analysis_aperture_pepperpot.py b/examples/aperture/analysis_aperture_pepperpot.py
new file mode 100755
index 000000000..57f094fd4
--- /dev/null
+++ b/examples/aperture/analysis_aperture_pepperpot.py
@@ -0,0 +1,108 @@
+#!/usr/bin/env python3
+#
+# Copyright 2022-2023 ImpactX contributors
+# Authors: Axel Huebl, Chad Mitchell
+# License: BSD-3-Clause-LBNL
+#
+
+import numpy as np
+import openpmd_api as io
+from scipy.stats import moment
+
+
+def get_moments(beam):
+    """Calculate standard deviations of beam position & momenta
+    and emittance values
+
+    Returns
+    -------
+    sigx, sigy, sigt, emittance_x, emittance_y, emittance_t
+    """
+    sigx = moment(beam["position_x"], moment=2) ** 0.5  # variance -> std dev.
+    sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
+    sigy = moment(beam["position_y"], moment=2) ** 0.5
+    sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
+    sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
+
+    epstrms = beam.cov(ddof=0)
+    emittance_x = (sigx**2 * sigpx**2 - epstrms["position_x"]["momentum_x"] ** 2) ** 0.5
+    emittance_y = (sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2) ** 0.5
+    emittance_t = (sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2) ** 0.5
+
+    return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
+
+
+# initial/final beam
+series = io.Series("diags/openPMD/monitor.h5", io.Access.read_only)
+last_step = list(series.iterations)[-1]
+initial = series.iterations[1].particles["beam"].to_df()
+final = series.iterations[last_step].particles["beam"].to_df()
+
+series_lost = io.Series("diags/openPMD/particles_lost.h5", io.Access.read_only)
+particles_lost = series_lost.iterations[0].particles["beam"].to_df()
+
+# compare number of particles
+num_particles = 1000000
+assert num_particles == len(initial)
+# we lost particles in apertures
+assert num_particles > len(final)
+assert num_particles == len(particles_lost) + len(final)
+
+print("Initial Beam:")
+sigx, sigy, sigt, emittance_x, emittance_y, emittance_t = get_moments(initial)
+print(f"  sigx={sigx:e} sigy={sigy:e} sigt={sigt:e}")
+print(
+    f"  emittance_x={emittance_x:e} emittance_y={emittance_y:e} emittance_t={emittance_t:e}"
+)
+
+atol = 0.0  # ignored
+rtol = 1.8 * num_particles**-0.5  # from random sampling of a smooth distribution
+print(f"  rtol={rtol} (ignored: atol~={atol})")
+
+assert np.allclose(
+    [sigx, sigy, sigt, emittance_x, emittance_y, emittance_t],
+    [
+        1.559531175539e-3,
+        2.205510139392e-3,
+        1.0e-3,
+        1.0e-6,
+        2.0e-6,
+        1.0e-6,
+    ],
+    rtol=rtol,
+    atol=atol,
+)
+
+# particle-wise comparison against the periodic rectangular aperture boundary
+xmax = 1.5e-4
+ymax = 1.0e-4
+repeat_x = 1.0e-3
+repeat_y = 1.0e-3
+
+# kept particles, shifted to the primary
+dx = abs(final["position_x"]) - xmax
+dy = abs(final["position_y"]) - ymax
+
+print()
+print(f"  x_max={final['position_x'].max()}")
+print(f"  x_min={final['position_x'].min()}")
+assert np.less_equal(dx.max(), 0.0)
+
+print(f"  y_max={final['position_y'].max()}")
+print(f"  y_min={final['position_y'].min()}")
+assert np.less_equal(dy.max(), 0.0)
+
+# lost particles
+dx = abs(particles_lost["position_x"]) - xmax
+dy = abs(particles_lost["position_y"]) - ymax
+
+print()
+print(f"  x_max={particles_lost['position_x'].max()}")
+print(f"  x_min={particles_lost['position_x'].min()}")
+assert np.greater_equal(dx.max(), 0.0)
+
+print(f"  y_max={particles_lost['position_y'].max()}")
+print(f"  y_min={particles_lost['position_y'].min()}")
+assert np.greater_equal(dy.max(), 0.0)
+
diff --git a/examples/aperture/input_aperture_pepperpot.in b/examples/aperture/input_aperture_pepperpot.in
new file mode 100644
index 000000000..30ef15229
--- /dev/null
+++ b/examples/aperture/input_aperture_pepperpot.in
@@ -0,0 +1,52 @@
+###############################################################################
+# Particle Beam(s)
+###############################################################################
+beam.npart = 1000000
+beam.units = static
+beam.kin_energy = 250.0
+beam.charge = 1.0e-9
+beam.particle = proton
+beam.distribution = waterbag
+beam.lambdaX = 1.559531175539e-3
+beam.lambdaY = 2.205510139392e-3
+beam.lambdaT = 1.0e-3
+beam.lambdaPx = 6.41218345413e-4
+beam.lambdaPy = 9.06819680526e-4
+beam.lambdaPt = 1.0e-3
+beam.muxpx = 0.0
+beam.muypy = 0.0
+beam.mutpt = 0.0
+
+
+###############################################################################
+# Beamline: lattice elements and segments
+###############################################################################
+lattice.elements = monitor drift pepperpot monitor
+lattice.nslice = 1
+
+monitor.type = beam_monitor
+monitor.backend = h5
+
+drift.type = drift
+drift.ds = 0.123
+
+pepperpot.type = aperture
+pepperpot.shape = rectangular
+pepperpot.xmax = 1.5e-4
+pepperpot.ymax = 1.0e-4
+pepperpot.repeat_x = 1.0e-3
+pepperpot.repeat_y = 1.0e-3
+
+
+###############################################################################
+# Algorithms
+###############################################################################
+algo.particle_shape = 2
+algo.space_charge = false
+
+
+###############################################################################
+# Diagnostics
+###############################################################################
+diag.slice_step_diagnostics = true
+diag.backend = h5
diff --git a/examples/aperture/run_aperture_pepperpot.py b/examples/aperture/run_aperture_pepperpot.py
new file mode 100755
index 000000000..7e80c11d0
--- /dev/null
+++ b/examples/aperture/run_aperture_pepperpot.py
@@ -0,0 +1,69 @@
+#!/usr/bin/env python3
+#
+# Copyright 2022-2023 ImpactX contributors
+# Authors: Axel Huebl, Chad Mitchell
+# License: BSD-3-Clause-LBNL
+#
+# -*- coding: utf-8 -*-
+
+import amrex.space3d as amr
+from impactx import ImpactX, distribution, elements
+
+# work-around for https://github.com/ECP-WarpX/impactx/issues/499
+pp_amrex = amr.ParmParse("amrex")
+pp_amrex.add("the_arena_is_managed", 1)
+
+sim = ImpactX()
+
+# set numerical parameters and IO control
+sim.particle_shape = 2  # B-spline order
+sim.space_charge = False
+# sim.diagnostics = False  # benchmarking
+sim.slice_step_diagnostics = True
+sim.particle_lost_diagnostics_backend = "h5"
+
+# domain decomposition & space charge mesh
+sim.init_grids()
+
+# load a 250 MeV proton beam with an initial
+# horizontal rms emittance of 1 um and an
+# initial vertical rms emittance of 2 um
+kin_energy_MeV = 250.0  # reference energy
+bunch_charge_C = 1.0e-9  # used with space charge
+npart = 1000000  # number of macro particles
+
+#   reference particle
+ref = sim.particle_container().ref_particle()
+ref.set_charge_qe(1.0).set_mass_MeV(938.27208816).set_kin_energy_MeV(kin_energy_MeV)
+
+#   particle bunch
+distr = distribution.Waterbag(
+    lambdaX=1.559531175539e-3,
+    lambdaY=2.205510139392e-3,
+    lambdaT=1.0e-3,
+    lambdaPx=6.41218345413e-4,
+    lambdaPy=9.06819680526e-4,
+    lambdaPt=1.0e-3,
+)
+sim.add_particles(bunch_charge_C, distr, npart)
+
+# add beam diagnostics
+monitor = elements.BeamMonitor("monitor", backend="h5")
+
+# design the accelerator lattice
+sim.lattice.extend(
+    [
+        monitor,
+        elements.Drift(name="drift", ds=0.123),
+        elements.Aperture(
+            name="pepperpot", xmax=1.5e-4, ymax=1.0e-4, repeat_x=1.0e-3, repeat_y=1.0e-3, shape="rectangular"
+        ),
+        monitor,
+    ]
+)
+
+# run simulation
+sim.evolve()
+
+# clean shutdown
+sim.finalize()

From 8b1ed60c97fdee9b838d9a5ac0a315b45d527c24 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Fri, 11 Oct 2024 21:45:55 +0000
Subject: [PATCH 03/14] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 examples/aperture/README.rst                     | 3 +--
 examples/aperture/analysis_aperture_pepperpot.py | 1 -
 examples/aperture/run_aperture_pepperpot.py      | 7 ++++++-
 3 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/examples/aperture/README.rst b/examples/aperture/README.rst
index bb21b4ced..d8844744e 100644
--- a/examples/aperture/README.rst
+++ b/examples/aperture/README.rst
@@ -57,7 +57,7 @@ We run the following script to analyze correctness:
 Aperture Collimation with Periodic Masking
 ===========================================
 
-Proton beam undergoing collimation by a periodic array of rectangular apertures, such as those used in a pepperpot emittance measurement. 
+Proton beam undergoing collimation by a periodic array of rectangular apertures, such as those used in a pepperpot emittance measurement.
 
 We use a 250 MeV proton beam with a horizontal rms beam size of 1.56 mm and a vertical rms beam size of 2.21 mm.
 
@@ -101,4 +101,3 @@ We run the following script to analyze correctness:
    .. literalinclude:: analysis_aperture_pepperpot.py
       :language: python3
       :caption: You can copy this file from ``examples/aperture/analysis_aperture_pepperpot.py``.
-
diff --git a/examples/aperture/analysis_aperture_pepperpot.py b/examples/aperture/analysis_aperture_pepperpot.py
index 57f094fd4..292cca90a 100755
--- a/examples/aperture/analysis_aperture_pepperpot.py
+++ b/examples/aperture/analysis_aperture_pepperpot.py
@@ -105,4 +105,3 @@ def get_moments(beam):
 print(f"  y_max={particles_lost['position_y'].max()}")
 print(f"  y_min={particles_lost['position_y'].min()}")
 assert np.greater_equal(dy.max(), 0.0)
-
diff --git a/examples/aperture/run_aperture_pepperpot.py b/examples/aperture/run_aperture_pepperpot.py
index 7e80c11d0..cda6938a5 100755
--- a/examples/aperture/run_aperture_pepperpot.py
+++ b/examples/aperture/run_aperture_pepperpot.py
@@ -56,7 +56,12 @@
         monitor,
         elements.Drift(name="drift", ds=0.123),
         elements.Aperture(
-            name="pepperpot", xmax=1.5e-4, ymax=1.0e-4, repeat_x=1.0e-3, repeat_y=1.0e-3, shape="rectangular"
+            name="pepperpot",
+            xmax=1.5e-4,
+            ymax=1.0e-4,
+            repeat_x=1.0e-3,
+            repeat_y=1.0e-3,
+            shape="rectangular",
         ),
         monitor,
     ]

From 3323f51e38a4b1fe64540561ba95c14428c5979f Mon Sep 17 00:00:00 2001
From: Chad Mitchell <ChadMitchell@lbl.gov>
Date: Fri, 11 Oct 2024 17:12:53 -0700
Subject: [PATCH 04/14] Fix modular arithmetic.

---
 src/particles/elements/Aperture.H | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/particles/elements/Aperture.H b/src/particles/elements/Aperture.H
index 8d185ca03..701d44a20 100644
--- a/src/particles/elements/Aperture.H
+++ b/src/particles/elements/Aperture.H
@@ -110,14 +110,16 @@ namespace impactx
             // scale horizontal and vertical coordinates
             amrex::ParticleReal u;
             amrex::ParticleReal v;
+            amrex::ParticleReal dx = m_xmax;
+            amrex::ParticleReal dy = m_ymax;
             // first case: a single aperture
             if (m_repeat_x == 0.0 || m_repeat_y == 0.0) {
                u = x / m_xmax;
                v = y / m_ymax;
             // second case: a periodically repeated aperture
             } else {
-               u = std::fmod(x,m_repeat_x) / m_xmax;
-               v = std::fmod(y,m_repeat_y) / m_ymax;
+               u = (std::fmod(std::abs(x)+dx,m_repeat_x)-dx) / m_xmax;
+               v = (std::fmod(std::abs(y)+dy,m_repeat_y)-dy) / m_ymax;
             }
 
             // compare against the aperture boundary

From 110e07bf1cc3e07a13fd3ab1346d15d8741ceeae Mon Sep 17 00:00:00 2001
From: Chad Mitchell <ChadMitchell@lbl.gov>
Date: Fri, 11 Oct 2024 18:00:32 -0700
Subject: [PATCH 05/14] Add analysis script for testing.

---
 examples/CMakeLists.txt                       | 16 +++++++
 .../aperture/analysis_aperture_pepperpot.py   | 42 ++++++++++++-------
 2 files changed, 43 insertions(+), 15 deletions(-)

diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 19c4d33ca..f81cd3a3e 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -992,3 +992,19 @@ add_impactx_test(coupled-optics.py
     examples/coupled_optics/analysis_coupled_optics.py
     OFF  # no plot script yet
 )
+
+# Aperture with Periodic Masking #########################################################
+#
+# w/o space charge
+add_impactx_test(aperture-pepperpot
+    examples/aperture/input_aperture_pepperpot.in
+    ON   # ImpactX MPI-parallel
+    examples/aperture/analysis_aperture_pepperpot.py
+    OFF  # no plot script yet
+)
+add_impactx_test(aperture-pepperpot.py
+    examples/aperture/run_aperture_pepperpot.py
+    OFF  # ImpactX MPI-parallel
+    examples/aperture/analysis_aperture_pepperpot.py
+    OFF  # no plot script yet
+)
diff --git a/examples/aperture/analysis_aperture_pepperpot.py b/examples/aperture/analysis_aperture_pepperpot.py
index 292cca90a..c01409695 100755
--- a/examples/aperture/analysis_aperture_pepperpot.py
+++ b/examples/aperture/analysis_aperture_pepperpot.py
@@ -9,7 +9,6 @@
 import openpmd_api as io
 from scipy.stats import moment
 
-
 def get_moments(beam):
     """Calculate standard deviations of beam position & momenta
     and emittance values
@@ -80,28 +79,41 @@ def get_moments(beam):
 repeat_x = 1.0e-3
 repeat_y = 1.0e-3
 
-# kept particles, shifted to the primary
-dx = abs(final["position_x"]) - xmax
-dy = abs(final["position_y"]) - ymax
+
+# kept particles, shifted to the fundamental domain
+xshifted = abs(final["position_x"]) + xmax
+yshifted = abs(final["position_y"]) + ymax
+u = np.fmod(xshifted,repeat_x) - xmax
+v = np.fmod(yshifted,repeat_y) - ymax
+
+# difference from maximum aperture
+dx = abs(u) - xmax
+dy = abs(v) - ymax 
 
 print()
-print(f"  x_max={final['position_x'].max()}")
-print(f"  x_min={final['position_x'].min()}")
+print(f"  fundamental x_max={u.max()}")
+print(f"  fundamental x_min={u.min()}")
 assert np.less_equal(dx.max(), 0.0)
 
-print(f"  y_max={final['position_y'].max()}")
-print(f"  y_min={final['position_y'].min()}")
+print(f"  fundamental y_max={v.max()}")
+print(f"  fundamental y_min={v.min()}")
 assert np.less_equal(dy.max(), 0.0)
 
-# lost particles
-dx = abs(particles_lost["position_x"]) - xmax
-dy = abs(particles_lost["position_y"]) - ymax
+# lost particles, shifted to the fundamental domain
+xshifted = abs(particles_lost["position_x"]) - xmax
+yshifted = abs(particles_lost["position_y"]) - ymax
+u = np.fmod(xshifted,repeat_x) - xmax
+v = np.fmod(yshifted,repeat_y) - ymax
+
+# difference from maximum aperture
+dx = abs(u) - xmax
+dy = abs(v) - ymax 
 
 print()
-print(f"  x_max={particles_lost['position_x'].max()}")
-print(f"  x_min={particles_lost['position_x'].min()}")
+print(f"  fundamental x_max={u.max()}")
+print(f"  fundamental x_min={u.min()}")
 assert np.greater_equal(dx.max(), 0.0)
 
-print(f"  y_max={particles_lost['position_y'].max()}")
-print(f"  y_min={particles_lost['position_y'].min()}")
+print(f"  fundamental y_max={v.max()}")
+print(f"  fundamental y_min={v.min()}")
 assert np.greater_equal(dy.max(), 0.0)

From 1412a5a1261fb55614917e19cd36214c79e8ab78 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Sat, 12 Oct 2024 01:01:00 +0000
Subject: [PATCH 06/14] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 examples/aperture/analysis_aperture_pepperpot.py | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/examples/aperture/analysis_aperture_pepperpot.py b/examples/aperture/analysis_aperture_pepperpot.py
index c01409695..2e0a011db 100755
--- a/examples/aperture/analysis_aperture_pepperpot.py
+++ b/examples/aperture/analysis_aperture_pepperpot.py
@@ -9,6 +9,7 @@
 import openpmd_api as io
 from scipy.stats import moment
 
+
 def get_moments(beam):
     """Calculate standard deviations of beam position & momenta
     and emittance values
@@ -83,12 +84,12 @@ def get_moments(beam):
 # kept particles, shifted to the fundamental domain
 xshifted = abs(final["position_x"]) + xmax
 yshifted = abs(final["position_y"]) + ymax
-u = np.fmod(xshifted,repeat_x) - xmax
-v = np.fmod(yshifted,repeat_y) - ymax
+u = np.fmod(xshifted, repeat_x) - xmax
+v = np.fmod(yshifted, repeat_y) - ymax
 
 # difference from maximum aperture
 dx = abs(u) - xmax
-dy = abs(v) - ymax 
+dy = abs(v) - ymax
 
 print()
 print(f"  fundamental x_max={u.max()}")
@@ -102,12 +103,12 @@ def get_moments(beam):
 # lost particles, shifted to the fundamental domain
 xshifted = abs(particles_lost["position_x"]) - xmax
 yshifted = abs(particles_lost["position_y"]) - ymax
-u = np.fmod(xshifted,repeat_x) - xmax
-v = np.fmod(yshifted,repeat_y) - ymax
+u = np.fmod(xshifted, repeat_x) - xmax
+v = np.fmod(yshifted, repeat_y) - ymax
 
 # difference from maximum aperture
 dx = abs(u) - xmax
-dy = abs(v) - ymax 
+dy = abs(v) - ymax
 
 print()
 print(f"  fundamental x_max={u.max()}")

From 2b7dfbfd897c100daf15235fb8cdbd4a56d0b135 Mon Sep 17 00:00:00 2001
From: Chad Mitchell <46825199+cemitch99@users.noreply.github.com>
Date: Sat, 12 Oct 2024 11:27:31 -0700
Subject: [PATCH 07/14] Reduce particle population for tests.

---
 examples/aperture/analysis_aperture_pepperpot.py | 4 ++--
 examples/aperture/input_aperture_pepperpot.in    | 2 +-
 examples/aperture/run_aperture_pepperpot.py      | 2 +-
 3 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/examples/aperture/analysis_aperture_pepperpot.py b/examples/aperture/analysis_aperture_pepperpot.py
index 2e0a011db..4c9c59cae 100755
--- a/examples/aperture/analysis_aperture_pepperpot.py
+++ b/examples/aperture/analysis_aperture_pepperpot.py
@@ -43,7 +43,7 @@ def get_moments(beam):
 particles_lost = series_lost.iterations[0].particles["beam"].to_df()
 
 # compare number of particles
-num_particles = 1000000
+num_particles = 100000
 assert num_particles == len(initial)
 # we lost particles in apertures
 assert num_particles > len(final)
@@ -57,7 +57,7 @@ def get_moments(beam):
 )
 
 atol = 0.0  # ignored
-rtol = 1.8 * num_particles**-0.5  # from random sampling of a smooth distribution
+rtol = 2.0 * num_particles**-0.5  # from random sampling of a smooth distribution
 print(f"  rtol={rtol} (ignored: atol~={atol})")
 
 assert np.allclose(
diff --git a/examples/aperture/input_aperture_pepperpot.in b/examples/aperture/input_aperture_pepperpot.in
index 30ef15229..eaec8f1ed 100644
--- a/examples/aperture/input_aperture_pepperpot.in
+++ b/examples/aperture/input_aperture_pepperpot.in
@@ -1,7 +1,7 @@
 ###############################################################################
 # Particle Beam(s)
 ###############################################################################
-beam.npart = 1000000
+beam.npart = 100000
 beam.units = static
 beam.kin_energy = 250.0
 beam.charge = 1.0e-9
diff --git a/examples/aperture/run_aperture_pepperpot.py b/examples/aperture/run_aperture_pepperpot.py
index cda6938a5..9523abcb5 100755
--- a/examples/aperture/run_aperture_pepperpot.py
+++ b/examples/aperture/run_aperture_pepperpot.py
@@ -30,7 +30,7 @@
 # initial vertical rms emittance of 2 um
 kin_energy_MeV = 250.0  # reference energy
 bunch_charge_C = 1.0e-9  # used with space charge
-npart = 1000000  # number of macro particles
+npart = 100000  # number of macro particles
 
 #   reference particle
 ref = sim.particle_container().ref_particle()

From a5714ec6882909f1c5676feb373ca418d6016d38 Mon Sep 17 00:00:00 2001
From: Chad Mitchell <46825199+cemitch99@users.noreply.github.com>
Date: Tue, 15 Oct 2024 15:58:06 -0700
Subject: [PATCH 08/14] Update examples/aperture/README.rst

Correct type in example README.
---
 examples/aperture/README.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/aperture/README.rst b/examples/aperture/README.rst
index d8844744e..daa65ff11 100644
--- a/examples/aperture/README.rst
+++ b/examples/aperture/README.rst
@@ -63,7 +63,7 @@ We use a 250 MeV proton beam with a horizontal rms beam size of 1.56 mm and a ve
 
 After a short drift of 0.123 m, the beam intercepts a periodic array of apertures of period 1 mm (in the horizontal and vertical), each of which is 0.15 mm x 0.1 mm in size.
 
-In this test, the initial values of :math:`\sigma_x`, :math:`\sigma_y`, :math:`\sigma_t`, :math:`\epsilon_x`, :math:`\epsilon_y`, and :math:`\epsilon_t` must agree with nomin>
+In this test, the initial values of :math:`\sigma_x`, :math:`\sigma_y`, :math:`\sigma_t`, :math:`\epsilon_x`, :math:`\epsilon_y`, and :math:`\epsilon_t` must agree with nominal values.
 The test fails if:
 
 * any of the final coordinates for the valid (not lost) particles lie outside the aperture boundary or

From 099e19ef4c8ac8984417b3a9ff3f06cb4f10f72b Mon Sep 17 00:00:00 2001
From: Chad Mitchell <46825199+cemitch99@users.noreply.github.com>
Date: Thu, 17 Oct 2024 09:47:19 -0700
Subject: [PATCH 09/14] Update Aperture.H

Generalize periodic aperture to allow periodicity in one direction only.
---
 src/particles/elements/Aperture.H | 19 +++++++++----------
 1 file changed, 9 insertions(+), 10 deletions(-)

diff --git a/src/particles/elements/Aperture.H b/src/particles/elements/Aperture.H
index 701d44a20..d60809aa2 100644
--- a/src/particles/elements/Aperture.H
+++ b/src/particles/elements/Aperture.H
@@ -107,20 +107,19 @@ namespace impactx
             // shift due to alignment errors of the element
             shift_in(x, y, px, py);
 
-            // scale horizontal and vertical coordinates
+            // define intermediate variables
             amrex::ParticleReal u;
             amrex::ParticleReal v;
             amrex::ParticleReal dx = m_xmax;
             amrex::ParticleReal dy = m_ymax;
-            // first case: a single aperture
-            if (m_repeat_x == 0.0 || m_repeat_y == 0.0) {
-               u = x / m_xmax;
-               v = y / m_ymax;
-            // second case: a periodically repeated aperture
-            } else {
-               u = (std::fmod(std::abs(x)+dx,m_repeat_x)-dx) / m_xmax;
-               v = (std::fmod(std::abs(y)+dy,m_repeat_y)-dy) / m_ymax;
-            }
+
+            // if the aperture is periodic, shift x,y coordinates to the fundamental domain
+            u = (m_repeat_x==0.0) ? x : (std::fmod(std::abs(x)+dx,m_repeat_x)-dx);
+            v = (m_repeat_y==0.0) ? y : (std::fmod(std::abs(y)+dy,m_repeat_y)-dy);
+
+            // scale horizontal and vertical coordinates
+            u = u / m_xmax;
+            v = v / m_ymax;
 
             // compare against the aperture boundary
             switch (m_shape)

From 5c67be09a19f7cd4b5c330b1dadce7b579d637f6 Mon Sep 17 00:00:00 2001
From: Chad Mitchell <46825199+cemitch99@users.noreply.github.com>
Date: Thu, 17 Oct 2024 10:04:30 -0700
Subject: [PATCH 10/14] Rename analysis_aperture_pepperpot.py to
 analysis_aperture_periodic.py

Rename example:  analysis script.
---
 ...alysis_aperture_pepperpot.py => analysis_aperture_periodic.py} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename examples/aperture/{analysis_aperture_pepperpot.py => analysis_aperture_periodic.py} (100%)

diff --git a/examples/aperture/analysis_aperture_pepperpot.py b/examples/aperture/analysis_aperture_periodic.py
similarity index 100%
rename from examples/aperture/analysis_aperture_pepperpot.py
rename to examples/aperture/analysis_aperture_periodic.py

From 520b5c909db87e054c907ff07f3f0d4a78007169 Mon Sep 17 00:00:00 2001
From: Chad Mitchell <46825199+cemitch99@users.noreply.github.com>
Date: Thu, 17 Oct 2024 10:05:16 -0700
Subject: [PATCH 11/14] Rename input_aperture_pepperpot.in to
 input_aperture_periodic.in

Rename example:  C++ input
---
 .../{input_aperture_pepperpot.in => input_aperture_periodic.in}   | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename examples/aperture/{input_aperture_pepperpot.in => input_aperture_periodic.in} (100%)

diff --git a/examples/aperture/input_aperture_pepperpot.in b/examples/aperture/input_aperture_periodic.in
similarity index 100%
rename from examples/aperture/input_aperture_pepperpot.in
rename to examples/aperture/input_aperture_periodic.in

From 90497edcffd3e499ba5c2329914e1743aa606338 Mon Sep 17 00:00:00 2001
From: Chad Mitchell <46825199+cemitch99@users.noreply.github.com>
Date: Thu, 17 Oct 2024 10:05:43 -0700
Subject: [PATCH 12/14] Rename run_aperture_pepperpot.py to
 run_aperture_periodic.py

Rename example:  Python input
---
 .../{run_aperture_pepperpot.py => run_aperture_periodic.py}       | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename examples/aperture/{run_aperture_pepperpot.py => run_aperture_periodic.py} (100%)

diff --git a/examples/aperture/run_aperture_pepperpot.py b/examples/aperture/run_aperture_periodic.py
similarity index 100%
rename from examples/aperture/run_aperture_pepperpot.py
rename to examples/aperture/run_aperture_periodic.py

From c05e2c439087677794b2c5435b8738d3ec03aca5 Mon Sep 17 00:00:00 2001
From: Chad Mitchell <46825199+cemitch99@users.noreply.github.com>
Date: Thu, 17 Oct 2024 10:07:25 -0700
Subject: [PATCH 13/14] Update CMakeLists.txt

Rename example:  CMakeLists
---
 examples/CMakeLists.txt | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index f81cd3a3e..8da70838a 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -996,15 +996,15 @@ add_impactx_test(coupled-optics.py
 # Aperture with Periodic Masking #########################################################
 #
 # w/o space charge
-add_impactx_test(aperture-pepperpot
-    examples/aperture/input_aperture_pepperpot.in
+add_impactx_test(aperture-periodic
+    examples/aperture/input_aperture_periodic.in
     ON   # ImpactX MPI-parallel
-    examples/aperture/analysis_aperture_pepperpot.py
+    examples/aperture/analysis_aperture_periodic.py
     OFF  # no plot script yet
 )
-add_impactx_test(aperture-pepperpot.py
-    examples/aperture/run_aperture_pepperpot.py
+add_impactx_test(aperture-periodic.py
+    examples/aperture/run_aperture_periodic.py
     OFF  # ImpactX MPI-parallel
-    examples/aperture/analysis_aperture_pepperpot.py
+    examples/aperture/analysis_aperture_periodic.py
     OFF  # no plot script yet
 )

From 914b0f7d357e1af3933f47a85c397bc0e1ee1a52 Mon Sep 17 00:00:00 2001
From: Chad Mitchell <46825199+cemitch99@users.noreply.github.com>
Date: Thu, 17 Oct 2024 10:12:03 -0700
Subject: [PATCH 14/14] Update README.rst

Rename example:  README documentation
---
 examples/aperture/README.rst | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/examples/aperture/README.rst b/examples/aperture/README.rst
index daa65ff11..4fac877e7 100644
--- a/examples/aperture/README.rst
+++ b/examples/aperture/README.rst
@@ -52,7 +52,7 @@ We run the following script to analyze correctness:
       :caption: You can copy this file from ``examples/aperture/analysis_aperture.py``.
 
 
-.. _examples-aperture-pepperpot:
+.. _examples-aperture-periodic:
 
 Aperture Collimation with Periodic Masking
 ===========================================
@@ -67,37 +67,37 @@ In this test, the initial values of :math:`\sigma_x`, :math:`\sigma_y`, :math:`\
 The test fails if:
 
 * any of the final coordinates for the valid (not lost) particles lie outside the aperture boundary or
-* any of the lost particles are inside the aperture boundary or
-* if the sum of lost and kept particles is not equal to the initial particles
+* any of the lost particles lie inside the aperture boundary or
+* if the sum of the numbers of lost and retained particles is not equal to the number of initial particles
 
 
 Run
 ---
 
-This example can be run as a Python script (``python3 run_aperture_pepperpot.py``) or with an app with an input file (``impactx input_aperture_pepperpot.in``).
+This example can be run as a Python script (``python3 run_aperture_periodic.py``) or with an app with an input file (``impactx input_aperture_periodic.in``).
 Each can also be prefixed with an `MPI executor <https://www.mpi-forum.org>`__, such as ``mpiexec -n 4 ...`` or ``srun -n 4 ...``, depending on the system.
 
 .. tab-set::
 
    .. tab-item:: Python Script
 
-       .. literalinclude:: run_aperture_pepperpot.py
+       .. literalinclude:: run_aperture_periodic.py
           :language: python3
-          :caption: You can copy this file from ``examples/aperture/run_aperture_pepperpot.py``.
+          :caption: You can copy this file from ``examples/aperture/run_aperture_periodic.py``.
 
    .. tab-item:: App Input File
 
-       .. literalinclude:: input_aperture_pepperpot.in
+       .. literalinclude:: input_aperture_periodic.in
           :language: ini
-          :caption: You can copy this file from ``examples/aperture/input_aperture_pepperpot.in``.
+          :caption: You can copy this file from ``examples/aperture/input_aperture_periodic.in``.
 
 Analyze
 -------
 
 We run the following script to analyze correctness:
 
-.. dropdown:: Script ``analysis_aperture_pepperpot.py``
+.. dropdown:: Script ``analysis_aperture_periodic.py``
 
-   .. literalinclude:: analysis_aperture_pepperpot.py
+   .. literalinclude:: analysis_aperture_periodic.py
       :language: python3
-      :caption: You can copy this file from ``examples/aperture/analysis_aperture_pepperpot.py``.
+      :caption: You can copy this file from ``examples/aperture/analysis_aperture_periodic.py``.