fix(phys): accumulate epitaxy strain per layer

CHANGELOG.md

@@ -37,6 +37,11 @@ Scientific and analysis additions:
 
 Scientific correctness and performance fixes:
 
+- Corrected ``EpitaxyInterface`` to accumulate the strain-coupled height of
+  every structural layer. The previous calculation advanced a complete unit
+  cell using only its final layer's strain, creating a layer-cycle-periodic
+  displacement and artificial quasi-Bragg intensity at systematic extinctions
+  for broad interfaces.
 - Optical profiles now preserve areal optical content under surface-normal
   strain, keep ionic forward scattering factors for charged species, avoid
   merging layers beyond the requested z tolerance, and remain finite at the

doc/source/ctr_structure_factors.rst

@@ -240,6 +240,21 @@ the linear interpolation
    z_m(\kappa,0)
    = z_{0,m} + \kappa\left(z_1-z_{0,m}\right).
 
+For the fully strain-coupled discrete profile, every structural layer supplies
+its own physical interval.  If :math:`d_{\mathrm{top},i}` and
+:math:`d_{\mathrm{bottom},i}` are the upper- and lower-lattice intervals for
+layer :math:`i`, the coupled interval is
+
+.. math::
+
+   d_{1,i}
+   = P_i d_{\mathrm{top},i}
+     + (1-P_i)d_{\mathrm{bottom},i},
+
+and the layer origins are the cumulative sum of these intervals.  The
+calculation must therefore accumulate all structural layers independently;
+one layer's strain cannot represent a complete multi-layer unit-cell cycle.
+
 The ``offset`` parameter :math:`o` is expressed in fractional lower-bulk
 :math:`c` coordinates.  Its physical displacement is
 

doc/source/release_notes.rst

@@ -9,11 +9,22 @@ Unreleased (2026-07-19)
 
 Scientific and analysis additions:
 
+- Added interface-based optical ``delta``/``beta`` profiles and layered wavefield calculations for s- and p-polarized X-rays. Scalar and multidimensional angle inputs are supported, and specular reflectivity can be calculated for s, p, or unpolarized incidence. This provides the optical wavefield foundation for future distorted-wave Born approximation (DWBA) calculations; the full distorted-wave scattering amplitude is not included yet.
 - Added optional CTR intensity-resolution modeling with constant or gamma-dependent box and Gaussian functions. Calculations can convolve irregular existing L points or sample the crystal structure factor with deterministic quadrature.
+- Added Poisson-distributed surface occupancies and coherent out-of-plane epitaxy/strain coupling for film-interface models.
 - Added py3Dmol atom-sphere rendering for Jupyter notebooks. ``plot3d`` now selects py3Dmol automatically in a notebook, can be directed to either py3Dmol or Mayavi explicitly, and can incrementally add unit cells to a shared viewer. Covalent radii are interpreted consistently by both backends and can be adjusted with the dimensionless ``radius_scale`` parameter.
 - ``PoissonSurface``'s basis now requires three parameters (``W``, ``alpha``, ``offset``) instead of two. The previous two-parameter form was only ever used in test code, so no migration path is provided; a saved ``.xtal``/ ``.xpr`` file with a two-parameter ``PoissonSurface`` basis will fail to load.
 - Wyckoff-parameterized fit values are now interpreted as absolute atomic positions instead of deltas relative to the Wyckoff position. Wyckoff parameter fitting was development-only and never part of a release, so no migration is provided.
 
+Scientific correctness and performance fixes:
+
+- Corrected ``EpitaxyInterface`` to accumulate the strain-coupled height of every structural layer. The previous calculation advanced a complete unit cell using only its final layer's strain, creating a layer-cycle-periodic displacement and artificial quasi-Bragg intensity at systematic extinctions for broad interfaces.
+- Optical profiles now preserve areal optical content under surface-normal strain, keep ionic forward scattering factors for charged species, avoid merging layers beyond the requested z tolerance, and remain finite at the exact p-polarized critical-angle limit.
+- Optical reflectivity and full wavefield calculations now evaluate angle arrays with NumPy-vectorized kernels while preserving the caller's input shape.
+- Corrected film and epitaxy-interface anchoring, stacking, and support ownership across strained interfaces.
+- Fixed missing support in one-dimensional electron-density calculations and corrected atomic-coordinate stacking when splitting unit cells into layers.
+- Added a C++ electron-density backend and bounded caches for atomic form factors, anomalous scattering factors, and accelerated form-factor lookup.
+
 A ***critical bug*** was fixed that affects bulk CTR calculations:
 
 - ``UnitCell.F_bulk``'s semi-infinite geometric lattice sum used the raw, untransformed ``l`` index instead of the index converted by ``refHKLTransform`` when computing the out-of-plane attenuation phase. This was only correct for the default case where a component uses its own bulk cell as the reference (no ``reference_uc`` set); any explicit ``reference_uc`` whose out-of-plane reciprocal axis differs from the bulk's — including a plain scale difference between the reference and bulk out-of-plane axis length, not only a rotated or reindexed reference — gave incorrect bulk structure-factor amplitudes. This bug was present in both the accelerated (numba/C++) and plain-Python code paths in all previous released versions, up to and including v1.5.0. See the CTR structure-factor documentation for details.

orgui/datautils/xrayutils/CTRfilm.py

@@ -763,51 +763,81 @@ def createInterfaceCells(self):
 
             mat_0 = np.vstack((np.identity(3).T, np.array([0, 0, 0]))).T
 
-            ratio_top = a3_bottom / a3_top
-            ratio_bottom = 1 / ratio_top
             coupled_top = [[] for _ in self.top_layers]
             coupled_bottom = [[] for _ in self.bottom_layers]
             ideal_top = [[] for _ in self.top_layers]
             ideal_bottom = [[] for _ in self.bottom_layers]
             coupled_height = 0.0
+            top_layer_space = np.diff(
+                self.layerpos,
+                append=self.layerpos[0] + 1.0,
+            )
+            bottom_layerpos = _unwrapped_layer_positions(
+                np.array(
+                    [
+                        self.uc_bottom.layerpos[layer_id]
+                        for layer_id in self._layer_ids
+                    ]
+                ),
+                self._layer_order_indices,
+            )
+            bottom_layer_space = np.diff(
+                bottom_layerpos,
+                append=bottom_layerpos[0] + 1.0,
+            )
 
             for cycle_index, (p_t, p_b) in enumerate(
                 zip(probability_top, probability_bottom)
             ):
-                top_strains = p_t + ratio_top * p_b
-                bottom_strains = ratio_bottom * p_t + p_b
-                physical_top_index = np.flatnonzero(
-                    self._layer_order_indices == n_layers - 1
-                )[0]
-                cell_height = a3_top * top_strains[physical_top_index]
+                # Accumulate the coherent height one structural layer at a
+                # time. Using a single layer's strain for the whole unit-cell
+                # cycle creates a displacement modulation with the layer-cycle
+                # period and leaks intensity into systematic extinctions.
+                physical_layer_heights = (
+                    p_t * a3_top * top_layer_space
+                    + p_b * a3_bottom * bottom_layer_space
+                )
+                top_strains = physical_layer_heights / (
+                    a3_top * top_layer_space
+                )
+                bottom_strains = physical_layer_heights / (
+                    a3_bottom * bottom_layer_space
+                )
 
                 for i, (uc_t, uc_b) in enumerate(
                     zip(self.top_layers, self.bottom_layers)
                 ):
                     mat_top_i = np.copy(mat_0)
                     top_strain_and_h = top_strains[i]
                     mat_top_i[2, 2] = top_strain_and_h
-                    relative_layer_position = self.layerpos[i] - self.layerpos[0]
+                    relative_top_position = (
+                        self.layerpos[i] - self.layerpos[0]
+                    )
                     top_layer_offset = (
-                        relative_layer_position
+                        relative_top_position
                         - (self.uc_top.layerpos[self.layer_order[i]])
                     )
                     mat_top_i[2, 3] = (
                         coupled_height / a3_top
-                        + top_strain_and_h * top_layer_offset
+                        - top_strain_and_h
+                        * self.uc_top.layerpos[self.layer_order[i]]
                     )
                     coupled_top[i].append(mat_top_i)
 
                     mat_bottom_i = np.copy(mat_0)
                     bottom_strain_and_h = bottom_strains[i]
                     mat_bottom_i[2, 2] = bottom_strain_and_h
+                    relative_bottom_position = (
+                        bottom_layerpos[i] - bottom_layerpos[0]
+                    )
                     bottom_layer_offset = (
-                        relative_layer_position
+                        relative_bottom_position
                         - (self.uc_bottom.layerpos[self.layer_order[i]])
                     )
                     mat_bottom_i[2, 3] = (
                         coupled_height / a3_bottom
-                        + bottom_strain_and_h * bottom_layer_offset
+                        - bottom_strain_and_h
+                        * self.uc_bottom.layerpos[self.layer_order[i]]
                     )
                     coupled_bottom[i].append(mat_bottom_i)
 
@@ -818,7 +848,7 @@ def createInterfaceCells(self):
                     ideal_bottom_i = np.copy(mat_0)
                     ideal_bottom_i[2, 3] = cycle_index + bottom_layer_offset
                     ideal_bottom[i].append(ideal_bottom_i)
-                coupled_height += cell_height
+                    coupled_height += physical_layer_heights[i]
 
             loc_rescaled = loc - unitcells[0]
             uc_no_loc = int(np.floor(loc_rescaled)) // n_layers

orgui/datautils/xrayutils/test/test_CTRcalc.py

@@ -2028,6 +2028,42 @@ def test_epitaxy_strain_coupled_field_is_anchored_to_deep_bulk(self):
             interface.height_absolute - independent.height_absolute,
         )
 
+    def test_epitaxy_accumulates_each_structural_layer_height(self):
+        interface = self.make_epitaxy_interface(strain_coupling=1.0)
+        n_layers = len(interface.layers)
+        unitcells = np.arange(-interface.fixed_ucs, interface.fixed_ucs)
+        probability_top = interface.profile.occupancy(
+            unitcells,
+            n_layers,
+        ).reshape((-1, n_layers))
+        probability_bottom = 1.0 - probability_top
+        layer_space = np.diff(
+            interface.layerpos,
+            append=interface.layerpos[0] + 1.0,
+        )
+        expected_layer_heights = (
+            probability_top * interface.uc_top.a[2] * layer_space
+            + probability_bottom * interface.uc_bottom.a[2] * layer_space
+        )
+
+        generated_origins = []
+        for cycle_index in range(probability_top.shape[0]):
+            for layer_index, layer in enumerate(interface.top_layers):
+                matrix = layer.coherentDomainMatrix[cycle_index]
+                layer_id = interface.layer_order[layer_index]
+                generated_origins.append(
+                    (
+                        matrix[2, 2] * interface.uc_top.layerpos[layer_id]
+                        + matrix[2, 3]
+                    )
+                    * interface.uc_top.a[2]
+                )
+
+        np.testing.assert_allclose(
+            np.diff(generated_origins),
+            expected_layer_heights.ravel()[:-1],
+        )
+
     def test_epitaxy_strain_coupling_interpolates_generated_positions(self):
         independent = self.make_epitaxy_interface(strain_coupling=0.0)
         halfway = self.make_epitaxy_interface(strain_coupling=0.5)
@@ -2373,40 +2409,70 @@ def test_epitaxy_split_bulk_domains_feed_all_profile_paths(self):
         np.testing.assert_allclose(
             structure_factor,
             [
-                12.290590935844 + 23.878872280467j,
-                -1.767926251260 - 0.197899416448j,
+                9.783621364160 + 25.055258758458j,
+                -0.401679480851 - 1.415519651389j,
             ],
             rtol=1e-8,
         )
         np.testing.assert_allclose(
             density[[0, 100, 200, 300, 400]],
             [
                 0.0,
-                -3.006212046217e-03 - 1.665286072977e-05j,
-                9.802746914877e-01 + 3.098533417278e-03j,
-                1.082893246501e-05,
+                -3.093100733396e-03 - 1.727080407634e-05j,
+                7.183801385126e-01 + 2.427528867636e-03j,
+                8.622439679978e-05,
                 0.0,
             ],
             atol=1e-9,
         )
         np.testing.assert_allclose(
             np.sum(density),
-            40.033115045982 + 0.034943462837j,
+            40.028295180092 + 0.034223926465j,
             rtol=1e-8,
         )
         mid_row = optical_profile.shape[0] // 2
         np.testing.assert_allclose(
             optical_profile[[0, mid_row, -1]],
             [
                 [-12.0, -7.72508074e-10, -7.18561204e-13],
-                [-0.209418853, 2.94939896e-06, 2.74343239e-09],
+                [0.0, 2.94939896e-06, 2.74343239e-09],
                 [10.0, 0.0, 0.0],
             ],
             atol=1e-9,
         )
 
 
 class TestLegacyLayeredCTR(unittest.TestCase):
+    def test_broad_interface_does_not_create_rutile_quasi_bragg_peak(self):
+        fixture_root = os.path.join(os.path.dirname(__file__), "testdata")
+        xtal = CTRcalc.SXRDCrystal.fromFile(
+            os.path.join(fixture_root, "RuO2_TiO2_Poisson_etching.xtal")
+        )
+        surface = xtal["RuO2surface"]
+        film = xtal["RuO2"]
+        interface = xtal["TiO2toRuO2"]
+        surface_width = 10.0 / (
+            film.unitcell.a[2] / len(film.layers)
+        )
+        interface_width = 10.0 / interface.uc_bottom.a[2]
+        surface.basis[:] = [surface_width, 1.0, 0.0]
+        film.basis[0] = 105.0
+        interface.basis[:] = [interface_width, 0.0, 1.0, 0.0]
+
+        l_values = np.array([3.9, 4.0, 4.1])
+        amplitudes = np.abs(
+            xtal.F(
+                np.ones(3, dtype=np.float64),
+                np.zeros(3, dtype=np.float64),
+                l_values,
+            )
+        )
+
+        self.assertLess(
+            amplitudes[1],
+            np.mean(amplitudes[[0, 2]]),
+        )
+
     def test_legacy_xtal_uses_corrected_interface_support(self):
         fixture_root = os.path.join(os.path.dirname(__file__), "testdata")
         xtal_path = os.path.join(
@@ -2440,8 +2506,9 @@ def test_legacy_xtal_uses_corrected_interface_support(self):
         # UNTESTED reference values: computed with this module's own
         # CTRcalc.SXRDCrystal.F() after the interfacial-strain fix in
         # EpitaxyInterface (commit beee738), the F_bulk lattice-sum
-        # denominator fix, and the deep-bulk strain-field anchoring fix, but
-        # not independently verified against real data or an analytic result
+        # denominator fix, the deep-bulk strain-field anchoring fix, and the
+        # cumulative per-structural-layer interface-height fix, but not
+        # independently verified against real data or an analytic result
         # -- unlike test_simple_unitcell_bulk_structure_factor_matches_
         # reference below, which does compare against real fit-derived
         # intensities for the un-layered unit-cell + bulk case.
@@ -2464,13 +2531,13 @@ def test_legacy_xtal_uses_corrected_interface_support(self):
         expected = np.array(
             [
                 17238.5440045948,
-                7001.3340402313,
-                3852.2553000918,
-                2726.9468772925,
-                2175.6251991885,
-                1855.7582644482,
-                1648.1329478054,
-                1501.1768061163,
+                7001.3031826720,
+                3852.1926461220,
+                2726.8587419552,
+                2175.5183472486,
+                1855.6392611730,
+                1648.0076716502,
+                1501.0502790582,
             ]
         )
         np.testing.assert_allclose(calculated, expected, rtol=2e-5, atol=2e-5)