From 3534864bb3e3e5973caf8ee836c40b19f80e46f2 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Thu, 12 Mar 2026 12:43:57 -0700
Subject: [PATCH 01/12] start at mismatch.py

---
 pvlib/mismatch.py | 97 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 97 insertions(+)
 create mode 100644 pvlib/mismatch.py

diff --git a/pvlib/mismatch.py b/pvlib/mismatch.py
new file mode 100644
index 0000000000..1083ceacd0
--- /dev/null
+++ b/pvlib/mismatch.py
@@ -0,0 +1,97 @@
+"""
+Contains functions for solving for DC power in arrays with mismatched conditions.
+
+"""
+import numpy as np
+import singlediode as _singlediode
+
+
+def _iv_series_lambertw(photocurrent, saturation_current, resistance_series,
+                        resistance_shunt, nNsVth, neg_v_limit=None,
+                        delta_i=0.001):
+    r'''Solve the IV curve for series-connected devices using a single diode
+    equivalent circuit model.
+
+    Uses a simplified model for reverse bias behavior, where current is
+    unbounded at a constant reverse bias voltage ``neg_v_limit``.
+
+    Input parameters photocurrent, saturation_current, resistance_series,
+    resistance_shunt, nNsVth may be arrays. If arrays, all must be
+    broadcastable to a common shape. The first dimension of each array
+    is time. The 2nd dimension is devices in series.
+
+    Parameters
+    ----------
+    photocurrent : numeric
+        photocurrent (A).
+    saturation_current : numeric
+        saturation current (A).
+    resistance_series : numeric
+        series resistance (ohm).
+    resistance_shunt : numeric
+        shunt resistance (ohm).
+    nNsVth : numeric
+        product of diode factor n, number of series cells Ns, and
+        thermal voltage (Vth), (V).
+    neg_v_limit : float, optional
+        Limit on reverse bias voltage, from cell breakdown voltage or reverse
+        bias diode activation voltage (V). Should be negative. For example,
+        if neg_v_limit=-5, then at V=-5 current is unbounded in the positive
+        direction.
+    delta_i : float, optional
+        Width of interval used to discretize current (A).
+
+    Returns
+    -------
+    None.
+
+    '''
+    # solve for Isc
+    isc = _singlediode._lambertw_i_from_v(
+        0., photocurrent, saturation_current, resistance_series,
+        resistance_shunt, nNsVth)
+
+    # discretize current from max(Isc) down to 0.
+    currents = np.arange(isc.max(), 0., step=-delta_i)
+
+    # shape all the arrays
+    # target shape is ntimes x ndevices x ncurrents
+    # use a dict so we can add axes using a loop
+    params = {'photocurrent': photocurrent,
+              'saturation_current': saturation_current,
+              'resistance_series': resistance_series,
+              'resistance_shunt': resistance_shunt,
+              'nNsVth': nNsVth}
+
+    for p in params:
+        if not isinstance(params[p], np.ndarray):
+            pass  # float, int
+        if len(params[p].shape) == 1:
+            params[p] = params[p][:, np.newaxis, np.newaxis]
+        elif len(params[p].shape) == 2:
+            params[p] = params[p][:, :, np.newaxis]
+        else:
+            pass  # already 3d
+
+    il, io, rs, rsh, a = (params[p] for p in params)
+
+    currents = currents[np.newaxis, np.newaxis, :]
+
+    il, io, rs, rsh, a, currents = np.broadcast_arrays(
+        il, io, rs, rsh, a, currents)
+
+    # solve voltages at each current for each IV curve
+    voltages = _singlediode._lambertw_v_from_i(
+        currents, il, io, rs, rsh, a)
+
+    # apply negative voltage limit
+    if neg_v_limit is not None:
+        voltages[voltages < neg_v_limit] = neg_v_limit
+
+    # add voltage at common current to get series voltage
+    voltage_sum = voltages.sum(axis=1)
+
+    # drop currents dimension for devices
+    currents = currents[:, 0, :]
+
+    return voltage_sum, currents

From 1c04fffe1cc4093a4c5a9efaf06b1d7e72ce6080 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Sun, 22 Mar 2026 16:51:56 -0700
Subject: [PATCH 02/12] work

---
 pvlib/mismatch.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pvlib/mismatch.py b/pvlib/mismatch.py
index 1083ceacd0..c8f33c7369 100644
--- a/pvlib/mismatch.py
+++ b/pvlib/mismatch.py
@@ -65,13 +65,13 @@ def _iv_series_lambertw(photocurrent, saturation_current, resistance_series,
 
     for p in params:
         if not isinstance(params[p], np.ndarray):
-            pass  # float, int
+            continue  # float, int
         if len(params[p].shape) == 1:
             params[p] = params[p][:, np.newaxis, np.newaxis]
         elif len(params[p].shape) == 2:
             params[p] = params[p][:, :, np.newaxis]
         else:
-            pass  # already 3d
+            continue  # already 3d
 
     il, io, rs, rsh, a = (params[p] for p in params)
 

From 824c74660eb6d155d8d02bb1c55b842c9c2bde22 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Mon, 17 Aug 2026 14:44:26 -0600
Subject: [PATCH 03/12] update series calculation

---
 pvlib/mismatch.py | 230 ++++++++++++++++++++++++++++++++++------------
 1 file changed, 171 insertions(+), 59 deletions(-)

diff --git a/pvlib/mismatch.py b/pvlib/mismatch.py
index c8f33c7369..4d20141e5d 100644
--- a/pvlib/mismatch.py
+++ b/pvlib/mismatch.py
@@ -1,97 +1,209 @@
 """
-Contains functions for solving for DC power in arrays with mismatched conditions.
+Contains functions for solving for DC power in arrays with mismatched
+conditions.
 
 """
 import numpy as np
+from scipy.optimize.elementwise import find_root
 import singlediode as _singlediode
 
 
+def _iv_series_lambert_v_from_i(I, il, io, rs, rsh, a, neg_v_limit,
+                                ndevices=None, idx=None):
+    # solve voltages at each current for each IV curve
+    if ndevices is not None:
+        # broadcast I to ndevices to apply same current each device
+        I = np.broadcast_to(I[np.newaxis, :], (ndevices, len(I)))
+    # slice each parameter on its ntimes dimension with idx
+    if idx is not None:
+        il, io, rs, rsh, a = (il[:, idx], io[:, idx], rs[:, idx], rsh[:, idx],
+                              a[:, idx])
+
+    voltages = _singlediode._lambertw_v_from_i(
+        I.flatten(), il.flatten(), io.flatten(), rs.flatten(), rsh.flatten(),
+        a.flatten())
+
+    # apply negative voltage limit
+    voltages[voltages < neg_v_limit] = neg_v_limit
+
+    # reshape
+    voltages = voltages.reshape(I.shape)
+    return voltages
+
+
+def _setup_currents(device_isc, string_isc, npts):
+    r''' Form array of currents. Array of currents will contain all values
+    from device_isc which are less than string_isc.
+
+    Parameters
+    ----------
+    device_isc : ndarray
+        Shape (ndevices, ntimes)
+    string_isc : ndarray
+        Shape (ntimes,)
+    npts : int
+        number of current points in the returned array
+
+    Returns
+    -------
+    ndarray
+        shape (ntimes, npts)
+
+    '''
+    ntimes = len(string_isc)
+    A = np.zeros((ntimes, npts))
+
+    u = device_isc < string_isc[np.newaxis, :]
+
+    # have to loop on ntimes since count of device_isc < string_isc may
+    # differ for each time
+    for i in range(ntimes):
+        vals = np.unique(device_isc[u[:, i], i])
+        # ensure string_isc and 0. are added
+        vals = np.append(vals, [string_isc[i], 0.])
+        k_i = len(vals)
+
+        if k_i == 0:
+            continue
+
+        # Copy original values
+        A[i, :k_i] = vals
+
+        n_fill = npts - k_i
+        if n_fill <= 0:
+            continue
+
+        # Build grid
+        grid = np.linspace(string_isc[i], 0., npts)
+
+        # Compute distance to nearest point in arr
+        # shape: (grid_size, k_i)
+        dists = np.abs(grid[:, None] - vals[None, :])
+
+        # nearest distance per grid point
+        nearest_dist = np.min(dists, axis=1)
+
+        # inverse-distance weights (higher near arr values)
+        weights = 1.0 / (nearest_dist + 1e-12)  # 1e-12 to avoid div by 0
+
+        # Avoid re-selecting original values exactly
+        mask_existing = np.isclose(nearest_dist, 0.0, atol=1e-12)
+        weights[mask_existing] = 0.0
+
+        # Select top-weighted grid points
+        idx = np.argpartition(weights, -n_fill)[-n_fill:]
+        selected = grid[idx]
+
+        # Combine and add to A
+        A[i, :] = np.concatenate([vals, selected])
+
+    # return sorted in descending order for each time
+    A = -np.sort(-A, axis=1)
+
+    return A
+
+
 def _iv_series_lambertw(photocurrent, saturation_current, resistance_series,
-                        resistance_shunt, nNsVth, neg_v_limit=None,
-                        delta_i=0.001):
-    r'''Solve the IV curve for series-connected devices using a single diode
-    equivalent circuit model.
+                        resistance_shunt, nNsVth, neg_v_limit=0.,
+                        npts=100):
+    r'''Solve the IV curve for series-connected devices where each device
+    is described by the single diode equation.
 
     Uses a simplified model for reverse bias behavior, where current is
     unbounded at a constant reverse bias voltage ``neg_v_limit``.
 
-    Input parameters photocurrent, saturation_current, resistance_series,
-    resistance_shunt, nNsVth may be arrays. If arrays, all must be
-    broadcastable to a common shape. The first dimension of each array
-    is time. The 2nd dimension is devices in series.
+    Input parameter ``photocurrent`` must have shape (devices, times).
+    Input parameters ``saturation_current``, ``resistance_series``,
+    ``resistance_shunt``, ``nNsVth`` may be arrays. If arrays, must be
+    broadcastable to the shape of ``photocurrent``.
 
     Parameters
     ----------
     photocurrent : numeric
-        photocurrent (A).
+        photocurrent (A). Must have shape (devices, times).
     saturation_current : numeric
-        saturation current (A).
+        saturation current (A). Must be broadcastable with photocurrent.
     resistance_series : numeric
-        series resistance (ohm).
+        series resistance (ohm). Must be broadcastable with photocurrent.
     resistance_shunt : numeric
-        shunt resistance (ohm).
+        shunt resistance (ohm). Must be broadcastable with photocurrent.
     nNsVth : numeric
         product of diode factor n, number of series cells Ns, and
-        thermal voltage (Vth), (V).
+        thermal voltage (Vth), (V). Must be broadcastable with photocurrent.
     neg_v_limit : float, optional
         Limit on reverse bias voltage, from cell breakdown voltage or reverse
         bias diode activation voltage (V). Should be negative. For example,
         if neg_v_limit=-5, then at V=-5 current is unbounded in the positive
         direction.
-    delta_i : float, optional
-        Width of interval used to discretize current (A).
+    npts : int, optional
+        Number of points used to discretize the returned IV curves.
 
     Returns
     -------
-    None.
+    voltages : numeric
+        Voltage points for the series IV curves (V), shape
+        (times, npts).
+    currents : numeric
+        Current points for the series IV curves (A), shape
+        (times, npts).
 
     '''
+    # target shape is ndevices x ntimes
+    IL, I0, Rs, Rsh, a = \
+        np.broadcast_arrays(photocurrent, saturation_current,
+                            resistance_series, resistance_shunt, nNsVth)
+
+    ndevices, ntimes = IL.shape
+
     # solve for Isc
-    isc = _singlediode._lambertw_i_from_v(
-        0., photocurrent, saturation_current, resistance_series,
-        resistance_shunt, nNsVth)
-
-    # discretize current from max(Isc) down to 0.
-    currents = np.arange(isc.max(), 0., step=-delta_i)
-
-    # shape all the arrays
-    # target shape is ntimes x ndevices x ncurrents
-    # use a dict so we can add axes using a loop
-    params = {'photocurrent': photocurrent,
-              'saturation_current': saturation_current,
-              'resistance_series': resistance_series,
-              'resistance_shunt': resistance_shunt,
-              'nNsVth': nNsVth}
-
-    for p in params:
-        if not isinstance(params[p], np.ndarray):
-            continue  # float, int
-        if len(params[p].shape) == 1:
-            params[p] = params[p][:, np.newaxis, np.newaxis]
-        elif len(params[p].shape) == 2:
-            params[p] = params[p][:, :, np.newaxis]
-        else:
-            continue  # already 3d
-
-    il, io, rs, rsh, a = (params[p] for p in params)
-
-    currents = currents[np.newaxis, np.newaxis, :]
-
-    il, io, rs, rsh, a, currents = np.broadcast_arrays(
-        il, io, rs, rsh, a, currents)
+    # Isc for each device
+    device_isc = _singlediode._lambertw_i_from_v(
+        neg_v_limit, IL, I0, Rs, Rsh, a)
+
+    # find Isc for series of devices
+    # 1d array contains bounds for each time
+    # add/substract eps in case max==min
+    max_isc = device_isc.max(axis=0) * 1.01
+    min_isc = device_isc.min(axis=0) * 0.99
+
+    # use index idx so that find_root can slice arguments
+    # internally find_root will slice current as each element converges
+    # as an argument, idx lets find_root also slice the other parameters
+    # remove idx and use preserve_shape once available in find_root
+    # https://github.com/scipy/scipy/issues/24869
+    idx = np.arange(ntimes)
+    def optfn(current, idx):
+        # current is ntimes only since it is common for all devices.
+        # other parameters are ntimes x ndevices
+        v = _iv_series_lambert_v_from_i(
+            current, IL, I0, Rs,
+            Rsh, a, neg_v_limit, ndevices, idx)
+        # return string voltage
+        return v.sum(axis=0)
+
+    isc_result = find_root(
+        optfn,
+        (min_isc, max_isc), args=(idx,))
+    string_isc = isc_result.x  # 1d in ntimes
+
+    # discretize current from max(Isc) down to 0. at each time step
+    # Include each device Isc (except the highest) so that the series IV curve
+    current_pts = _setup_currents(device_isc, string_isc, npts)
+
+    # shape all arrays to be ndevices x ntimes x ncurrents
+    I = np.repeat(current_pts[np.newaxis, :, :], ndevices, axis=0)
+    I, IL, I0, Rs, Rsh, a = np.broadcast_arrays(
+        I, IL[:, :, np.newaxis], I0[:, :, np.newaxis], Rs[:, :, np.newaxis],
+        Rsh[:, :, np.newaxis], a[:, :, np.newaxis])
 
     # solve voltages at each current for each IV curve
-    voltages = _singlediode._lambertw_v_from_i(
-        currents, il, io, rs, rsh, a)
-
-    # apply negative voltage limit
-    if neg_v_limit is not None:
-        voltages[voltages < neg_v_limit] = neg_v_limit
+    voltages = _iv_series_lambert_v_from_i(
+        I, IL, I0, Rs, Rsh, a, neg_v_limit)
 
-    # add voltage at common current to get series voltage
-    voltage_sum = voltages.sum(axis=1)
+    # add voltage across devices to get series voltage
+    voltage_sum = voltages.sum(axis=0)
 
     # drop currents dimension for devices
-    currents = currents[:, 0, :]
+    I = I[0, :, :]
 
-    return voltage_sum, currents
+    return voltage_sum, I

From 80b952050017e2302926b54dfd37bd2abe6644ed Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Mon, 17 Aug 2026 16:05:00 -0600
Subject: [PATCH 04/12] rename some 1 letter variables

---
 pvlib/mismatch.py | 21 +++++++++++----------
 1 file changed, 11 insertions(+), 10 deletions(-)

diff --git a/pvlib/mismatch.py b/pvlib/mismatch.py
index 4d20141e5d..d70751b071 100644
--- a/pvlib/mismatch.py
+++ b/pvlib/mismatch.py
@@ -8,26 +8,27 @@
 import singlediode as _singlediode
 
 
-def _iv_series_lambert_v_from_i(I, il, io, rs, rsh, a, neg_v_limit,
+def _iv_series_lambert_v_from_i(current, il, io, rs, rsh, a, neg_v_limit,
                                 ndevices=None, idx=None):
     # solve voltages at each current for each IV curve
     if ndevices is not None:
         # broadcast I to ndevices to apply same current each device
-        I = np.broadcast_to(I[np.newaxis, :], (ndevices, len(I)))
+        current = np.broadcast_to(current[np.newaxis, :],
+                                  (ndevices, len(current)))
     # slice each parameter on its ntimes dimension with idx
     if idx is not None:
         il, io, rs, rsh, a = (il[:, idx], io[:, idx], rs[:, idx], rsh[:, idx],
                               a[:, idx])
 
     voltages = _singlediode._lambertw_v_from_i(
-        I.flatten(), il.flatten(), io.flatten(), rs.flatten(), rsh.flatten(),
-        a.flatten())
+        current.flatten(), il.flatten(), io.flatten(), rs.flatten(),
+        rsh.flatten(), a.flatten())
 
     # apply negative voltage limit
     voltages[voltages < neg_v_limit] = neg_v_limit
 
     # reshape
-    voltages = voltages.reshape(I.shape)
+    voltages = voltages.reshape(current.shape)
     return voltages
 
 
@@ -51,7 +52,7 @@ def _setup_currents(device_isc, string_isc, npts):
 
     '''
     ntimes = len(string_isc)
-    A = np.zeros((ntimes, npts))
+    currents = np.zeros((ntimes, npts))
 
     u = device_isc < string_isc[np.newaxis, :]
 
@@ -67,7 +68,7 @@ def _setup_currents(device_isc, string_isc, npts):
             continue
 
         # Copy original values
-        A[i, :k_i] = vals
+        currents[i, :k_i] = vals
 
         n_fill = npts - k_i
         if n_fill <= 0:
@@ -95,12 +96,12 @@ def _setup_currents(device_isc, string_isc, npts):
         selected = grid[idx]
 
         # Combine and add to A
-        A[i, :] = np.concatenate([vals, selected])
+        currents[i, :] = np.concatenate([vals, selected])
 
     # return sorted in descending order for each time
-    A = -np.sort(-A, axis=1)
+    currents = -np.sort(-currents, axis=1)
 
-    return A
+    return currents
 
 
 def _iv_series_lambertw(photocurrent, saturation_current, resistance_series,

From 0918ad16ce27f51e26f32d2951d45c96481b6baa Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Tue, 18 Aug 2026 13:46:03 -0600
Subject: [PATCH 05/12] tests for series combinations

---
 tests/test_mismatch.py | 77 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 77 insertions(+)
 create mode 100644 tests/test_mismatch.py

diff --git a/tests/test_mismatch.py b/tests/test_mismatch.py
new file mode 100644
index 0000000000..8914678c35
--- /dev/null
+++ b/tests/test_mismatch.py
@@ -0,0 +1,77 @@
+import numpy as np
+from pvlib import singlediode as _singlediode
+from pvlib.mismatch import _setup_currents, _iv_series_lambertw
+
+
+def test__setup_currents():
+
+    # 2 devices, 1 time
+    cur_bkpts = np.array([[4.1, 5.2]]).T
+    string_isc = np.array([6.])
+    curs = _setup_currents(cur_bkpts, string_isc, 10)
+    assert np.isin(cur_bkpts, curs).all()
+    assert (np.diff(curs) < 0.).all()  # strictly decreasing
+    assert curs[:, -1] == 0.
+    assert curs[:, 0] == string_isc
+
+
+def test__iv_series_lambertw_isc():
+
+    # 2 devices, 1 time
+    # Isc should be equal to the current on the higher curve where voltage is
+    # +neg_v_limit
+    IL = np.array([[1.0], [6.01]])
+    Io = 1e-9
+    nNsVth = 2.5
+    Rs = 0.5
+    Rsh = 1000.
+    npts = 5
+    neg_v_limit = -5.
+    expected_isc = _singlediode._lambertw_i_from_v(-neg_v_limit, IL[1], Io, Rs,
+                                                   Rsh, nNsVth)
+    vs, cs = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth, neg_v_limit=neg_v_limit,
+                                 npts=npts)
+    assert vs.shape == cs.shape
+    assert vs.shape == (1, npts)  # ntimes x npts
+    assert np.isclose(cs[0, 0], expected_isc)
+    assert np.isclose(vs[0, 0], 0.)
+
+
+def test__iv_series_lambertw_voc():
+
+    # 2 devices, 1 time
+    # Voc should be equal to the sum of Voc for each device
+    IL = np.array([[1.0], [6.01]])
+    Io = 1e-9
+    nNsVth = 2.5
+    Rs = 0.5
+    Rsh = 1000.
+    npts = 5
+    neg_v_limit = -5.
+    expected_voc = _singlediode._lambertw_v_from_i(0, IL, Io, Rs,
+                                                   Rsh, nNsVth)
+    vs, cs = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth, neg_v_limit=neg_v_limit,
+                                 npts=npts)
+    
+    assert np.isclose(vs[0, -1], expected_voc.sum())
+    assert np.isclose(cs[0, -1], 0.)
+
+
+def test__iv_series_lambertw_breakpoints():
+    
+    # 2 devices, 1 time
+    # Voc should be equal to the sum of Voc for each device
+    IL = np.array([[1.0], [6.01]])
+    Io = 1e-9
+    nNsVth = 2.5
+    Rs = 0.5
+    Rsh = 1000.
+    npts = 10
+    neg_v_limit = -5.
+    
+    cur_breakpoints = _singlediode._lambertw_i_from_v(neg_v_limit, IL, Io, Rs,
+                                                      Rsh, nNsVth)
+    vs, cs = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth, neg_v_limit=neg_v_limit,
+                                 npts=npts)
+    # current cs should contain all breakpoints less than string_isc
+    assert np.isin(cur_breakpoints[cur_breakpoints < cs[0, 0]], cs).all()

From d51a6a3483373a640cd87704b0b039eb39e951d9 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Tue, 18 Aug 2026 13:54:05 -0600
Subject: [PATCH 06/12] fix import

---
 pvlib/mismatch.py | 43 +++++++++++++++++++++++--------------------
 1 file changed, 23 insertions(+), 20 deletions(-)

diff --git a/pvlib/mismatch.py b/pvlib/mismatch.py
index d70751b071..a9d2cb06e9 100644
--- a/pvlib/mismatch.py
+++ b/pvlib/mismatch.py
@@ -5,11 +5,13 @@
 """
 import numpy as np
 from scipy.optimize.elementwise import find_root
-import singlediode as _singlediode
+from pvlib import singlediode as _singlediode
 
 
 def _iv_series_lambert_v_from_i(current, il, io, rs, rsh, a, neg_v_limit,
                                 ndevices=None, idx=None):
+    # wrapper for pvlib._singlediode._lambertw_v_from_i, handles
+    # dimensions expected in series calculation
     # solve voltages at each current for each IV curve
     if ndevices is not None:
         # broadcast I to ndevices to apply same current each device
@@ -32,9 +34,10 @@ def _iv_series_lambert_v_from_i(current, il, io, rs, rsh, a, neg_v_limit,
     return voltages
 
 
-def _setup_currents(device_isc, string_isc, npts):
-    r''' Form array of currents. Array of currents will contain all values
-    from device_isc which are less than string_isc.
+def _setup_currents(current_bkpts, string_isc, npts):
+    r''' Form array of currents from string_isc down to 0.
+    The array of currents will contain all values
+    from current_bkpts which are less than string_isc.
 
     Parameters
     ----------
@@ -54,12 +57,12 @@ def _setup_currents(device_isc, string_isc, npts):
     ntimes = len(string_isc)
     currents = np.zeros((ntimes, npts))
 
-    u = device_isc < string_isc[np.newaxis, :]
+    u = current_bkpts < string_isc[np.newaxis, :]
 
     # have to loop on ntimes since count of device_isc < string_isc may
     # differ for each time
     for i in range(ntimes):
-        vals = np.unique(device_isc[u[:, i], i])
+        vals = np.unique(current_bkpts[u[:, i], i])
         # ensure string_isc and 0. are added
         vals = np.append(vals, [string_isc[i], 0.])
         k_i = len(vals)
@@ -156,21 +159,20 @@ def _iv_series_lambertw(photocurrent, saturation_current, resistance_series,
 
     ndevices, ntimes = IL.shape
 
-    # solve for Isc
-    # Isc for each device
-    device_isc = _singlediode._lambertw_i_from_v(
+    # solve for current at negative voltage limit for each device.
+    # these currents create breakpoints in the series IV curve
+    current_bkpts = _singlediode._lambertw_i_from_v(
         neg_v_limit, IL, I0, Rs, Rsh, a)
 
-    # find Isc for series of devices
-    # 1d array contains bounds for each time
-    # add/substract eps in case max==min
-    max_isc = device_isc.max(axis=0) * 1.01
-    min_isc = device_isc.min(axis=0) * 0.99
+    # find Isc for string IV curve
+    # bounds, 1d array for each time
+    max_isc = current_bkpts.max(axis=0) * 1.01
+    min_isc = current_bkpts.min(axis=0) * 0.99
 
-    # use index idx so that find_root can slice arguments
-    # internally find_root will slice current as each element converges
-    # as an argument, idx lets find_root also slice the other parameters
-    # remove idx and use preserve_shape once available in find_root
+    # Use an index idx so that find_root can slice arguments
+    # Internally find_root will slice current as each element converges
+    # As an argument, idx lets find_root also slice the other parameters
+    # Remove idx and use preserve_shape once available in find_root
     # https://github.com/scipy/scipy/issues/24869
     idx = np.arange(ntimes)
     def optfn(current, idx):
@@ -187,9 +189,10 @@ def optfn(current, idx):
         (min_isc, max_isc), args=(idx,))
     string_isc = isc_result.x  # 1d in ntimes
 
-    # discretize current from max(Isc) down to 0. at each time step
+    # discretize current from string_isc down to 0 at each time step
     # Include each device Isc (except the highest) so that the series IV curve
-    current_pts = _setup_currents(device_isc, string_isc, npts)
+    # includes the breakpoints
+    current_pts = _setup_currents(current_bkpts, string_isc, npts)
 
     # shape all arrays to be ndevices x ntimes x ncurrents
     I = np.repeat(current_pts[np.newaxis, :, :], ndevices, axis=0)

From 8998ef591f0fe1f292799d6b12d84eb07456c643 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Tue, 18 Aug 2026 13:59:18 -0600
Subject: [PATCH 07/12] flake8

---
 pvlib/mismatch.py      | 13 +++++++------
 tests/test_mismatch.py | 15 +++++++++------
 2 files changed, 16 insertions(+), 12 deletions(-)

diff --git a/pvlib/mismatch.py b/pvlib/mismatch.py
index a9d2cb06e9..2af7f0148a 100644
--- a/pvlib/mismatch.py
+++ b/pvlib/mismatch.py
@@ -175,6 +175,7 @@ def _iv_series_lambertw(photocurrent, saturation_current, resistance_series,
     # Remove idx and use preserve_shape once available in find_root
     # https://github.com/scipy/scipy/issues/24869
     idx = np.arange(ntimes)
+
     def optfn(current, idx):
         # current is ntimes only since it is common for all devices.
         # other parameters are ntimes x ndevices
@@ -195,19 +196,19 @@ def optfn(current, idx):
     current_pts = _setup_currents(current_bkpts, string_isc, npts)
 
     # shape all arrays to be ndevices x ntimes x ncurrents
-    I = np.repeat(current_pts[np.newaxis, :, :], ndevices, axis=0)
-    I, IL, I0, Rs, Rsh, a = np.broadcast_arrays(
-        I, IL[:, :, np.newaxis], I0[:, :, np.newaxis], Rs[:, :, np.newaxis],
+    curs = np.repeat(current_pts[np.newaxis, :, :], ndevices, axis=0)
+    curs, IL, I0, Rs, Rsh, a = np.broadcast_arrays(
+        curs, IL[:, :, np.newaxis], I0[:, :, np.newaxis], Rs[:, :, np.newaxis],
         Rsh[:, :, np.newaxis], a[:, :, np.newaxis])
 
     # solve voltages at each current for each IV curve
     voltages = _iv_series_lambert_v_from_i(
-        I, IL, I0, Rs, Rsh, a, neg_v_limit)
+        curs, IL, I0, Rs, Rsh, a, neg_v_limit)
 
     # add voltage across devices to get series voltage
     voltage_sum = voltages.sum(axis=0)
 
     # drop currents dimension for devices
-    I = I[0, :, :]
+    curs = curs[0, :, :]
 
-    return voltage_sum, I
+    return voltage_sum, curs
diff --git a/tests/test_mismatch.py b/tests/test_mismatch.py
index 8914678c35..4c097df070 100644
--- a/tests/test_mismatch.py
+++ b/tests/test_mismatch.py
@@ -29,7 +29,8 @@ def test__iv_series_lambertw_isc():
     neg_v_limit = -5.
     expected_isc = _singlediode._lambertw_i_from_v(-neg_v_limit, IL[1], Io, Rs,
                                                    Rsh, nNsVth)
-    vs, cs = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth, neg_v_limit=neg_v_limit,
+    vs, cs = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth,
+                                 neg_v_limit=neg_v_limit,
                                  npts=npts)
     assert vs.shape == cs.shape
     assert vs.shape == (1, npts)  # ntimes x npts
@@ -50,15 +51,16 @@ def test__iv_series_lambertw_voc():
     neg_v_limit = -5.
     expected_voc = _singlediode._lambertw_v_from_i(0, IL, Io, Rs,
                                                    Rsh, nNsVth)
-    vs, cs = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth, neg_v_limit=neg_v_limit,
+    vs, cs = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth,
+                                 neg_v_limit=neg_v_limit,
                                  npts=npts)
-    
+
     assert np.isclose(vs[0, -1], expected_voc.sum())
     assert np.isclose(cs[0, -1], 0.)
 
 
 def test__iv_series_lambertw_breakpoints():
-    
+
     # 2 devices, 1 time
     # Voc should be equal to the sum of Voc for each device
     IL = np.array([[1.0], [6.01]])
@@ -68,10 +70,11 @@ def test__iv_series_lambertw_breakpoints():
     Rsh = 1000.
     npts = 10
     neg_v_limit = -5.
-    
+
     cur_breakpoints = _singlediode._lambertw_i_from_v(neg_v_limit, IL, Io, Rs,
                                                       Rsh, nNsVth)
-    vs, cs = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth, neg_v_limit=neg_v_limit,
+    vs, cs = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth,
+                                 neg_v_limit=neg_v_limit,
                                  npts=npts)
     # current cs should contain all breakpoints less than string_isc
     assert np.isin(cur_breakpoints[cur_breakpoints < cs[0, 0]], cs).all()

From b0c9e71070d473c2ef4dac45d166650dfe65c06b Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Wed, 19 Aug 2026 06:32:32 -0600
Subject: [PATCH 08/12] add to docs

---
 .../effects_on_pv_system_output/electrical-mismatch.rst  | 9 +++++++++
 .../reference/effects_on_pv_system_output/index.rst      | 1 +
 docs/sphinx/source/whatsnew/v0.15.3.rst                  | 3 +++
 3 files changed, 13 insertions(+)
 create mode 100644 docs/sphinx/source/reference/effects_on_pv_system_output/electrical-mismatch.rst

diff --git a/docs/sphinx/source/reference/effects_on_pv_system_output/electrical-mismatch.rst b/docs/sphinx/source/reference/effects_on_pv_system_output/electrical-mismatch.rst
new file mode 100644
index 0000000000..7723287917
--- /dev/null
+++ b/docs/sphinx/source/reference/effects_on_pv_system_output/electrical-mismatch.rst
@@ -0,0 +1,9 @@
+.. currentmodule:: pvlib
+
+Electrical mismatch
+-------------------
+
+.. autosummary::
+   :toctree: ../generated/
+
+   mismatch._iv_series_lambertw
diff --git a/docs/sphinx/source/reference/effects_on_pv_system_output/index.rst b/docs/sphinx/source/reference/effects_on_pv_system_output/index.rst
index cb046de7db..a2e0dc2de3 100644
--- a/docs/sphinx/source/reference/effects_on_pv_system_output/index.rst
+++ b/docs/sphinx/source/reference/effects_on_pv_system_output/index.rst
@@ -11,3 +11,4 @@ Effects on PV System Output
    soiling
    shading
    spectrum
+   electrical-mismatch
diff --git a/docs/sphinx/source/whatsnew/v0.15.3.rst b/docs/sphinx/source/whatsnew/v0.15.3.rst
index 8cf8515577..a5b4339406 100644
--- a/docs/sphinx/source/whatsnew/v0.15.3.rst
+++ b/docs/sphinx/source/whatsnew/v0.15.3.rst
@@ -68,6 +68,8 @@ Enhancements
 * Accelerate :py:func:`~pvlib.bifacial.utils.vf_ground_sky_2d_integ` by one or
   two orders of magnitude. This also makes :py:mod:`pvlib.bifacial.infinite_sheds` faster.
   (:pull:`2740`)
+* Add `:py:mod:`~pvlib.mismatch` to contain functions for computing DC power
+  and IV curves accounting for electrical mismatch (:pull:`2718`)
 
 
 Documentation
@@ -120,3 +122,4 @@ Contributors
 * Leonardo Scappatura (:ghuser:`Leonard013`)
 * Carolina Crespo (:ghuser:`cbcrespo`)
 * Ioannis Sifnaios (:ghuser:`IoannisSifnaios`)
+* Cliff Hansen (:ghuser:`cwhanse`)

From 9fcf6a9d2a02a081af74db833da15b95ac5a8d15 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Wed, 19 Aug 2026 06:36:02 -0600
Subject: [PATCH 09/12] oops to v0.16.0

---
 docs/sphinx/source/whatsnew/v0.16.0.rst | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/docs/sphinx/source/whatsnew/v0.16.0.rst b/docs/sphinx/source/whatsnew/v0.16.0.rst
index 3ea266e6cb..8b7ea2b9ec 100644
--- a/docs/sphinx/source/whatsnew/v0.16.0.rst
+++ b/docs/sphinx/source/whatsnew/v0.16.0.rst
@@ -66,7 +66,8 @@ Enhancements
   (:issue:`2828`, :pull:`2832`)
 * Allow variables from multiple datasets to be requested at once in
   :py:func:`~pvlib.iotools.get_merra2`. (:pull:`2839`)
-
+* Add `:py:mod:`~pvlib.mismatch` to contain functions for computing DC power
+  and IV curves accounting for electrical mismatch (:pull:`2718`)
 
 Documentation
 ~~~~~~~~~~~~~
@@ -97,3 +98,4 @@ Contributors
 * Sai Asish Y (:ghuser:`SAY-5`)
 * Kevin Anderson (:ghuser:`kandersolar`)
 * Johann Loux (:ghuser:`JoLo90`)
+* Cliff Hansen (:ghuser:`cwhanse`)

From 64dcfbf93c99865828ed6098a6121274a765fe58 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Wed, 19 Aug 2026 06:38:13 -0600
Subject: [PATCH 10/12] not to v0.15.3

---
 docs/sphinx/source/whatsnew/v0.15.3.rst | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/docs/sphinx/source/whatsnew/v0.15.3.rst b/docs/sphinx/source/whatsnew/v0.15.3.rst
index a5b4339406..8cf8515577 100644
--- a/docs/sphinx/source/whatsnew/v0.15.3.rst
+++ b/docs/sphinx/source/whatsnew/v0.15.3.rst
@@ -68,8 +68,6 @@ Enhancements
 * Accelerate :py:func:`~pvlib.bifacial.utils.vf_ground_sky_2d_integ` by one or
   two orders of magnitude. This also makes :py:mod:`pvlib.bifacial.infinite_sheds` faster.
   (:pull:`2740`)
-* Add `:py:mod:`~pvlib.mismatch` to contain functions for computing DC power
-  and IV curves accounting for electrical mismatch (:pull:`2718`)
 
 
 Documentation
@@ -122,4 +120,3 @@ Contributors
 * Leonardo Scappatura (:ghuser:`Leonard013`)
 * Carolina Crespo (:ghuser:`cbcrespo`)
 * Ioannis Sifnaios (:ghuser:`IoannisSifnaios`)
-* Cliff Hansen (:ghuser:`cwhanse`)

From 62925ce05f675549a43290ca02eb481b89d00253 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Mon, 24 Aug 2026 14:59:45 -0600
Subject: [PATCH 11/12] add mpp calculation

---
 pvlib/mismatch.py      | 87 +++++++++++++++++++++++++++++-------------
 tests/test_mismatch.py | 18 ++++-----
 2 files changed, 69 insertions(+), 36 deletions(-)

diff --git a/pvlib/mismatch.py b/pvlib/mismatch.py
index 2af7f0148a..e41294f51a 100644
--- a/pvlib/mismatch.py
+++ b/pvlib/mismatch.py
@@ -4,7 +4,7 @@
 
 """
 import numpy as np
-from scipy.optimize.elementwise import find_root
+from scipy.optimize.elementwise import find_root, find_minimum
 from pvlib import singlediode as _singlediode
 
 
@@ -13,10 +13,7 @@ def _iv_series_lambert_v_from_i(current, il, io, rs, rsh, a, neg_v_limit,
     # wrapper for pvlib._singlediode._lambertw_v_from_i, handles
     # dimensions expected in series calculation
     # solve voltages at each current for each IV curve
-    if ndevices is not None:
-        # broadcast I to ndevices to apply same current each device
-        current = np.broadcast_to(current[np.newaxis, :],
-                                  (ndevices, len(current)))
+
     # slice each parameter on its ntimes dimension with idx
     if idx is not None:
         il, io, rs, rsh, a = (il[:, idx], io[:, idx], rs[:, idx], rsh[:, idx],
@@ -37,7 +34,9 @@ def _iv_series_lambert_v_from_i(current, il, io, rs, rsh, a, neg_v_limit,
 def _setup_currents(current_bkpts, string_isc, npts):
     r''' Form array of currents from string_isc down to 0.
     The array of currents will contain all values
-    from current_bkpts which are less than string_isc.
+    from current_bkpts which are less than string_isc. Remaining points are
+    selected from a linear spacing [0, string_isc], avoiding points that are
+    closest to values in current_bkpts.
 
     Parameters
     ----------
@@ -62,32 +61,34 @@ def _setup_currents(current_bkpts, string_isc, npts):
     # have to loop on ntimes since count of device_isc < string_isc may
     # differ for each time
     for i in range(ntimes):
+        # start with values in current_bkpts
         vals = np.unique(current_bkpts[u[:, i], i])
-        # ensure string_isc and 0. are added
+        # add string_isc and 0.
         vals = np.append(vals, [string_isc[i], 0.])
         k_i = len(vals)
 
         if k_i == 0:
             continue
 
-        # Copy original values
+        # Store these values
         currents[i, :k_i] = vals
 
+        # how many more we need
         n_fill = npts - k_i
         if n_fill <= 0:
             continue
 
-        # Build grid
+        # Build linear spacing to draw from
         grid = np.linspace(string_isc[i], 0., npts)
 
-        # Compute distance to nearest point in arr
+        # Compute distances to nearest points in current_bkpts etc.
         # shape: (grid_size, k_i)
         dists = np.abs(grid[:, None] - vals[None, :])
 
         # nearest distance per grid point
         nearest_dist = np.min(dists, axis=1)
 
-        # inverse-distance weights (higher near arr values)
+        # inverse-distance weights (higher near must-have values)
         weights = 1.0 / (nearest_dist + 1e-12)  # 1e-12 to avoid div by 0
 
         # Avoid re-selecting original values exactly
@@ -175,40 +176,72 @@ def _iv_series_lambertw(photocurrent, saturation_current, resistance_series,
     # Remove idx and use preserve_shape once available in find_root
     # https://github.com/scipy/scipy/issues/24869
     idx = np.arange(ntimes)
-
-    def optfn(current, idx):
+    def isc_optfn(current, idx):
         # current is ntimes only since it is common for all devices.
         # other parameters are ntimes x ndevices
+        # broadcast current to ntimes x ndevices
+        cur2 = np.broadcast_to(current[np.newaxis, :], (ndevices, len(current)))
         v = _iv_series_lambert_v_from_i(
-            current, IL, I0, Rs,
+            cur2, IL, I0, Rs,
             Rsh, a, neg_v_limit, ndevices, idx)
         # return string voltage
         return v.sum(axis=0)
 
     isc_result = find_root(
-        optfn,
+        isc_optfn,
         (min_isc, max_isc), args=(idx,))
     string_isc = isc_result.x  # 1d in ntimes
 
     # discretize current from string_isc down to 0 at each time step
+    # Include each device's current at neg_v_limit since these will be the
+    # curvature breakpoints in the series IV curve.
+    # discretize current from max(Isc) down to 0. at each time step
     # Include each device Isc (except the highest) so that the series IV curve
-    # includes the breakpoints
-    current_pts = _setup_currents(current_bkpts, string_isc, npts)
+    currents = _setup_currents(current_bkpts, string_isc, npts)
 
     # shape all arrays to be ndevices x ntimes x ncurrents
-    curs = np.repeat(current_pts[np.newaxis, :, :], ndevices, axis=0)
-    curs, IL, I0, Rs, Rsh, a = np.broadcast_arrays(
-        curs, IL[:, :, np.newaxis], I0[:, :, np.newaxis], Rs[:, :, np.newaxis],
+    cur3 = np.repeat(currents[np.newaxis, :, :], ndevices, axis=0)
+    cur3, il, io, rs, rsh, a3 = np.broadcast_arrays(
+        cur3, IL[:, :, np.newaxis], I0[:, :, np.newaxis], Rs[:, :, np.newaxis],
         Rsh[:, :, np.newaxis], a[:, :, np.newaxis])
 
     # solve voltages at each current for each IV curve
-    voltages = _iv_series_lambert_v_from_i(
-        curs, IL, I0, Rs, Rsh, a, neg_v_limit)
+    device_voltages = _iv_series_lambert_v_from_i(
+        cur3, il, io, rs, rsh, a3, neg_v_limit)
+
+    # add voltage across devices to get string voltage
+    voltages = device_voltages.sum(axis=0)
+
+    # objective function for MPP
+    idx = np.arange(ntimes)
+    def mpp_optfn(current, idx):
+        # solve power at each current for each device
+        # current is ntimes only since it is common for all devices.
+        # other parameters are ntimes x ndevices
+        # broadcast current to ntimes x ndevices
+        cur2 = np.broadcast_to(current[np.newaxis, :], (ndevices, len(current)))
+        voltage = _iv_series_lambert_v_from_i(
+            cur2, IL, I0, Rs, Rsh, a, neg_v_limit, ndevices, idx)
+
+        # return negative of string power
+        return current * -voltage.sum(axis=0)
 
-    # add voltage across devices to get series voltage
-    voltage_sum = voltages.sum(axis=0)
+    # mpp calculations
+    idxs = np.argmax(currents * voltages, axis=1)
+    # reversed since currents is decreasing
+    idcs = idxs[:, np.newaxis] + [1, 0, -1]
+    intervals = np.take_along_axis(currents, idcs, axis=1)
+    init = (intervals[:, 0], intervals[:, 1], intervals[:, 2])
 
-    # drop currents dimension for devices
-    curs = curs[0, :, :]
+    idx = np.arange(ntimes)
+    imp_result = find_minimum(mpp_optfn, init, args=(idx,))
+    string_imp = imp_result.x
+    
+    # use string_imp to calculate string_vmp
+    idx = np.arange(ntimes)
+    cur2 = np.broadcast_to(string_imp[np.newaxis, :], (ndevices, len(string_imp)))
+    device_vmp = _iv_series_lambert_v_from_i(
+        cur2, IL, I0, Rs, Rsh, a, ndevices, neg_v_limit)
+    string_vmp = device_vmp.sum(axis=0)
 
-    return voltage_sum, curs
+    return voltages, currents, string_vmp, string_imp
diff --git a/tests/test_mismatch.py b/tests/test_mismatch.py
index 4c097df070..ae1f4da8d2 100644
--- a/tests/test_mismatch.py
+++ b/tests/test_mismatch.py
@@ -29,9 +29,9 @@ def test__iv_series_lambertw_isc():
     neg_v_limit = -5.
     expected_isc = _singlediode._lambertw_i_from_v(-neg_v_limit, IL[1], Io, Rs,
                                                    Rsh, nNsVth)
-    vs, cs = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth,
-                                 neg_v_limit=neg_v_limit,
-                                 npts=npts)
+    vs, cs, _, _ = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth,
+                                       neg_v_limit=neg_v_limit,
+                                       npts=npts)
     assert vs.shape == cs.shape
     assert vs.shape == (1, npts)  # ntimes x npts
     assert np.isclose(cs[0, 0], expected_isc)
@@ -51,9 +51,9 @@ def test__iv_series_lambertw_voc():
     neg_v_limit = -5.
     expected_voc = _singlediode._lambertw_v_from_i(0, IL, Io, Rs,
                                                    Rsh, nNsVth)
-    vs, cs = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth,
-                                 neg_v_limit=neg_v_limit,
-                                 npts=npts)
+    vs, cs, _, _ = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth,
+                                       neg_v_limit=neg_v_limit,
+                                       npts=npts)
 
     assert np.isclose(vs[0, -1], expected_voc.sum())
     assert np.isclose(cs[0, -1], 0.)
@@ -73,8 +73,8 @@ def test__iv_series_lambertw_breakpoints():
 
     cur_breakpoints = _singlediode._lambertw_i_from_v(neg_v_limit, IL, Io, Rs,
                                                       Rsh, nNsVth)
-    vs, cs = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth,
-                                 neg_v_limit=neg_v_limit,
-                                 npts=npts)
+    vs, cs, _, _ = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth,
+                                       neg_v_limit=neg_v_limit,
+                                       npts=npts)
     # current cs should contain all breakpoints less than string_isc
     assert np.isin(cur_breakpoints[cur_breakpoints < cs[0, 0]], cs).all()

From 7c16703fd309de119dab79cc07322e987f5d9487 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Wed, 26 Aug 2026 14:24:10 -0600
Subject: [PATCH 12/12] refactor to parallel singlediode, v_from_i functions

---
 pvlib/mismatch.py      | 130 ++++++++++++++++++++++++++++++++++++-----
 tests/test_mismatch.py |  65 ++++-----------------
 2 files changed, 125 insertions(+), 70 deletions(-)

diff --git a/pvlib/mismatch.py b/pvlib/mismatch.py
index e41294f51a..4b40692cc6 100644
--- a/pvlib/mismatch.py
+++ b/pvlib/mismatch.py
@@ -108,7 +108,24 @@ def _setup_currents(current_bkpts, string_isc, npts):
     return currents
 
 
-def _iv_series_lambertw(photocurrent, saturation_current, resistance_series,
+def _insert_mpp(vol, cur, vmp, imp):
+    # insert MPP into current and voltage arrays
+    cur_out = np.empty((vol.shape[0], vol.shape[1] + 1))
+    vol_out = np.empty((vol.shape[0], vol.shape[1] + 1))
+    # insert mpp value
+    for i in range(vol.shape[0]):
+        idx = vol[i, :].searchsorted(vmp[i], side='left')
+        vol_out[i, :idx] = vol[i, :idx]
+        vol_out[i, idx] = vmp[i]
+        vol_out[i, idx+1:] = vol[i, idx:]
+        cur_out[i, :idx] = cur[i, :idx]
+        cur_out[i, idx] = imp[i]
+        cur_out[i, idx+1:] = cur[i, idx:]
+
+    return vol_out, cur_out
+
+
+def _singlediode_mismatch(photocurrent, saturation_current, resistance_series,
                         resistance_shunt, nNsVth, neg_v_limit=0.,
                         npts=100):
     r'''Solve the IV curve for series-connected devices where each device
@@ -145,12 +162,14 @@ def _iv_series_lambertw(photocurrent, saturation_current, resistance_series,
 
     Returns
     -------
-    voltages : numeric
-        Voltage points for the series IV curves (V), shape
-        (times, npts).
-    currents : numeric
-        Current points for the series IV curves (A), shape
-        (times, npts).
+    string_isc : numeric
+        Short-circuit current for the string. [A]
+    string_voc : numeric
+        Open-circuit voltage for the string. [V]
+    string_imp : numeric
+        Current at maximum power for the string. [A]
+    string_vmp : numeric
+        Voltage at maximum power for the string. [V]
 
     '''
     # target shape is ndevices x ntimes
@@ -192,12 +211,17 @@ def isc_optfn(current, idx):
         (min_isc, max_isc), args=(idx,))
     string_isc = isc_result.x  # 1d in ntimes
 
+    # find Voc for string-level curves
+    device_voc = _singlediode._lambertw_v_from_i(0., IL, I0, Rs, Rsh, a)
+    string_voc = device_voc.sum(axis=0)
+
+    # prepare for MPP calculation by computing voltages on a grid of currents
     # discretize current from string_isc down to 0 at each time step
     # Include each device's current at neg_v_limit since these will be the
     # curvature breakpoints in the series IV curve.
-    # discretize current from max(Isc) down to 0. at each time step
-    # Include each device Isc (except the highest) so that the series IV curve
-    currents = _setup_currents(current_bkpts, string_isc, npts)
+    # Leave a space for inserting MPP
+    # currents is ntimes x npts-1, decreasing
+    currents = _setup_currents(current_bkpts, string_isc, npts-1)
 
     # shape all arrays to be ndevices x ntimes x ncurrents
     cur3 = np.repeat(currents[np.newaxis, :, :], ndevices, axis=0)
@@ -210,6 +234,7 @@ def isc_optfn(current, idx):
         cur3, il, io, rs, rsh, a3, neg_v_limit)
 
     # add voltage across devices to get string voltage
+    # voltages is ntimes x ncurrents
     voltages = device_voltages.sum(axis=0)
 
     # objective function for MPP
@@ -219,14 +244,15 @@ def mpp_optfn(current, idx):
         # current is ntimes only since it is common for all devices.
         # other parameters are ntimes x ndevices
         # broadcast current to ntimes x ndevices
-        cur2 = np.broadcast_to(current[np.newaxis, :], (ndevices, len(current)))
+        cur2 = np.broadcast_to(current[np.newaxis, :],
+                               (ndevices, len(current)))
         voltage = _iv_series_lambert_v_from_i(
             cur2, IL, I0, Rs, Rsh, a, neg_v_limit, ndevices, idx)
 
         # return negative of string power
         return current * -voltage.sum(axis=0)
 
-    # mpp calculations
+    # mpp calculation
     idxs = np.argmax(currents * voltages, axis=1)
     # reversed since currents is decreasing
     idcs = idxs[:, np.newaxis] + [1, 0, -1]
@@ -236,12 +262,84 @@ def mpp_optfn(current, idx):
     idx = np.arange(ntimes)
     imp_result = find_minimum(mpp_optfn, init, args=(idx,))
     string_imp = imp_result.x
-    
+
     # use string_imp to calculate string_vmp
     idx = np.arange(ntimes)
-    cur2 = np.broadcast_to(string_imp[np.newaxis, :], (ndevices, len(string_imp)))
+    cur2 = np.broadcast_to(string_imp[np.newaxis, :],
+                           (ndevices, ntimes))
     device_vmp = _iv_series_lambert_v_from_i(
-        cur2, IL, I0, Rs, Rsh, a, ndevices, neg_v_limit)
+        cur2, IL, I0, Rs, Rsh, a, neg_v_limit, ndevices)
     string_vmp = device_vmp.sum(axis=0)
 
-    return voltages, currents, string_vmp, string_imp
+    # in case we decide that this function should return (voltage, current)
+    voltages, currents = _insert_mpp(
+        voltages, currents, string_vmp, string_imp)
+
+    return string_isc, string_voc, string_imp, string_vmp
+    
+
+def _v_from_i_mismatch(currents, photocurrent, saturation_current,
+                       resistance_series,
+                       resistance_shunt, nNsVth, neg_v_limit=0.):
+    r'''Solve the IV curve for series-connected devices where each device
+    is described by the single diode equation.
+
+    Uses a simplified model for reverse bias behavior, where current is
+    unbounded at a constant reverse bias voltage ``neg_v_limit``.
+
+    Input parameter ``photocurrent`` must have shape (devices, times).
+    Input parameters ``saturation_current``, ``resistance_series``,
+    ``resistance_shunt``, ``nNsVth`` may be arrays. If arrays, must be
+    broadcastable to the shape of ``photocurrent``.
+
+    Parameters
+    ----------
+    currents : numeric
+        String current (A) at which string voltage is to be computed.
+        Must have shape (times, currents).
+    photocurrent : numeric
+        photocurrent (A). Must have shape (devices, times).
+    saturation_current : numeric
+        saturation current (A). Must be broadcastable with photocurrent.
+    resistance_series : numeric
+        series resistance (ohm). Must be broadcastable with photocurrent.
+    resistance_shunt : numeric
+        shunt resistance (ohm). Must be broadcastable with photocurrent.
+    nNsVth : numeric
+        product of diode factor n, number of series cells Ns, and
+        thermal voltage (Vth), (V). Must be broadcastable with photocurrent.
+    neg_v_limit : float, optional
+        Limit on reverse bias voltage, from cell breakdown voltage or reverse
+        bias diode activation voltage (V). Should be negative. For example,
+        if neg_v_limit=-5, then at V=-5 current is unbounded in the positive
+        direction.
+
+    Returns
+    -------
+    voltages : numeric
+        String voltage (V) at the current points, shape (times).
+
+    '''
+    # target shape is ndevices x ntimes
+    IL, I0, Rs, Rsh, a = \
+        np.broadcast_arrays(photocurrent, saturation_current,
+                            resistance_series, resistance_shunt, nNsVth)
+
+    ndevices, ntimes = IL.shape
+
+    # shape all arrays to be ndevices x ntimes x ncurrents
+    cur3 = np.repeat(currents[np.newaxis, :, :], ndevices, axis=0)
+    cur3, il, io, rs, rsh, a3 = np.broadcast_arrays(
+        cur3, IL[:, :, np.newaxis], I0[:, :, np.newaxis], Rs[:, :, np.newaxis],
+        Rsh[:, :, np.newaxis], a[:, :, np.newaxis])
+
+    # solve voltages at each current for each IV curve
+    # applies neg_v_limit
+    device_voltages = _iv_series_lambert_v_from_i(
+        cur3, il, io, rs, rsh, a3, neg_v_limit)
+
+    # sum voltage across devices to get string voltage
+    # voltages is ntimes x ncurrents
+    voltages = device_voltages.sum(axis=0)
+
+    return voltages
diff --git a/tests/test_mismatch.py b/tests/test_mismatch.py
index ae1f4da8d2..c43a6863af 100644
--- a/tests/test_mismatch.py
+++ b/tests/test_mismatch.py
@@ -1,6 +1,6 @@
 import numpy as np
 from pvlib import singlediode as _singlediode
-from pvlib.mismatch import _setup_currents, _iv_series_lambertw
+from pvlib.mismatch import _setup_currents, _singlediode_mismatch
 
 
 def test__setup_currents():
@@ -15,66 +15,23 @@ def test__setup_currents():
     assert curs[:, 0] == string_isc
 
 
-def test__iv_series_lambertw_isc():
+def test__singlediode_mismatch_isc_voc():
 
     # 2 devices, 1 time
     # Isc should be equal to the current on the higher curve where voltage is
     # +neg_v_limit
+    # Voc should be equal to the sum of Voc from each device
     IL = np.array([[1.0], [6.01]])
     Io = 1e-9
     nNsVth = 2.5
     Rs = 0.5
     Rsh = 1000.
-    npts = 5
     neg_v_limit = -5.
-    expected_isc = _singlediode._lambertw_i_from_v(-neg_v_limit, IL[1], Io, Rs,
-                                                   Rsh, nNsVth)
-    vs, cs, _, _ = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth,
-                                       neg_v_limit=neg_v_limit,
-                                       npts=npts)
-    assert vs.shape == cs.shape
-    assert vs.shape == (1, npts)  # ntimes x npts
-    assert np.isclose(cs[0, 0], expected_isc)
-    assert np.isclose(vs[0, 0], 0.)
-
-
-def test__iv_series_lambertw_voc():
-
-    # 2 devices, 1 time
-    # Voc should be equal to the sum of Voc for each device
-    IL = np.array([[1.0], [6.01]])
-    Io = 1e-9
-    nNsVth = 2.5
-    Rs = 0.5
-    Rsh = 1000.
-    npts = 5
-    neg_v_limit = -5.
-    expected_voc = _singlediode._lambertw_v_from_i(0, IL, Io, Rs,
-                                                   Rsh, nNsVth)
-    vs, cs, _, _ = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth,
-                                       neg_v_limit=neg_v_limit,
-                                       npts=npts)
-
-    assert np.isclose(vs[0, -1], expected_voc.sum())
-    assert np.isclose(cs[0, -1], 0.)
-
-
-def test__iv_series_lambertw_breakpoints():
-
-    # 2 devices, 1 time
-    # Voc should be equal to the sum of Voc for each device
-    IL = np.array([[1.0], [6.01]])
-    Io = 1e-9
-    nNsVth = 2.5
-    Rs = 0.5
-    Rsh = 1000.
-    npts = 10
-    neg_v_limit = -5.
-
-    cur_breakpoints = _singlediode._lambertw_i_from_v(neg_v_limit, IL, Io, Rs,
-                                                      Rsh, nNsVth)
-    vs, cs, _, _ = _iv_series_lambertw(IL, Io, Rs, Rsh, nNsVth,
-                                       neg_v_limit=neg_v_limit,
-                                       npts=npts)
-    # current cs should contain all breakpoints less than string_isc
-    assert np.isin(cur_breakpoints[cur_breakpoints < cs[0, 0]], cs).all()
+    expected_isc = _singlediode._lambertw_i_from_v(
+        -neg_v_limit, IL[1], Io, Rs, Rsh, nNsVth)
+    expected_voc = _singlediode._lambertw_v_from_i(
+        0, IL, Io, Rs, Rsh, nNsVth)
+    isc, voc, _, _ = _singlediode_mismatch(
+        IL, Io, Rs, Rsh, nNsVth, neg_v_limit=neg_v_limit)
+    assert np.isclose(isc[0], expected_isc)
+    assert np.isclose(voc[0], expected_voc.sum())