From 067843cd6009f8b7640eb92d256809edad40ce07 Mon Sep 17 00:00:00 2001
From: Morgan Hayward <mh653@leicester.ac.uk>
Date: Wed, 14 May 2025 23:13:09 +0100
Subject: [PATCH] General: Moved peak width at half height values into the
 peaklist variables to allow for peak-wise (and therein multiplet-wise) width
 control. Peaklist handling methods and functions have been adapted to handle
 the new values. Width default values have been consistently set to 0.5Hz
 throughout the package.

plt.py:
Divided mplplot into two functions (mplplot in plt.py and create_lineshape in math.py) to separate the matplotlib object creation from the NMR line. This allows users to skip matplotlib if they want to use the line for other uses.
Fixed a bug in mplplot_stick which uses the outdated use_line_collection argument in ax.stem.

qm.py:
Removed the print lines in qm.hamiltonian_sparse.
Added width arguments to _compile_peaklist permit the new peaklist structure.

math.py:
Added the new create_lineshape function that is taken from plt.mplplot so it can be called independently.
Added width arguments to many of the functions.
Added gaussian functions to complement lorentz functions.
Converted add_lorentzians to add_lorentzians_limitable which includes a new and optional parameter: cutoff. This argument allows the restriction of lorentzian curve calculation to only above the given intensity cutoff. Uses the (also new) inverse_lorentz and get_index_of_nearest functions to quickly calculate the restriction zones without any element-wise calculation by numpy. This permits much faster simulation speeds if speed is more valuable than 100% accuracy. Default is None so tutorial/older examples behave the same.

Other files:
Added support for the new width parameter.
Refactoring for add_lorentzians_limitable function.

Tests:
All executed successfully after modifying expected values to accommodate the above changes.
---
 src/nmrsim/_classes.py   |  14 +--
 src/nmrsim/discrete.py   |  40 ++++----
 src/nmrsim/firstorder.py |  18 ++--
 src/nmrsim/math.py       | 192 +++++++++++++++++++++++++++++++++++----
 src/nmrsim/plt.py        |  31 +++----
 src/nmrsim/qm.py         |  11 +--
 tests/accepted_data.py   |  24 ++---
 tests/test_firstorder.py |  34 ++++---
 tests/test_math.py       |  83 ++++++++++++-----
 tests/test_nmrsim.py     |  26 +++---
 tests/test_plt.py        |  22 +++--
 11 files changed, 350 insertions(+), 145 deletions(-)

diff --git a/src/nmrsim/_classes.py b/src/nmrsim/_classes.py
index a02088d..233e2b1 100644
--- a/src/nmrsim/_classes.py
+++ b/src/nmrsim/_classes.py
@@ -7,7 +7,7 @@
 import numpy as np
 
 from nmrsim.firstorder import first_order_spin_system, multiplet
-from nmrsim.math import reduce_peaks, add_lorentzians
+from nmrsim.math import reduce_peaks, add_lorentzians_limitable
 from nmrsim.qm import qm_spinsystem
 from nmrsim._utils import low_high
 
@@ -73,7 +73,7 @@ def __init__(self, v, I, J, w=0.5):
         self.I = I
         self.J = J
         self.w = w
-        self._peaklist = multiplet((v, I), J)
+        self._peaklist = multiplet((v, I, w), J)
 
     def __eq__(self, other):
         if hasattr(other, "peaklist") and callable(other.peaklist):
@@ -87,7 +87,7 @@ def __add__(self, other):
 
     def __mul__(self, scalar):
         if isinstance(scalar, numbers.Real):
-            return Multiplet(self.v, self.I * scalar, self.J)
+            return Multiplet(self.v, self.I * scalar, self.J, self.w)
         else:
             return NotImplemented
 
@@ -105,7 +105,7 @@ def __itruediv__(self, scalar):
         return self.__imul__(1 / scalar)
 
     def _refresh(self):
-        self._peaklist = multiplet((self.v, self.I), self.J)
+        self._peaklist = multiplet((self.v, self.I, self.w), self.J)
 
     def peaklist(self):
         """
@@ -245,9 +245,9 @@ def peaklist(self):
             Array of (frequency, intensity) signals.
         """
         if self._second_order:
-            return qm_spinsystem(self._v, self._J)
+            return qm_spinsystem(self._v, self._J, width=self.w)
         else:
-            return first_order_spin_system(self._v, self._J)
+            return first_order_spin_system(self._v, self._J, self.w)
 
     def __eq__(self, other):
         if hasattr(other, "peaklist") and callable(other.peaklist):
@@ -368,7 +368,7 @@ def lineshape(self, points=800):
         """
         vmin, vmax = low_high((self.vmin, self.vmax))
         x = np.linspace(vmin, vmax, points)
-        y = [add_lorentzians(x, c.peaklist(), c.w) for c in self._components]
+        y = [add_lorentzians_limitable(x, c.peaklist()) for c in self._components]
         y_sum = np.sum(y, 0)
         return x, y_sum
 
diff --git a/src/nmrsim/discrete.py b/src/nmrsim/discrete.py
index eb1c6fe..a0d785a 100644
--- a/src/nmrsim/discrete.py
+++ b/src/nmrsim/discrete.py
@@ -55,8 +55,8 @@ def AB(Jab, Vab, Vcentr, normalize=True):
 
     Returns
     -------
-    [(float, float)...]
-        A list of four (frequency, intensity) tuples.
+    [(float, float, float)...]
+        A list of four (frequency, intensity, width) tuples.
     """
     J = Jab
     dv = Vab
@@ -73,9 +73,10 @@ def AB(Jab, Vab, Vcentr, normalize=True):
     I4 = I1
     vList = [v1, v2, v3, v4]
     IList = [I1, I2, I3, I4]
+    wList = [0.5 for _ in range(len(IList))]
     if normalize:
         _normalize(IList, 2)
-    return list(zip(vList, IList))
+    return list(zip(vList, IList, wList))
 
 
 def AB2(Jab, Vab, Vcentr, normalize=True):
@@ -96,8 +97,8 @@ def AB2(Jab, Vab, Vcentr, normalize=True):
 
     Returns
     -------
-    [(float, float)...]
-        a list of (frequency, intensity) tuples.
+    [(float, float, float)...]
+        a list of (frequency, intensity, width) tuples.
     """
     # There is a disconnect between the variable names in the WINDNMR GUI and
     # the variable names in this function.
@@ -148,10 +149,11 @@ def AB2(Jab, Vab, Vcentr, normalize=True):
     I9 = (sqrt(2) * sin_dtheta + sintheta_plus * sintheta_minus) ** 2
     vList = [V1, V2, V3, V4, V5, V6, V7, V8, V9]
     IList = [I1, I2, I3, I4, I5, I6, I7, I8, I9]
+    wList = [0.5 for _ in range(len(IList))]
 
     if normalize:
         _normalize(IList, 3)
-    return list(zip(vList, IList))
+    return list(zip(vList, IList, wList))
 
 
 def ABX(Jab, Jax, Jbx, Vab, Vcentr, vx, normalize=True):
@@ -181,8 +183,8 @@ def ABX(Jab, Jax, Jbx, Vab, Vcentr, vx, normalize=True):
 
     Returns
     -------
-    [(float, float)...]
-        a list of (frequency, intensity) tuples.
+    [(float, float, float)...]
+        a list of (frequency, intensity, width) tuples.
     """
     # Contradictions in naming between WINDNMR's interface, internal code, and
     # Pople/Schneider/Bernstein "fixed" with these reassignments:
@@ -239,9 +241,10 @@ def ABX(Jab, Jax, Jbx, Vab, Vcentr, vx, normalize=True):
     I14 = I13
     VList = [V1, V2, V3, V4, V5, V6, V7, V8, V9, V10, V11, V12, V13, V14]
     IList = [I1, I2, I3, I4, I5, I6, I7, I8, I9, I10, I11, I12, I13, I14]
+    wList = [0.5 for _ in range(len(IList))]
     if normalize:
         _normalize(IList, 3)
-    return list(zip(VList, IList))
+    return list(zip(VList, IList, wList))
 
 
 def ABX3(Jab, Jax, Jbx, Vab, Vcentr):
@@ -264,14 +267,14 @@ def ABX3(Jab, Jax, Jbx, Vab, Vcentr):
 
     Returns
     -------
-    [(float, float)...]
-        a list of (frequency, intensity) tuples.
+    [(float, float, float)...]
+        a list of (frequency, intensity, width) tuples.
     """
     # First: simulate two quartets for va and vb ("Jab turned off")
     va = Vcentr - Vab / 2
     vb = Vcentr + Vab / 2
-    a_quartet = multiplet((va, 1), [(Jax, 3)])
-    b_quartet = multiplet((vb, 1), [(Jbx, 3)])
+    a_quartet = multiplet((va, 1, 0.5), [(Jax, 3)])
+    b_quartet = multiplet((vb, 1, 0.5), [(Jbx, 3)])
     res = []
     # Then: for each pair of a and b singlets in the quartets, calculate an
     # AB quartet ("Turn Jab on").
@@ -280,7 +283,7 @@ def ABX3(Jab, Jax, Jbx, Vab, Vcentr):
         abcenter = (b_quartet[i][0] + a_quartet[i][0]) / 2
         sub_abq = AB(Jab, dv, abcenter, normalize=True)
         scale_factor = a_quartet[i][1]
-        scaled_sub_abq = [(v, i * scale_factor) for v, i in sub_abq]
+        scaled_sub_abq = [(v, i * scale_factor, w) for v, i, w in sub_abq]
         res.extend(scaled_sub_abq)
     return res
 
@@ -305,8 +308,8 @@ def AAXX(Jaa, Jxx, Jax, Jax_prime, Vcentr, normalize=True):
 
     Returns
     -------
-    [(float, float)...]
-        a list of (frequency, intensity) tuples.
+    [(float, float, float)...]
+        a list of (frequency, intensity, width) tuples.
     """
     # Define the constants required to calculate frequencies and intensities
 
@@ -347,9 +350,10 @@ def AAXX(Jaa, Jxx, Jax, Jax_prime, Vcentr, normalize=True):
 
     VList = [V1, V2, V3, V4, V5, V6, V7, V8, V9, V10]
     IList = [I1, I2, I3, I4, I5, I6, I7, I8, I9, I10]
+    wList = [0.5 for _ in range(len(IList))]
     if normalize:
         _normalize(IList, 2)
-    return list(zip(VList, IList))
+    return list(zip(VList, IList, wList))
 
 
 def AABB(Vab, Jaa, Jbb, Jab, Jab_prime, Vcentr, normalize=True, **kwargs):
@@ -374,7 +378,7 @@ def AABB(Vab, Jaa, Jbb, Jab, Jab_prime, Vcentr, normalize=True, **kwargs):
 
     Returns
     -------
-    [(float, float)...]
+    [(float, float, float)...]
         a list of (frequency, intensity) tuples.
     """
     from nmrsim.qm import qm_spinsystem
diff --git a/src/nmrsim/firstorder.py b/src/nmrsim/firstorder.py
index 414df9a..bb821f0 100644
--- a/src/nmrsim/firstorder.py
+++ b/src/nmrsim/firstorder.py
@@ -30,9 +30,9 @@ def _doublet(plist, J):
         a list of (frequency, intensity) tuples.
     """
     res = []
-    for v, i in plist:
-        res.append((v - J / 2, i / 2))
-        res.append((v + J / 2, i / 2))
+    for v, i, w in plist:
+        res.append((v - J / 2, i / 2, w))
+        res.append((v + J / 2, i / 2, w))
     return res
 
 
@@ -43,8 +43,8 @@ def multiplet(signal, couplings):
 
     Parameters
     ---------
-    signal : (float, float)
-        a (frequency (Hz), intensity) tuple;
+    signal : (float, float, float)
+        a (frequency (Hz), intensity, width (Hz)) tuple;
     couplings : [(float, int)...]
         A list of (*J*, # of nuclei) tuples. The order of the tuples in
         couplings does not matter.
@@ -64,7 +64,7 @@ def multiplet(signal, couplings):
     return sorted(reduce_peaks(res))
 
 
-def first_order_spin_system(v, J):
+def first_order_spin_system(v, J, w):
     """
     Create a first-order peaklist of several multiplets from the same v/J
     arguments used for qm calculations.
@@ -80,16 +80,18 @@ def first_order_spin_system(v, J):
         an array of frequencies
     J : 2D array-like (square)
         a matrix of J coupling constants
+    w : float or int
+        peak width at half height (Hz)
 
     Returns
     -------
-    [(float, float)...]
+    [(float, float, float)...]
         a combined peaklist of signals for all the multiplets in the spin
         system.
     """
     result = []
     for i, v_ in enumerate(v):
         couplings = ((j, 1) for j in J[i] if j != 0)
-        signal = multiplet((v_, 1), couplings)
+        signal = multiplet((v_, 1, w), couplings)
         result += signal
     return reduce_peaks(sorted(result))
diff --git a/src/nmrsim/math.py b/src/nmrsim/math.py
index 90883b0..d1511a5 100644
--- a/src/nmrsim/math.py
+++ b/src/nmrsim/math.py
@@ -20,7 +20,51 @@
 * get_intensity: given a lineshape and a frequency, find the intensity at the
   datapoint closest to that frequency.
 """
+import math
 import numpy as np
+from nmrsim._utils import low_high
+
+
+def create_lineshape(peaklist, points=800, limits=None, function='lorentzian', cutoff=None):
+    """
+    Numpy Arrays of the simulated lineshape for a peaklist.
+    Parameters
+    ----------
+    peaklist : [(float, float, float)...]
+        A list of (frequency, intensity, width) tuples.
+    y_min : float or int
+        Minimum intensity for the plot.
+    y_max : float or int
+        Maximum intensity for the plot.
+    points : int
+        Number of data points.
+    limits : (float, float)
+        Frequency limits for the plot.
+    function: string
+        Plotting function for the peak shape, either lorentzian or gaussian.
+    cutoff: float
+        The minimum height to calculate for the peak curves.
+    Returns
+    -------
+    x, y : numpy.array
+        Arrays for frequency (x) and intensity (y) for the simulated lineshape.
+    """
+    peaklist.sort()
+    if limits:
+        l_limit, r_limit = low_high(limits)
+    else:
+        l_limit = peaklist[0][0] - 50
+        r_limit = peaklist[-1][0] + 50
+    x_inv = np.linspace(l_limit, r_limit, points)
+    x = x_inv[::-1]  # reverses the x axis.
+    if len(peaklist) > 0:
+        if function == 'lorentzian':
+            y = add_lorentzians_limitable(x, peaklist, cutoff)
+        elif function == 'gaussian':
+            y = add_gaussians(x, peaklist)
+    else:
+        y = np.zeros(points)
+    return x, y
 
 
 def add_peaks(plist):
@@ -40,10 +84,10 @@ def add_peaks(plist):
     """
     v_total = 0
     i_total = 0
-    for v, i in plist:
+    for v, i, w in plist:
         v_total += v
         i_total += i
-    return v_total / len(plist), i_total
+    return v_total / len(plist), i_total, w
 
 
 def reduce_peaks(plist_, tolerance=0):
@@ -110,9 +154,9 @@ def normalize_peaklist(peaklist, n=1):
     n : int or float (optional)
         total intensity to normalize to (default = 1).
     """
-    freq, int_ = [x for x, y in peaklist], [y for x, y in peaklist]
+    freq, int_, width = [x for x, y, w in peaklist], [y for x, y, w in peaklist], [w for x, y, w in peaklist]
     _normalize(int_, n)
-    return list(zip(freq, int_))
+    return list(zip(freq, int_, width))
 
 
 def lorentz(v, v0, I, w):
@@ -146,32 +190,148 @@ def lorentz(v, v0, I, w):
     return scaling_factor * I * ((0.5 * w) ** 2 / ((0.5 * w) ** 2 + (v - v0) ** 2))
 
 
-def add_lorentzians(linspace, peaklist, w):
+def inverse_lorentz(v0, I, w, cutoff):
+    """
+    A rearranged version of the lorentz function.
+
+    Returns two values of frequency where intensity is equal to the cutoff.
+
+    Arguments
+    ---------
+    v0 : float
+        The center of the distribution, aka the peak chemical shift.
+    I : float
+        the relative intensity of the signal, aka the peak height.
+    w : float
+        the peak width at half maximum intensity
+    cutoff: float
+        The intensity value to determine what frequency values are calculated.
+        (WARNING: if 0.0 >= cutoff >= I, it will return complex numbers)
+
+    Returns
+    -------
+    float, float
+        The two frequency values where the intensity equals cutoff.
+    """
+    val = ((((I * ((0.5 * w)**2)) / cutoff) - ((0.5 * w)**2))**0.5)
+    return v0 + val, v0 - val
+
+
+def get_index_of_nearest(linspace, v):
     """
-    Given a numpy linspace, a peaklist of (frequency, intensity)
-    tuples, and a linewidth, returns an array of y coordinates for the
-    total line shape.
+    Quickly returns the nearest index of an evenly distributed array.
+    Does not need to calculate element-wise so is much faster than .argmin() approaches.
+    Can not return any negative values or values greater than the length of the array.
 
     Arguments
     ---------
     linspace : array-like
         Normally a numpy.linspace of x coordinates corresponding to frequency
         in Hz.
-    peaklist : [(float, float)...]
-        A list of (frequency, intensity) tuples.
+    v: float
+        The frequency (x coordinate) in Hz at which to evaluate the index.
+
+    Returns
+    -------
+    int
+        The index of the closest value to v in linspace.
+    """
+    if v >= max(linspace):
+        return 0
+    elif v <= min(linspace):
+        return len(linspace)
+    else:
+        return len(linspace) - round((v - min(linspace)) * (len(linspace) - 1) / (max(linspace) - min(linspace))) - 1
+
+
+def add_lorentzians_limitable(linspace, peaklist, cutoff=None):
+    """
+        Given a numpy linspace, a peaklist of (frequency, intensity, width)
+        tuples and a cutoff for the y..
+        Returns an array of y coordinates for all points > cutoff, or zeros otherwise.
+
+        Arguments
+        ---------
+        linspace : array-like
+            Normally a numpy.linspace of x coordinates corresponding to frequency
+            in Hz.
+        peaklist : [(float, float, float)...]
+            A list of (frequency, intensity, width at half maximum intensity) tuples. Frequency and width in Hz.
+        cutoff : float
+            The intensity value, below which no values are calculated.
+
+        Returns
+        -------
+        [float...]
+            an array of y coordinates corresponding to intensity, if intensity > cutoff, otherwise zero.
+        """
+    result = np.zeros(len(linspace))
+    for v, i, w in peaklist:
+        if cutoff is None:
+            high_index = 0
+            low_index = len(linspace)
+        elif cutoff <= 0:
+            raise ValueError(f"Cutoff must be greater than 0. Got {cutoff}")
+        elif i <= cutoff:
+            continue
+        else:
+            high, low = inverse_lorentz(v, i, w, cutoff)
+            high_index = get_index_of_nearest(linspace, high)
+            low_index = get_index_of_nearest(linspace, low)
+        result[high_index:low_index] += lorentz(linspace[high_index:low_index], v, i, w)
+    return result
+
+
+def gauss(v, v0, I, w):
+    """
+    A gaussian function that takes linewidth at half intensity (w) as a
+    parameter.
+
+    When `v` = `v0`, and `w` = 0.5 (Hz), the function returns intensity I.
+
+    Arguments
+    ---------
+    v : float
+        The frequency (x coordinate) in Hz at which to evaluate intensity (y
+        coordinate).
+    v0 : float
+        The center of the distribution.
+    I : float
+        the relative intensity of the signal
     w : float
-        Peak width at half maximum intensity.
+        the peak width at half maximum intensity
+
+    Returns
+    -------
+    float
+        the intensity (y coordinate) for the Lorentzian distribution
+        evaluated at frequency `v`.
+    """
+    wf = 0.4246609
+    return I * (math.e**((-(v - v0)**2) / (2 * ((w * wf)**2))))
+
+
+def add_gaussians(linspace, peaklist):
+    """
+    Given a numpy linspace, a peaklist of (frequency, intensity, width) tuples.
+    Returns an array of y coordinates using gaussian line calculator.
+
+    Arguments
+    ---------
+    linspace : array-like
+        Normally a numpy.linspace of x coordinates corresponding to frequency
+        in Hz.
+    peaklist : [(float, float, float)...]
+        A list of (frequency, intensity, width at half maximum intensity) tuples. Frequency and width in Hz.
 
     Returns
     -------
     [float...]
         an array of y coordinates corresponding to intensity.
     """
-    # TODO: consider naming, and confusion with .math.add_peaks
-    # TODO: function looks clunky. Refactor?
-    result = lorentz(linspace, peaklist[0][0], peaklist[0][1], w)
-    for v, i in peaklist[1:]:
-        result += lorentz(linspace, v, i, w)
+    result = gauss(linspace, peaklist[0][0], peaklist[0][1], peaklist[0][2])
+    for v, i, w in peaklist[1:]:
+        result += gauss(linspace, v, i, w)
     return result
 
 
diff --git a/src/nmrsim/plt.py b/src/nmrsim/plt.py
index db02c10..165329e 100644
--- a/src/nmrsim/plt.py
+++ b/src/nmrsim/plt.py
@@ -14,7 +14,7 @@
 """
 import numpy as np
 
-from nmrsim.math import add_lorentzians
+from nmrsim.math import create_lineshape
 from nmrsim._utils import low_high
 
 # Pyplot assumes a TkAgg backend as a default. This can cause problems in
@@ -37,16 +37,14 @@
 # TODO: possibly refactor plot routines to avoid repetitive code
 
 
-def mplplot(peaklist, w=1, y_min=-0.01, y_max=1, points=800, limits=None, hidden=False):
+def mplplot(peaklist, y_min=-0.01, y_max=1, points=800, limits=None, hidden=False):
     """
     A matplotlib plot of the simulated lineshape for a peaklist.
 
     Parameters
     ----------
-    peaklist : [(float, float)...]
-        A list of (frequency, intensity) tuples.
-    w : float
-        Peak width at half height
+    peaklist : [(float, float, float)...]
+        A list of (frequency, intensity, width) tuples.
     y_min : float or int
         Minimum intensity for the plot.
     y_max : float or int
@@ -63,18 +61,14 @@ def mplplot(peaklist, w=1, y_min=-0.01, y_max=1, points=800, limits=None, hidden
     x, y : numpy.array
         Arrays for frequency (x) and intensity (y) for the simulated lineshape.
     """
-    peaklist.sort()
-    if limits:
-        l_limit, r_limit = low_high(limits)
-    else:
-        l_limit = peaklist[0][0] - 50
-        r_limit = peaklist[-1][0] + 50
-    x = np.linspace(l_limit, r_limit, points)
-    plt.ylim(y_min, y_max)
-    plt.gca().invert_xaxis()  # reverses the x axis
-    y = add_lorentzians(x, peaklist, w)
+    # Enforce peaklist has 3 values.
+    # Allows uses with the previous peaklist parameter of only height/width.
+    peaklist = [(peak[0], peak[1], peak[2]) if len(peak) == 3 else (peak[0], peak[1], 0.5) for peak in peaklist]
+    x, y = create_lineshape(peaklist, points=points, limits=limits, function='lorentzian')
     # noinspection PyTypeChecker
     lines = plt.plot(x, y)
+    plt.ylim(y_min, y_max)
+    plt.gca().invert_xaxis()  # reverses the x axis
     print(lines)
     if not hidden:
         plt.show()
@@ -102,6 +96,9 @@ def mplplot_stick(peaklist, y_min=-0.01, y_max=1, limits=None, hidden=False):
     numpy.array, numpy.array
         The arrays of x and y coordinates used for the plot.
     """
+    # Enforce peaklist has 2 values.
+    # Allows users to provide the same peaklist parameter as mplplot without breaking.
+    peaklist = [(peak[0], peak[1]) if len(peak) == 3 else peak for peak in peaklist]
     fig, ax = plt.subplots()
     if limits:
         l_limit, r_limit = low_high(limits)
@@ -117,7 +114,7 @@ def mplplot_stick(peaklist, y_min=-0.01, y_max=1, limits=None, hidden=False):
     y = np.append(y, [0.001, 0.001])
     plt.xlim(r_limit, l_limit)
     plt.ylim(y_min, y_max)
-    ax.stem(x, y, markerfmt=" ", basefmt="C0-", use_line_collection=True)  # suppress warning until mpl 3.3
+    ax.stem(x, y, markerfmt=" ", basefmt="C0-")  # suppress warning until mpl 3.3
     if not hidden:
         plt.show()
     return x, y
diff --git a/src/nmrsim/qm.py b/src/nmrsim/qm.py
index 5cff0a4..eecd243 100644
--- a/src/nmrsim/qm.py
+++ b/src/nmrsim/qm.py
@@ -222,10 +222,6 @@ def hamiltonian_sparse(v, J):
     """
     nspins = len(v)
     Lz, Lproduct = _so_sparse(nspins)  # noqa
-    # TODO: remove the following lines once tests pass
-    print("From hamiltonian_sparse:")
-    print("Lz is type: ", type(Lz))
-    print("Lproduct is type: ", type(Lproduct))
     assert isinstance(Lz, (sparse.COO, np.ndarray, scipy.sparse.spmatrix))
     # On large spin systems, converting v and J to sparse improved speed of
     # sparse.tensordot calls with them.
@@ -393,7 +389,7 @@ def _intensity_and_energy(H, nspins):
     return I, E
 
 
-def _compile_peaklist(I, E, cutoff=0.001):
+def _compile_peaklist(I, E, cutoff=0.001, width=0.5):
     """
     Generate a peaklist from intensity and energy matrices.
 
@@ -414,8 +410,9 @@ def _compile_peaklist(I, E, cutoff=0.001):
     I_upper = np.triu(I)
     E_matrix = np.abs(E[:, np.newaxis] - E)
     E_upper = np.triu(E_matrix)
-    combo = np.stack([E_upper, I_upper])
-    iv = combo.reshape(2, I.shape[0] ** 2).T
+    W = np.full(I_upper.shape, width)
+    combo = np.stack([E_upper, I_upper, W])
+    iv = combo.reshape(3, I.shape[0] ** 2).T
     return iv[iv[:, 1] >= cutoff]
 
 
diff --git a/tests/accepted_data.py b/tests/accepted_data.py
index 8fe7975..33f210b 100644
--- a/tests/accepted_data.py
+++ b/tests/accepted_data.py
@@ -8,18 +8,18 @@
     [0., 0., 0., 7.5, 0., 0.75, -199.875, 0.],
     [0., 0., 0., 0., 0., 0., 0., -491.625]]
 
-SPECTRUM_RIOUX = [(260.66152857482973, 0.23011249131787825),
-                  (291.3191136690316, 0.22882003310401844),
-                  (419.5193577561387, 0.2910724454555947),
-                  (292.8468885409387, 0.21381231157251776),
-                  (426.4877446901903, 0.2662986769673388),
-                  (262.1893034467369, 0.2487606172641341),
-                  (434.5231959501799, 0.23004586196807364),
-                  (267.62991550888137, 0.24855578963215685),
-                  (306.3229518630728, 0.2925168028407963),
-                  (441.49158288423155, 0.2125727818192929),
-                  (307.85072673497996, 0.2648680204713063),
-                  (269.1576903807885, 0.27256416758689167)]
+SPECTRUM_RIOUX = [(260.66152857482973, 0.23011249131787825, 0.5),
+                  (291.3191136690316, 0.22882003310401844, 0.5),
+                  (419.5193577561387, 0.2910724454555947, 0.5),
+                  (292.8468885409387, 0.21381231157251776, 0.5),
+                  (426.4877446901903, 0.2662986769673388, 0.5),
+                  (262.1893034467369, 0.2487606172641341, 0.5),
+                  (434.5231959501799, 0.23004586196807364, 0.5),
+                  (267.62991550888137, 0.24855578963215685, 0.5),
+                  (306.3229518630728, 0.2925168028407963, 0.5),
+                  (441.49158288423155, 0.2125727818192929, 0.5),
+                  (307.85072673497996, 0.2648680204713063, 0.5),
+                  (269.1576903807885, 0.27256416758689167, 0.5)]
 
 # accepted OLD output [(x, y)...] for nmrsim.plt.add_lorentzians, using:
 #  linspace = np.linspace(390, 410, 200)
diff --git a/tests/test_firstorder.py b/tests/test_firstorder.py
index 86fe471..e128fcb 100644
--- a/tests/test_firstorder.py
+++ b/tests/test_firstorder.py
@@ -6,22 +6,30 @@
 
 
 def test_multiplet_allows_singlet():
-    refspec = [(1200.0, 2.0)]
+    refspec = [(1200.0, 2.0, 0.5)]
     # GIVEN a signal
     # WHEN multiplet is called with an empty J list
-    testspec = multiplet((1200.0, 2.0), [])
+    testspec = multiplet((1200.0, 2.0, 0.5), [])
     # THEN the returned peaklist only contains the original signal
     assert np.allclose(refspec, testspec)
 
 
 def test_multiplet():
-    refspec = [(293.0, 0.75), (300.0, 1.5), (307.0, 0.75),
-               (432.5, 0.0625), (439.5, 0.3125), (446.5, 0.625),
-               (453.5, 0.625), (460.5, 0.3125), (467.5, 0.0625),
-               (1193.0, 0.5), (1200.0, 1.0), (1207.0, 0.5)]
-    v1 = (1200, 2)
-    v2 = (450, 2)
-    v3 = (300, 3)
+    refspec = [(293.0, 0.75, 1.0),
+               (300.0, 1.5, 1.0),
+               (307.0, 0.75, 1.0),
+               (432.5, 0.0625, 0.5),
+               (439.5, 0.3125, 0.5),
+               (446.5, 0.625, 0.5),
+               (453.5, 0.625, 0.5),
+               (460.5, 0.3125, 0.5),
+               (467.5, 0.0625, 0.5),
+               (1193.0, 0.5, 0.5),
+               (1200.0, 1.0, 0.5),
+               (1207.0, 0.5, 0.5)]
+    v1 = (1200, 2, 0.5)
+    v2 = (450, 2, 0.5)
+    v3 = (300, 3, 1.0)
     J12 = 7
     J23 = 7
     m1 = multiplet(v1, [(J12, 2)])
@@ -33,10 +41,10 @@ def test_multiplet():
 
 def test_first_order_spin_system():
     v, J = rioux()
-    spectrum = first_order_spin_system(v, J)
-    m1 = multiplet((430, 1), [(7, 1), (15, 1)])
-    m2 = multiplet((265, 1), [(7, 1), (1.5, 1)])
-    m3 = multiplet((300, 1), [(15, 1), (1.5, 1)])
+    spectrum = first_order_spin_system(v, J, 0.5)
+    m1 = multiplet((430, 1, 0.5), [(7, 1), (15, 1)])
+    m2 = multiplet((265, 1, 0.5), [(7, 1), (1.5, 1)])
+    m3 = multiplet((300, 1, 0.5), [(15, 1), (1.5, 1)])
     m = reduce_peaks(sorted(m1 + m2 + m3))
     # x = np.array([i[0] for i in spectrum])
     # y = np.array([i[1] for i in spectrum])
diff --git a/tests/test_math.py b/tests/test_math.py
index dac1670..02828d9 100644
--- a/tests/test_math.py
+++ b/tests/test_math.py
@@ -7,32 +7,65 @@
 
 
 def test_add_peaks():
-    peaklist = [(100, 1.1), (110, 1.2), (150, 1.6)]
-    expected_result = (120.0, 3.9)
+    peaklist = [(100, 1.1, 0.5), (110, 1.2, 0.5), (150, 1.6, 0.5)]
+    expected_result = (120.0, 3.9, 0.5)
     result = add_peaks(peaklist)
     assert result == expected_result
 
 
 def test_reduce_peaks():
-    refspec = [(293.0, 0.75), (300.0, 1.5), (307.0, 0.75),
-               (432.5, 0.0625), (439.5, 0.3125), (446.5, 0.625),
-               (453.5, 0.625), (460.5, 0.3125), (467.5, 0.0625),
-               (1193.0, 0.5), (1200.0, 1.0), (1207.0, 0.5)]
-    tobereduced = [
-        (1193.0, 0.5), (1200.0, 0.5), (1200.0, 0.5), (1207.0, 0.5),
-        (432.5, 0.0625), (439.5, 0.0625), (439.5, 0.0625),
-        (446.5, 0.0625), (439.5, 0.0625), (446.5, 0.0625),
-        (446.5, 0.0625), (453.5, 0.0625), (439.5, 0.0625),
-        (446.5, 0.0625), (446.5, 0.0625), (453.5, 0.0625),
-        (446.5, 0.0625), (453.5, 0.0625), (453.5, 0.0625),
-        (460.5, 0.0625), (439.5, 0.0625), (446.5, 0.0625),
-        (446.5, 0.0625), (453.5, 0.0625), (446.5, 0.0625),
-        (453.5, 0.0625), (453.5, 0.0625), (460.5, 0.0625),
-        (446.5, 0.0625), (453.5, 0.0625), (453.5, 0.0625),
-        (460.5, 0.0625), (453.5, 0.0625), (460.5, 0.0625),
-        (460.5, 0.0625), (467.5, 0.0625),
-        (293.0, 0.75), (300.0, 0.75), (300.0, 0.75), (307.0, 0.75)
-    ]
+    refspec = [(293.0, 0.75, 0.5),
+               (300.0, 1.5, 0.5),
+               (307.0, 0.75, 0.5),
+               (432.5, 0.0625, 0.5),
+               (439.5, 0.3125, 0.5),
+               (446.5, 0.625, 0.5),
+               (453.5, 0.625, 0.5),
+               (460.5, 0.3125, 0.5),
+               (467.5, 0.0625, 0.5),
+               (1193.0, 0.5, 0.5),
+               (1200.0, 1.0, 0.5),
+               (1207.0, 0.5, 0.5)]
+    tobereduced = [(1193.0, 0.5, 0.5),
+                   (1200.0, 0.5, 0.5),
+                   (1200.0, 0.5, 0.5),
+                   (1207.0, 0.5, 0.5),
+                   (432.5, 0.0625, 0.5),
+                   (439.5, 0.0625, 0.5),
+                   (439.5, 0.0625, 0.5),
+                   (446.5, 0.0625, 0.5),
+                   (439.5, 0.0625, 0.5),
+                   (446.5, 0.0625, 0.5),
+                   (446.5, 0.0625, 0.5),
+                   (453.5, 0.0625, 0.5),
+                   (439.5, 0.0625, 0.5),
+                   (446.5, 0.0625, 0.5),
+                   (446.5, 0.0625, 0.5),
+                   (453.5, 0.0625, 0.5),
+                   (446.5, 0.0625, 0.5),
+                   (453.5, 0.0625, 0.5),
+                   (453.5, 0.0625, 0.5),
+                   (460.5, 0.0625, 0.5),
+                   (439.5, 0.0625, 0.5),
+                   (446.5, 0.0625, 0.5),
+                   (446.5, 0.0625, 0.5),
+                   (453.5, 0.0625, 0.5),
+                   (446.5, 0.0625, 0.5),
+                   (453.5, 0.0625, 0.5),
+                   (453.5, 0.0625, 0.5),
+                   (460.5, 0.0625, 0.5),
+                   (446.5, 0.0625, 0.5),
+                   (453.5, 0.0625, 0.5),
+                   (453.5, 0.0625, 0.5),
+                   (460.5, 0.0625, 0.5),
+                   (453.5, 0.0625, 0.5),
+                   (460.5, 0.0625, 0.5),
+                   (460.5, 0.0625, 0.5),
+                   (467.5, 0.0625, 0.5),
+                   (293.0, 0.75, 0.5),
+                   (300.0, 0.75, 0.5),
+                   (300.0, 0.75, 0.5),
+                   (307.0, 0.75, 0.5)]
     testspec = reduce_peaks(tobereduced)
     np.testing.assert_array_almost_equal(testspec, refspec, decimal=2)
 
@@ -47,8 +80,8 @@ def test_normalize():
 
 
 def test_normalize_spectrum():
-    unnormalized = [(1200.0, 1.0), (500.0, 2.0)]
-    expected = [(1200.0, 2.0), (500.0, 4.0)]
+    unnormalized = [(1200.0, 1.0, 0.5), (500.0, 2.0, 0.5)]
+    expected = [(1200.0, 2.0, 0.5), (500.0, 4.0, 0.5)]
     result = normalize_peaklist(unnormalized, n=6)
     assert np.allclose(result, expected)
 
@@ -59,8 +92,8 @@ def test_lorentz_width():
     I = 1
     w = 2
     max_height = lorentz(v0, v0, I, w)
-    low_width_height = lorentz(v0 - w/2, v0, I, w)  # noqa: E226
-    high_width_height = lorentz(v0 + w/2, v0, I, w)  # noqa: E226
+    low_width_height = lorentz(v0 - w / 2, v0, I, w)  # noqa: E226
+    high_width_height = lorentz(v0 + w / 2, v0, I, w)  # noqa: E226
     assert low_width_height / max_height == approx(0.5)
     assert high_width_height / max_height == approx(0.5)
 
diff --git a/tests/test_nmrsim.py b/tests/test_nmrsim.py
index c404f31..0b2464f 100644
--- a/tests/test_nmrsim.py
+++ b/tests/test_nmrsim.py
@@ -4,7 +4,7 @@
 from nmrsim import Multiplet, SpinSystem, Spectrum
 from nmrsim._classes import extract_components
 from nmrsim.firstorder import first_order_spin_system
-from nmrsim.math import add_lorentzians
+from nmrsim.math import add_lorentzians_limitable
 from tests.accepted_data import SPECTRUM_RIOUX
 from tests.qm_arguments import rioux
 
@@ -239,8 +239,8 @@ def test_Spectrum_instantiates_with_multiplet(self):
         m1 = Multiplet(100, 1, [(10, 2)])
         m2 = Multiplet(80, 1, [(10, 2)])
         s = Spectrum([m1, m2])
-        expected_peaklist = sorted([(110, 0.25), (100, 0.5), (90, 0.5), (80, 0.5),
-                                    (70, 0.25)])
+        expected_peaklist = sorted([(110, 0.25, 0.5), (100, 0.5, 0.5), (90, 0.5, 0.5), (80, 0.5, 0.5),
+                                    (70, 0.25, 0.5)])
         result = s.peaklist()
         assert np.array_equal(expected_peaklist, result)
 
@@ -251,15 +251,15 @@ def test_add_and_eq(self):
         s2 = m1 + m2
         assert s == s2
         s3 = m1 + m1 + m2 + m2  # test for more than two objects being added
-        expected_peaklist = sorted([(110, 0.5), (100, 1), (90, 1), (80, 1),
-                                    (70, 0.5)])
+        expected_peaklist = sorted([(110, 0.5, 0.5), (100, 1, 0.5), (90, 1, 0.5), (80, 1, 0.5),
+                                    (70, 0.5, 0.5)])
         assert np.allclose(s3.peaklist(), expected_peaklist)
         with pytest.raises(TypeError):
             _ = s + 1
 
     def test_add_appends_to_components(self):
-        m1 = Multiplet(100, 1, [(10, 2)])
-        m2 = Multiplet(80, 1, [(10, 2)])
+        m1 = Multiplet(100, 1, [(10, 2)], 0.5)
+        m2 = Multiplet(80, 1, [(10, 2)], 0.5)
         s = Spectrum([m1])
         s2 = s + m2
         assert s2._components == [m1, m2]
@@ -270,8 +270,8 @@ def test_lineshape(self):
         # also test vmin/vmax work
         spectrum = Spectrum([m1, m2], vmin=0.0, vmax=200.0)
         x = np.linspace(0.0, 200.0, 1000)
-        y1 = add_lorentzians(x, m1.peaklist(), m1.w)
-        y2 = add_lorentzians(x, m2.peaklist(), m2.w)
+        y1 = add_lorentzians_limitable(x, m1.peaklist())
+        y2 = add_lorentzians_limitable(x, m2.peaklist())
         y_sum = [sum(i) for i in zip(y1, y2)]
         spec_x, spec_y = spectrum.lineshape(points=1000)
         assert np.allclose(spec_x, x)
@@ -294,8 +294,8 @@ def test_add(self):
         spectrum2 = m3 + subspectrum
         assert spectrum._components == [m1, m2, m3]
         assert spectrum2._components == [m3, m1, m2]
-        expected_peaklist = sorted([(110, 0.25), (100, 0.5), (90, 0.5), (80, 0.5),
-                                    (70, 0.25), (50, 0.25), (40, 0.5), (30, 0.25)])
+        expected_peaklist = sorted([(110, 0.25, 0.5), (100, 0.5, 0.5), (90, 0.5, 0.5), (80, 0.5, 0.5),
+                                    (70, 0.25, 0.5), (50, 0.25, 0.5), (40, 0.5, 0.5), (30, 0.25, 0.5)])
         assert np.allclose(expected_peaklist, spectrum.peaklist())
         assert np.allclose(expected_peaklist, spectrum2.peaklist())
 
@@ -305,8 +305,8 @@ def test_iadd(self):
         m3 = Multiplet(40, 1, [(10, 2)])
         spectrum = Spectrum([m1, m2])
         spectrum += m3
-        expected_peaklist = sorted([(110, 0.25), (100, 0.5), (90, 0.5), (80, 0.5),
-                                    (70, 0.25), (50, 0.25), (40, 0.5), (30, 0.25)])
+        expected_peaklist = sorted([(110, 0.25, 0.5), (100, 0.5, 0.5), (90, 0.5, 0.5), (80, 0.5, 0.5),
+                                    (70, 0.25, 0.5), (50, 0.25, 0.5), (40, 0.5, 0.5), (30, 0.25, 0.5)])
         assert np.allclose(expected_peaklist, spectrum.peaklist())
         spectrum2 = Spectrum([m1])
         spectrum += spectrum2
diff --git a/tests/test_plt.py b/tests/test_plt.py
index 9bd92e0..b200c68 100644
--- a/tests/test_plt.py
+++ b/tests/test_plt.py
@@ -2,7 +2,7 @@
 import pytest
 
 from nmrsim.plt import mplplot, mplplot_stick, mplplot_lineshape
-from nmrsim.math import add_lorentzians
+from nmrsim.math import add_lorentzians_limitable
 from tests.accepted_data import ADD_SIGNALS_DATASET
 from tests.dnmr_standards import TWOSPIN_SLOW
 
@@ -18,8 +18,8 @@ def test_add_signals():
     # test was written before normalization of height vs width was built into
     #  lorentz(). Fudge-factor added to scale old accepted data.
     x = np.linspace(390, 410, 200)
-    doublet = [(399, 1), (401, 1)]
-    y = add_lorentzians(x, doublet, 1)
+    doublet = [(399, 1, 1), (401, 1, 1)]
+    y = add_lorentzians_limitable(x, doublet, None)
     X = np.array([x for x, _ in ADD_SIGNALS_DATASET])
     Y = np.array([y * 0.5 for _, y in ADD_SIGNALS_DATASET])
 
@@ -28,24 +28,28 @@ def test_add_signals():
 
 
 def test_mplplot_defaults():
-    doublet = [(399, 1), (401, 1)]
-    x, y = mplplot(doublet, hidden=True)
+    doublet = [(399, 1, 1), (401, 1, 1)]
+    x, y = mplplot(doublet, hidden=False)
     assert len(x) == 800
-    assert x[0] == 399 - 50
-    assert x[-1] == 401 + 50
+    assert x[0] == 401 + 50
+    assert x[-1] == 399 - 50
 
 
 @pytest.mark.parametrize('limits', ['foo', (1,), (1, 'foo'), (1, 2, 3)])
 def test_mplplot_limit_error(limits):
-    doublet = [(399, 1), (401, 1)]
+    doublet = [(399, 1, 1), (401, 1, 1)]
     with pytest.raises((AttributeError, TypeError, ValueError)):
         mplplot(doublet, limits=limits, hidden=True)  # noqa
 
 
 def test_mplplot():
-    doublet = [(399, 1), (401, 1)]
+    doublet = [(399, 1, 1), (401, 1, 1)]
     limits = (410, 390)  # deliberately opposite order
     x, y = mplplot(doublet, points=200, limits=limits, hidden=True)
+    # Added an inversion to fit the ADD_SIGNALS_DATASET because create_lineshape
+    # now inverts the line by default.
+    x = x[::-1]
+    y = y[::-1]
     # test was written before normalization of height vs width was built into
     #  lorentz(). Fudge-factor added to scale old accepted data.
     y = y * 2