Fix Repeated Plasma Lens: Start (ECP-WarpX#4220)

* Fix Repeated Plasma Lens: Start

The current implementation has problems finding the right
index of the current plasma lens if:
- the lens started before z=0
- the lens did not start between z=0...first-period

This should fix it.

* Document Start, Define an End

Docs/source/usage/parameters.rst

@@ -1474,8 +1474,8 @@ Applied to Particles
 
       Note that the position is defined in Cartesian coordinates, as a function of (x,y,z), even for RZ.
 
-    * ``repeated_plasma_lens``: apply a series of plasma lenses. The properties of the lenses are defined in the
-      lab frame by the input parameters:
+    * ``repeated_plasma_lens``: apply a series of plasma lenses.
+      The properties of the lenses are defined in the lab frame by the input parameters:
 
         * ``repeated_plasma_lens_period``, the period length of the repeat, a single float number,
 
@@ -1489,6 +1489,9 @@ Applied to Particles
         * ``repeated_plasma_lens_strengths_B``, the magnetic focusing strength of each lens, an array of floats, when
           ``particles.B_ext_particle_init_style`` is set to ``repeated_plasma_lens``.
 
+      The repeated lenses are only defined for :math:`z > 0`.
+      Once the number of lenses specified in the input are exceeded, the repeated lens stops.
+
       The applied field is uniform longitudinally (along z) with a hard edge,
       where residence corrections are used for more accurate field calculation. On the time step when a particle enters
       or leaves each lens, the field applied is scaled by the fraction of the time step spent within the lens.

Source/Particles/Gather/GetExternalFields.H

@@ -145,6 +145,7 @@ struct GetExternalEBField
             const amrex::ParticleReal gamma = std::sqrt(1._prt + (uxp*uxp + uyp*uyp + uzp*uzp)*inv_c2);
             const amrex::ParticleReal vzp = uzp/gamma;
 
+            // the current slice in z between now and the next time step
             amrex::ParticleReal zl = z;
             amrex::ParticleReal zr = z + vzp*m_dt;
 
@@ -153,27 +154,31 @@ struct GetExternalEBField
                 zr = m_gamma_boost*zr + m_uz_boost*(m_time + m_dt);
             }
 
-            // This assumes that zl > 0.
-            int i_lens = static_cast<int>(std::floor(zl/m_repeated_plasma_lens_period));
-            i_lens = i_lens % m_n_lenses;
-            amrex::ParticleReal const lens_start = m_repeated_plasma_lens_starts[i_lens] + i_lens*m_repeated_plasma_lens_period;
-            amrex::ParticleReal const lens_end = lens_start + m_repeated_plasma_lens_lengths[i_lens];
-
-            // Calculate the residence correction
-            // frac will be 1 if the step is completely inside the lens, between 0 and 1
-            // when entering or leaving the lens, and otherwise 0.
-            // This accounts for the case when particles step over the element without landing in it.
-            // This assumes that vzp > 0.
-            amrex::ParticleReal const zl_bounded = std::min(std::max(zl, lens_start), lens_end);
-            amrex::ParticleReal const zr_bounded = std::min(std::max(zr, lens_start), lens_end);
-            amrex::ParticleReal const frac = ((zr - zl) == 0._rt ? 1._rt : (zr_bounded - zl_bounded)/(zr - zl));
-
-            // Note that "+=" is used since the fields may have been set above
-            // if a different E or Btype was specified.
-            Ex += x*frac*m_repeated_plasma_lens_strengths_E[i_lens];
-            Ey += y*frac*m_repeated_plasma_lens_strengths_E[i_lens];
-            Bx += +y*frac*m_repeated_plasma_lens_strengths_B[i_lens];
-            By += -x*frac*m_repeated_plasma_lens_strengths_B[i_lens];
+            // the plasma lens periods do not start before z=0
+            if (zl > 0) {
+                // find which is the next lens
+                int i_lens = static_cast<int>(std::floor(zl/m_repeated_plasma_lens_period));
+                if (i_lens < m_n_lenses) {
+                    amrex::ParticleReal const lens_start = m_repeated_plasma_lens_starts[i_lens] + i_lens*m_repeated_plasma_lens_period;
+                    amrex::ParticleReal const lens_end = lens_start + m_repeated_plasma_lens_lengths[i_lens];
+
+                    // Calculate the residence correction
+                    // frac will be 1 if the step is completely inside the lens, between 0 and 1
+                    // when entering or leaving the lens, and otherwise 0.
+                    // This accounts for the case when particles step over the element without landing in it.
+                    // This assumes that vzp > 0.
+                    amrex::ParticleReal const zl_bounded = std::min(std::max(zl, lens_start), lens_end);
+                    amrex::ParticleReal const zr_bounded = std::min(std::max(zr, lens_start), lens_end);
+                    amrex::ParticleReal const frac = ((zr - zl) == 0._rt ? 1._rt : (zr_bounded - zl_bounded)/(zr - zl));
+
+                    // Note that "+=" is used since the fields may have been set above
+                    // if a different E or Btype was specified.
+                    Ex += x*frac*m_repeated_plasma_lens_strengths_E[i_lens];
+                    Ey += y*frac*m_repeated_plasma_lens_strengths_E[i_lens];
+                    Bx += +y*frac*m_repeated_plasma_lens_strengths_B[i_lens];
+                    By += -x*frac*m_repeated_plasma_lens_strengths_B[i_lens];
+                }
+            }
 
         }