diff --git a/webbpsf/constants.py b/webbpsf/constants.py
index 4222b398..6351b894 100644
--- a/webbpsf/constants.py
+++ b/webbpsf/constants.py
@@ -401,3 +401,11 @@
 # This is a rough approximation of a detector-position-dependent phenomenon
 MIRI_CRUCIFORM_INNER_RADIUS_PIX = 12
 MIRI_CRUCIFORM_RADIAL_SCALEFACTOR = 0.005   # Brightness factor for the diffuse circular halo
+
+# Various parameters to model the more complex parts of the cruciform
+# There's a fair amount of ad-hoc tuning still needed here to match real data
+MIRI_CRUCIFORM_PEAK_REFWAVE = 5.5  #  # See Gaspar et al. sections 4.2 and 4.3
+MIRI_CRUCIFORM_PEAKS_LOC = [12/1.1, 16.5,   19.5, 23.5, 27.9, 33, 39.6]  # pixels, at 5.5 microns, for the centered PSF; slightly tuned to in-flight ePSFs
+MIRI_CRUCIFORM_PEAKS_AMP = [1.8e-4, 1.0e-4, 1.0e-4, 9e-5, 4e-5, 2e-5, 1e-5]  # read off from fig 12 in Gaspar et al.  Slightly tuned to in-flight ePSFs.
+MIRI_CRUCIFORM_PEAKS_AMP_ADJUST = 20        # ad hoc scale factor, for units conversion from fig 12 to how it's normalized here 
+MIRI_CRUCIFORM_PEAKS_SIGMA = 1              # how wide are the peaks/bumps along the crucifor
diff --git a/webbpsf/detectors.py b/webbpsf/detectors.py
index fe4aa86e..631757a0 100644
--- a/webbpsf/detectors.py
+++ b/webbpsf/detectors.py
@@ -328,11 +328,33 @@ def _make_miri_scattering_kernel_2d(in_psf, kernel_amp, oversample=1, wavelength
     kernel_x[np.abs(x) < constants.MIRI_CRUCIFORM_INNER_RADIUS_PIX] *= 0.5
     kernel_y[np.abs(y) < constants.MIRI_CRUCIFORM_INNER_RADIUS_PIX] *= 0.5
 
-    # Add in the offset copies of the main 1d kernels
+    # add in the extra diffraction peaks into each 1d kernel, as seen in Gaspar et al. 2021
+    # but first, save the total flux for normalization later
+    normfactor = kernel_x.sum()
+
     # Empirically, the 'center' of the cruciform shifts inwards towards the center of the detector
     # i.e. for the upper right corner, the cruciform shifts down and left a bit, etc. 
+    # Turns out the centers of the peaks in each cruciform arm also seem to shift a bit, so
+    # let's use the same shifts to try to model those, with some scale factor
     yshift = cruciform_yshifts(detector_position[1])
     xshift = cruciform_xshifts(detector_position[0])
+
+    for loc, amp in zip(constants.MIRI_CRUCIFORM_PEAKS_LOC, constants.MIRI_CRUCIFORM_PEAKS_AMP):
+        # Empirically, the locations of the peaks shift slightly around the FOV, in the opposite sign as the cruciform itself shifts
+        # we model this ad hoc based on comparisons iwth the ePSF data.
+        # The scale factors here are a bit of a handwave by eye, not yet a rigorous fit...
+        scaled_loc_x = loc * wavelength / constants.MIRI_CRUCIFORM_PEAK_REFWAVE * 1.1
+        scaled_loc_y = loc * wavelength / constants.MIRI_CRUCIFORM_PEAK_REFWAVE * 1.1
+
+        scaled_amp = amp * constants.MIRI_CRUCIFORM_PEAKS_AMP_ADJUST  # ad hoc, to make it work; basically a units scale factor from the plot in Andras' paper 
+        peak_x = scaled_amp * np.exp(-(x-scaled_loc_x + xshift/2)**2) + scaled_amp * np.exp(-(x+scaled_loc_x + xshift/2)**2)
+        peak_y = scaled_amp * np.exp(-(y-scaled_loc_y + yshift/2)**2) + scaled_amp * np.exp(-(y+scaled_loc_y + yshift/2)**2)
+        kernel_x += peak_x
+        kernel_y += peak_y
+    kernel_x *= normfactor/kernel_x.sum()
+    kernel_y *= normfactor/kernel_y.sum()
+
+    # Add in the offset copies of the main 1d kernels
     kernel_2d[cen + int(round(yshift*oversample))] = kernel_x
     kernel_2d[:, cen + int(round(xshift*oversample))] = kernel_y