diff --git a/dev-docs/source/AbinsImplementation.rst b/dev-docs/source/AbinsImplementation.rst
index 3f83de0c944fa6c1b291014d37cf89c71d44c41b..0fa88a4669a87921075d6d32e7b621541916acb8 100644
--- a/dev-docs/source/AbinsImplementation.rst
+++ b/dev-docs/source/AbinsImplementation.rst
@@ -42,6 +42,13 @@ named _freq_points_.
Broadening
----------
+Instrumental broadening in Abins involves an energy-dependent
+resolution function; in the implemented case (the TOSCA instrument)
+this is convolution with a Gaussian kernel with the energy-dependent
+width parameter (sigma) determined by a quadratic polynomial.
+Convolution with a varying kernel is not a common operation and is
+implemented within AbinsModules.
+
Earlier versions of Abins implemented a Gaussian kernel with a
fixed number of points spread over a range scaled to the peak width,
set by *pkt_per_peak* and *fwhm* parameters in
@@ -63,8 +70,56 @@ is made available so that this can be overruled, but it is up to the
Instrument to interpret this information. 'auto' should select an
intelligent scheme and inappropriate methods can be forbidden.
+Fast approximate broadening with "interpolate"
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The "interpolate" scheme estimates broadened spectra using a limited
+number of kernel evaluations or convolution steps in order to reduce
+the computation time. The method appears to be novel, so some
+explanation is needed.
+
+Consider first that we can approximate a Gaussian function by a linear
+combination of two other Gaussians; one narrower and one wider. If the
+mixing parameters are linearly based on the relationship between the
+widths the results are not impressive:
+
+.. image:: ../images/gaussian_mix_linear.png
+
+But if we optimise the mixing parameter at each width then the
+magnitudes improve significantly, even if the shapes remain distinctly non-Gaussian:
+
+.. image:: ../images/gaussian_mix_optimal_scale4.png
+
+This error is closely related to the width difference between the
+endpoints. Here the range is reduced from a factor 4 to a factor 2,
+and the resulting functions are visually quite convincing
+
+.. image:: ../images/gaussian_mix_optimal_scale2.png
+
+while a gap of :math:`\sqrt{2}` is practically indistinguishable with error below 1% of the peak maximum.
+
+.. image:: ../images/gaussian_mix_optimal_scale_sqrt2.png
+
+For TOSCA :math:`\sigma = a f^2 + b f + c` where :math:`a, b, c$ = $10^{-7}, 0.005, 2.5`. For an energy range of 32 cm\ :sup:`-1` to 4100 cm\ :sup:`-1` sigma ranges from 2.66 to 24.68, which could covered by 5 Gaussians separated by width factor 2 or 9 Gaussians seperated by width factor :math:`\sqrt{2}`.
+This could present a significant cost saving compared to full evaluation of ~4000 convolution kernels (one per convolution bin).
+Alternatively we can perform convolution of the full spectrum with each of the sampled kernels, and then interpolate *between the spectra* using the predetermined mixing weights. The convolution is performed efficiently using FFTs, and relatively little memory is required to hold this limited number of spectra and interpolate between them.
+NEW FIGURE NEEDED: ILLUSTRATE INTERP BETWEEN SPECTRA
+
+This procedure is not strictly equivalent to a summation over frequency-dependent functions, even if there is no interpolation error.
+At each energy coordinate :math:`\epsilon` we "see" a fragment of full spectrum convolved at the same width as any points at :math:`\epsilon` would be.
+In a typical indirect INS spectrum which becomes broader at high energy, this would overestimate the contribution from peaks originating below this :math:`\epsilon` and underestimate the contribution from peaks originating above :math:`\epsilon`.
+As a result, peaks will appear asymmetric.
+In practice, the magnitude of this error depends on the rate of change of :math:`\sigma` relative to the size of :math:`\sigma`.
+In the case of the TOSCA parameters, the error is very small. This should be re-evaluated for other instruments with different energy-dependent broadening functions.
+
+.. image:: ../images/abins-interpolation-benzene.png
+
+We can see the artefacts of this approach more clearly if we use a wider sample spacing (factor 2) and zoom in on the spectrum. The interpolation method has a tendency to show small peaks at turning points; this may be related to the imperfection in the shape of the smooth bell.
+
+.. image:: ../images/abins-interpolation-zoom.png
+
Deprecation plans
-----------------
@@ -80,3 +135,117 @@ Deprecation plans
- The *frequencies_threshold* parameter in
``AbinsParameters.sampling`` is currently non-functional and should
be removed until it *is* functional.
+
+
+
+.. Source code for Gaussian interpolation plots
+
+ import numpy as np
+ from scipy.optimize import curve_fit
+ import matplotlib.pyplot as plt
+ from matplotlib.lines import Line2D
+
+ def gaussian(x, sigma=2, center=0):
+ g = np.exp(-0.5 * ((x - center) / sigma)**2) / (sigma * np.sqrt(2 * np.pi))
+ return g
+
+ margin = 0.05
+
+ def plot_linear_interp(filename='gaussian_mix_linear.png'):
+ """Plot linearly-interpolated Gaussians with wide sigma range"""
+
+ g1_center = 0
+ g2_center = 40
+ sigma_max = 4
+ sigma_min = 1
+
+ x = np.linspace(-10, 50, 401)
+
+ fig, ax = plt.subplots()
+ for sigma, color in zip(np.linspace(sigma_min, sigma_max, 5),
+ ['C0', 'C1', 'C2', 'C3', 'C4']):
+ center = (g1_center
+ + ((sigma - sigma_min)
+ * (g2_center - g1_center) / (sigma_max - sigma_min)))
+ ax.plot(x, gaussian(x, sigma=sigma, center=center), c=color)
+
+ low_ref = gaussian(x, sigma=sigma_min, center=center)
+ high_ref = gaussian(x, sigma=sigma_max, center=center)
+ mix = (sigma - sigma_min) / (sigma_max - sigma_min)
+ ax.plot(x, (1 - mix) * low_ref + mix * high_ref,
+ c=color, linestyle='--')
+
+ ax.set_xlim([-10, 50])
+ ax.set_ylim([0, None])
+ ax.tick_params(labelbottom=False, labelleft=False)
+
+ custom_lines = [Line2D([0], [0], color='k', linestyle='-', lw=2),
+ Line2D([0], [0], color='k', linestyle='--', lw=2)]
+
+ ax.legend(custom_lines, ['Exact', 'Linear interpolation'])
+ fig.subplots_adjust(left=margin, bottom=margin,
+ right=(1 - margin), top=(1 - margin))
+
+ fig.savefig(filename)
+
+ def plot_optimised_interp(filename='gaussian_mix_optimal_scale4.png',
+ sigma_max=4):
+ g1_center = 0
+ g2_center = 40
+ sigma_min = 1
+
+ x = np.linspace(-10, 10, 101)
+ npts = 7
+
+ fig, [ax1, ax2, ax3] = plt.subplots(nrows=3,
+ sharex=True,
+ gridspec_kw={
+ 'height_ratios': [3, 1, 1]})
+ mix1_list, mix2_list = [], []
+
+ def gaussian_mix(x, w1, w2):
+ """Return a linear combination of two Gaussians with weights"""
+ return (w1 * gaussian(x, sigma=sigma_min)
+ + w2 * gaussian(x, sigma=sigma_max))
+
+
+ for sigma, color in zip(np.linspace(sigma_min, sigma_max, npts),
+ ['C0', 'C1', 'C2', 'C3', 'C4', 'C5', 'C6']):
+ ydata = gaussian(x, sigma=sigma)
+ (mix1, mix2), _ = curve_fit(gaussian_mix, x, ydata, p0=[0.5, 0.5])
+
+ x_offset = (g1_center
+ + ((sigma - sigma_min)
+ * (g2_center - g1_center) / (sigma_max - sigma_min)))
+ actual = gaussian(x, sigma=sigma)
+ est = gaussian_mix(x, mix1, mix2)
+ rms = np.sqrt(np.mean((actual - est)**2))
+ ax1.plot(x + x_offset, actual, color=color)
+ ax1.plot(x + x_offset, est, color=color, linestyle='--')
+ ax2.plot([x_offset], [rms], 'o', c='C0')
+
+ mix1_list.append(mix1)
+ mix2_list.append(mix2)
+
+
+ custom_lines = [Line2D([0], [0], color='k', linestyle='-', lw=2),
+ Line2D([0], [0], color='k', linestyle='--', lw=2)]
+
+ ax1.legend(custom_lines, ['Exact', 'Optimised interpolation'])
+
+ ax2.set_ylabel('RMS error')
+
+ ax3.plot(np.linspace(g1_center, g2_center, npts), mix1_list)
+ ax3.plot(np.linspace(g1_center, g2_center, npts), mix2_list)
+ ax3.set_ylabel('Weights')
+ ax3.set_ylim([0, 1])
+
+ fig.savefig(filename)
+
+ if __name__ == '__main__':
+ plot_linear_interp()
+ plot_optimised_interp()
+ plot_optimised_interp(filename='gaussian_mix_optimal_scale2.png',
+ sigma_max=2)
+ plot_optimised_interp(filename='gaussian_mix_optimal_scale_sqrt2.png',
+ sigma_max=np.sqrt(2))
diff --git a/dev-docs/source/images/abins-interpolation-benzene.png b/dev-docs/source/images/abins-interpolation-benzene.png
new file mode 100644
index 0000000000000000000000000000000000000000..48a15807d22a304e4d26d948e91080701e7ae0e1
Binary files /dev/null and b/dev-docs/source/images/abins-interpolation-benzene.png differ
diff --git a/dev-docs/source/images/abins-interpolation-zoom.png b/dev-docs/source/images/abins-interpolation-zoom.png
new file mode 100644
index 0000000000000000000000000000000000000000..77cd2a6ec83035b76046669326bc62a2deb5d13e
Binary files /dev/null and b/dev-docs/source/images/abins-interpolation-zoom.png differ
diff --git a/dev-docs/source/images/gaussian_mix_linear.png b/dev-docs/source/images/gaussian_mix_linear.png
new file mode 100644
index 0000000000000000000000000000000000000000..9041c9431df6a0ea4103a56cd6ef109fc0570664
Binary files /dev/null and b/dev-docs/source/images/gaussian_mix_linear.png differ
diff --git a/dev-docs/source/images/gaussian_mix_optimal_scale2.png b/dev-docs/source/images/gaussian_mix_optimal_scale2.png
new file mode 100644
index 0000000000000000000000000000000000000000..04d291b1ba9384f0dab14f21bdc43dce4d162215
Binary files /dev/null and b/dev-docs/source/images/gaussian_mix_optimal_scale2.png differ
diff --git a/dev-docs/source/images/gaussian_mix_optimal_scale4.png b/dev-docs/source/images/gaussian_mix_optimal_scale4.png
new file mode 100644
index 0000000000000000000000000000000000000000..0080fc734f77ac1bc8da1aa8dd7a7fa4012117c2
Binary files /dev/null and b/dev-docs/source/images/gaussian_mix_optimal_scale4.png differ
diff --git a/dev-docs/source/images/gaussian_mix_optimal_scale_sqrt2.png b/dev-docs/source/images/gaussian_mix_optimal_scale_sqrt2.png
new file mode 100644
index 0000000000000000000000000000000000000000..52eb2afe2b3ad48ecd80a7262d906782497b767e
Binary files /dev/null and b/dev-docs/source/images/gaussian_mix_optimal_scale_sqrt2.png differ