From c95ac2af95ee873832ffbf7eb9acf4103871452b Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Tue, 31 Mar 2020 14:14:53 -0600
Subject: [PATCH 01/12] add capability for reverse bias

---
 pvlib/singlediode.py            | 180 +++++++++++++++++++++++---------
 pvlib/tests/test_singlediode.py |  88 ++++++++++++++++
 2 files changed, 218 insertions(+), 50 deletions(-)

diff --git a/pvlib/singlediode.py b/pvlib/singlediode.py
index 707042bd7a..e489664f64 100644
--- a/pvlib/singlediode.py
+++ b/pvlib/singlediode.py
@@ -69,42 +69,65 @@ def estimate_voc(photocurrent, saturation_current, nNsVth):
 
 def bishop88(diode_voltage, photocurrent, saturation_current,
              resistance_series, resistance_shunt, nNsVth, d2mutau=0,
-             NsVbi=np.Inf, gradients=False):
-    """
+             NsVbi=np.Inf, breakdown_factor=0., breakdown_voltage=-5.5,
+             breakdown_exp=3.28, gradients=False):
+    r"""
     Explicit calculation of points on the IV curve described by the single
     diode equation.  Values are calculated as described in [1]_.
 
+    The single diode equation with recombination current and reverse bias
+    breakdown is
+
+    .. math::
+
+        I = I_{L} - I_{0} (\exp \frac{V_{d}}{nNsVth} - 1)
+        - \frac{V_{d}}{R_{sh}}
+        - \frac{I_{L} \frac{d^{2}}{\mu \tau}{N_{s} V_{bi} - V_{d}}
+        - a \frac{V_{d}{R_{sh}} (1 - \frac{V_{d}}{V_{br}})^-m
+
+    The input `diode_voltage` must be :math:`V + I R_{s}`.
+
+
     .. warning::
+       * Usage of ``d2mutau`` is required with PVSyst
+         coefficients for cadmium-telluride (CdTe) and amorphous-silicon
+         (a:Si) PV modules only.
        * Do not use ``d2mutau`` with CEC coefficients.
-       * Usage of ``d2mutau`` with PVSyst coefficients is required for cadmium-
-         telluride (CdTe) and amorphous-silicon (a:Si) PV modules only.
 
     Parameters
     ----------
     diode_voltage : numeric
         diode voltages [V]
     photocurrent : numeric
-        photo-generated current [A]
+        photo-generated current :math:`I_{L}` [A]
     saturation_current : numeric
-        diode reverse saturation current [A]
+        diode reverse saturation current :math:`I_{0}` [A]
     resistance_series : numeric
-        series resistance [ohms]
+        series resistance :math:`R_{s}` [ohms]
     resistance_shunt: numeric
-        shunt resistance [ohms]
+        shunt resistance :math:`R_{sh}` [ohms]
     nNsVth : numeric
-        product of thermal voltage ``Vth`` [V], diode ideality factor ``n``,
-        and number of series cells ``Ns``
+        product of thermal voltage :math:`V_{th}` [V], diode ideality factor
+        ``n``, and number of series cells :math:`N_{s}`
     d2mutau : numeric, default 0
         PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon
         (a-Si) modules that accounts for recombination current in the
         intrinsic layer. The value is the ratio of intrinsic layer thickness
         squared :math:`d^2` to the diffusion length of charge carriers
-        :math:`\\mu \\tau`. [V]
+        :math:`\mu \tau`. [V]
     NsVbi : numeric, default np.inf
         PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon
         (a-Si) modules that is the product of the PV module number of series
-        cells ``Ns`` and the builtin voltage ``Vbi`` of the intrinsic layer.
-        [V].
+        cells :math:`N_{s}` and the builtin voltage :math:`V_{bi}` of the
+        intrinsic layer. [V].
+    breakdown_factor : numeric, default 0
+        fraction of ohmic current involved in avalanche breakdown :math:`a`.
+        Default of 0 excludes the reverse bias terms from the model. [unitless]
+    breakdown_voltage : numeric, default -5.5
+        reverse breakdown voltage of the photovoltaic junction :math:`V_{br}`
+        [V]
+    breakdown_exp : numeric, default 3.28
+        avalanche breakdown exponent :math:`m` [unitless]
     gradients : bool
         False returns only I, V, and P. True also returns gradients
 
@@ -150,8 +173,14 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
     # calculate temporary values to simplify calculations
     v_star = diode_voltage / nNsVth  # non-dimensional diode voltage
     g_sh = 1.0 / resistance_shunt  # conductance
+    if breakdown_factor > 0:  # reverse bias is considered
+        brk_term = 1 - diode_voltage / breakdown_voltage
+        brk_pwr = np.power(brk_term, -breakdown_exp)
+        i_breakdown = breakdown_factor * diode_voltage * g_sh * brk_pwr
+    else:
+        i_breakdown = 0.
     i = (photocurrent - saturation_current * np.expm1(v_star)
-         - diode_voltage * g_sh - i_recomb)
+         - diode_voltage * g_sh - i_recomb - i_breakdown)
     v = diode_voltage - i * resistance_series
     retval = (i, v, i*v)
     if gradients:
@@ -159,12 +188,23 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
         grad_i_recomb = np.where(is_recomb, i_recomb / v_recomb, 0)
         grad_2i_recomb = np.where(is_recomb, 2 * grad_i_recomb / v_recomb, 0)
         g_diode = saturation_current * np.exp(v_star) / nNsVth  # conductance
-        grad_i = -g_diode - g_sh - grad_i_recomb  # di/dvd
+        if breakdown_factor > 0:  # reverse bias is considered
+            brk_pwr_1 = np.power(brk_term, -breakdown_exp - 1)
+            brk_pwr_2 = np.power(brk_term, -breakdown_exp - 2)
+            brk_fctr = breakdown_factor * g_sh
+            grad_i_brk = brk_fctr * (brk_pwr + diode_voltage *
+                                     -breakdown_exp * brk_pwr_1)
+            grad2i_brk = brk_fctr * -breakdown_exp * (2 * brk_pwr_1 + \
+                diode_voltage * (-breakdown_exp - 1) * brk_pwr_2)
+        else:
+            grad_i_brk = 0.
+            grad2i_brk = 0.
+        grad_i = -g_diode - g_sh - grad_i_recomb - grad_i_brk  # di/dvd
         grad_v = 1.0 - grad_i * resistance_series  # dv/dvd
         # dp/dv = d(iv)/dv = v * di/dv + i
         grad = grad_i / grad_v  # di/dv
         grad_p = v * grad + i  # dp/dv
-        grad2i = -g_diode / nNsVth - grad_2i_recomb  # d2i/dvd
+        grad2i = -g_diode / nNsVth - grad_2i_recomb - grad2i_brk  # d2i/dvd
         grad2v = -grad2i * resistance_series  # d2v/dvd
         grad2p = (
             grad_v * grad + v * (grad2i/grad_v - grad_i*grad2v/grad_v**2)
@@ -176,7 +216,9 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
 
 def bishop88_i_from_v(voltage, photocurrent, saturation_current,
                       resistance_series, resistance_shunt, nNsVth,
-                      d2mutau=0, NsVbi=np.Inf, method='newton'):
+                      d2mutau=0, NsVbi=np.Inf, breakdown_factor=1.e-4,
+                      breakdown_voltage=-5.5, breakdown_exp=3.28,
+                      method='newton'):
     """
     Find current given any voltage.
 
@@ -185,13 +227,13 @@ def bishop88_i_from_v(voltage, photocurrent, saturation_current,
     voltage : numeric
         voltage (V) in volts [V]
     photocurrent : numeric
-        photogenerated current (Iph or IL) in amperes [A]
+        photogenerated current (Iph or IL) [A]
     saturation_current : numeric
-        diode dark or saturation current (Io or Isat) in amperes [A]
+        diode dark or saturation current (Io or Isat) [A]
     resistance_series : numeric
-        series resistance (Rs) in ohms
+        series resistance (Rs) in [Ohm]
     resistance_shunt : numeric
-        shunt resistance (Rsh) in ohms
+        shunt resistance (Rsh) [Ohm]
     nNsVth : numeric
         product of diode ideality factor (n), number of series cells (Ns), and
         thermal voltage (Vth = k_b * T / q_e) in volts [V]
@@ -206,18 +248,27 @@ def bishop88_i_from_v(voltage, photocurrent, saturation_current,
         (a-Si) modules that is the product of the PV module number of series
         cells ``Ns`` and the builtin voltage ``Vbi`` of the intrinsic layer.
         [V].
-    method : str
-        one of two optional search methods: either ``'brentq'``, a reliable and
-        bounded method or ``'newton'`` which is the default.
+    breakdown_factor : numeric, default 1.0e-4
+        fraction of ohmic current involved in avalanche breakdown :math:`a`
+        [unitless]
+    breakdown_voltage : numeric, default -5.5
+        reverse breakdown voltage of the photovoltaic junction :math:`V_{br}`
+        [V]
+    breakdown_exp : numeric, default 3.28
+        avalanche breakdown exponent :math:`m` [unitless]
+    method : str, default 'newton'
+        must be either ``'newton'`` or ``'brentq'``. If ``'breakdown_factor'``
+        is not 0, must be ``'newton'``.
 
     Returns
     -------
     current : numeric
-        current (I) at the specified voltage (V) in amperes [A]
+        current (I) at the specified voltage (V). [A]
     """
     # collect args
     args = (photocurrent, saturation_current, resistance_series,
-            resistance_shunt, nNsVth, d2mutau, NsVbi)
+            resistance_shunt, nNsVth, d2mutau, NsVbi,
+            breakdown_factor, breakdown_voltage, breakdown_exp)
 
     def fv(x, v, *a):
         # calculate voltage residual given diode voltage "x"
@@ -230,9 +281,12 @@ def fv(x, v, *a):
         # brentq only works with scalar inputs, so we need a set up function
         # and np.vectorize to repeatedly call the optimizer with the right
         # arguments for possible array input
-        def vd_from_brent(voc, v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi):
+        def vd_from_brent(voc, v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
+                          breakdown_factor, breakdown_voltage, breakdown_exp):
             return brentq(fv, 0.0, voc,
-                          args=(v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi))
+                          args=(v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
+                                breakdown_factor, breakdown_voltage,
+                                breakdown_exp))
 
         vd_from_brent_vectorized = np.vectorize(vd_from_brent)
         vd = vd_from_brent_vectorized(voc_est, voltage, *args)
@@ -250,7 +304,9 @@ def vd_from_brent(voc, v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi):
 
 def bishop88_v_from_i(current, photocurrent, saturation_current,
                       resistance_series, resistance_shunt, nNsVth,
-                      d2mutau=0, NsVbi=np.Inf, method='newton'):
+                      d2mutau=0, NsVbi=np.Inf, breakdown_factor=1.e-4,
+                      breakdown_voltage=-5.5, breakdown_exp=3.28,
+                      method='newton'):
     """
     Find voltage given any current.
 
@@ -259,13 +315,13 @@ def bishop88_v_from_i(current, photocurrent, saturation_current,
     current : numeric
         current (I) in amperes [A]
     photocurrent : numeric
-        photogenerated current (Iph or IL) in amperes [A]
+        photogenerated current (Iph or IL) [A]
     saturation_current : numeric
-        diode dark or saturation current (Io or Isat) in amperes [A]
+        diode dark or saturation current (Io or Isat) [A]
     resistance_series : numeric
-        series resistance (Rs) in ohms
+        series resistance (Rs) in [Ohm]
     resistance_shunt : numeric
-        shunt resistance (Rsh) in ohms
+        shunt resistance (Rsh) [Ohm]
     nNsVth : numeric
         product of diode ideality factor (n), number of series cells (Ns), and
         thermal voltage (Vth = k_b * T / q_e) in volts [V]
@@ -280,9 +336,17 @@ def bishop88_v_from_i(current, photocurrent, saturation_current,
         (a-Si) modules that is the product of the PV module number of series
         cells ``Ns`` and the builtin voltage ``Vbi`` of the intrinsic layer.
         [V].
-    method : str
-        one of two optional search methods: either ``'brentq'``, a reliable and
-        bounded method or ``'newton'`` which is the default.
+    breakdown_factor : numeric, default 1.0e-4
+        fraction of ohmic current involved in avalanche breakdown :math:`a`
+        [unitless]
+    breakdown_voltage : numeric, default -5.5
+        reverse breakdown voltage of the photovoltaic junction :math:`V_{br}`
+        [V]
+    breakdown_exp : numeric, default 3.28
+        avalanche breakdown exponent :math:`m` [unitless]
+    method : str, default 'newton'
+        must be either ``'newton'`` or ``'brentq'``. If ``'breakdown_factor'``
+        is not 0, must be ``'newton'``.
 
     Returns
     -------
@@ -291,7 +355,8 @@ def bishop88_v_from_i(current, photocurrent, saturation_current,
     """
     # collect args
     args = (photocurrent, saturation_current, resistance_series,
-            resistance_shunt, nNsVth, d2mutau, NsVbi)
+            resistance_shunt, nNsVth, d2mutau, NsVbi, breakdown_factor,
+            breakdown_voltage, breakdown_exp)
     # first bound the search using voc
     voc_est = estimate_voc(photocurrent, saturation_current, nNsVth)
 
@@ -303,9 +368,12 @@ def fi(x, i, *a):
         # brentq only works with scalar inputs, so we need a set up function
         # and np.vectorize to repeatedly call the optimizer with the right
         # arguments for possible array input
-        def vd_from_brent(voc, i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi):
+        def vd_from_brent(voc, i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
+                          breakdown_factor, breakdown_voltage, breakdown_exp):
             return brentq(fi, 0.0, voc,
-                          args=(i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi))
+                          args=(i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
+                                breakdown_factor, breakdown_voltage,
+                                breakdown_exp))
 
         vd_from_brent_vectorized = np.vectorize(vd_from_brent)
         vd = vd_from_brent_vectorized(voc_est, current, *args)
@@ -323,20 +391,21 @@ def vd_from_brent(voc, i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi):
 
 def bishop88_mpp(photocurrent, saturation_current, resistance_series,
                  resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.Inf,
-                 method='newton'):
+                 breakdown_factor=1.e-4, breakdown_voltage=-5.5,
+                 breakdown_exp=3.28, method='newton'):
     """
     Find max power point.
 
     Parameters
     ----------
     photocurrent : numeric
-        photogenerated current (Iph or IL) in amperes [A]
+        photogenerated current (Iph or IL) [A]
     saturation_current : numeric
-        diode dark or saturation current (Io or Isat) in amperes [A]
+        diode dark or saturation current (Io or Isat) [A]
     resistance_series : numeric
-        series resistance (Rs) in ohms
+        series resistance (Rs) in [Ohm]
     resistance_shunt : numeric
-        shunt resistance (Rsh) in ohms
+        shunt resistance (Rsh) [Ohm]
     nNsVth : numeric
         product of diode ideality factor (n), number of series cells (Ns), and
         thermal voltage (Vth = k_b * T / q_e) in volts [V]
@@ -351,9 +420,17 @@ def bishop88_mpp(photocurrent, saturation_current, resistance_series,
         (a-Si) modules that is the product of the PV module number of series
         cells ``Ns`` and the builtin voltage ``Vbi`` of the intrinsic layer.
         [V].
-    method : str
-        one of two optional search methods: either ``'brentq'``, a reliable and
-        bounded method or ``'newton'`` which is the default.
+    breakdown_factor : numeric, default 1.0e-4
+        fraction of ohmic current involved in avalanche breakdown :math:`a`
+        [unitless]
+    breakdown_voltage : numeric, default -5.5
+        reverse breakdown voltage of the photovoltaic junction :math:`V_{br}`
+        [V]
+    breakdown_exp : numeric, default 3.28
+        avalanche breakdown exponent :math:`m` [unitless]
+    method : str, default 'newton'
+        must be either ``'newton'`` or ``'brentq'``. If ``'breakdown_factor'``
+        is not 0, must be ``'newton'``.
 
     Returns
     -------
@@ -363,7 +440,8 @@ def bishop88_mpp(photocurrent, saturation_current, resistance_series,
     """
     # collect args
     args = (photocurrent, saturation_current, resistance_series,
-            resistance_shunt, nNsVth, d2mutau, NsVbi)
+            resistance_shunt, nNsVth, d2mutau, NsVbi, breakdown_factor,
+            breakdown_voltage, breakdown_exp)
     # first bound the search using voc
     voc_est = estimate_voc(photocurrent, saturation_current, nNsVth)
 
@@ -373,9 +451,11 @@ def fmpp(x, *a):
     if method.lower() == 'brentq':
         # break out arguments for numpy.vectorize to handle broadcasting
         vec_fun = np.vectorize(
-            lambda voc, iph, isat, rs, rsh, gamma, d2mutau, NsVbi:
+            lambda voc, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
+                vbr_a, vbr, vbr_exp:
                 brentq(fmpp, 0.0, voc,
-                       args=(iph, isat, rs, rsh, gamma, d2mutau, NsVbi))
+                       args=(iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
+                             vbr_a, vbr, vbr_exp))
         )
         vd = vec_fun(voc_est, *args)
     elif method.lower() == 'newton':
diff --git a/pvlib/tests/test_singlediode.py b/pvlib/tests/test_singlediode.py
index 2d5074c884..41470a3762 100644
--- a/pvlib/tests/test_singlediode.py
+++ b/pvlib/tests/test_singlediode.py
@@ -185,3 +185,91 @@ def test_pvsyst_recombination_loss(method, poa, temp_cell, expected, tol):
 
     vsc_88 = bishop88_v_from_i(isc_88, *x, **y, method=method)
     assert np.isclose(vsc_88, 0.0, *tol)
+
+
+@requires_scipy
+@pytest.mark.parametrize(
+    'poa, temp_cell, expected, tol', [
+        # reference conditions
+        (
+            get_pvsyst_fs_495()['irrad_ref'],
+            get_pvsyst_fs_495()['temp_ref'],
+            {
+                'pmp': (get_pvsyst_fs_495()['I_mp_ref'] *
+                        get_pvsyst_fs_495()['V_mp_ref']),
+                'isc': get_pvsyst_fs_495()['I_sc_ref'],
+                'voc': get_pvsyst_fs_495()['V_oc_ref']
+            },
+            (5e-4, 0.04)
+        ),
+        # other conditions
+        (
+            POA,
+            TCELL,
+            {
+                'pmp': 79.723,
+                'isc': 1.4071,
+                'voc': 79.646
+            },
+            (1e-4, 1e-4)
+        )
+    ]
+)
+@pytest.mark.parametrize('method', ['newton', 'brentq'])
+def test_pvsyst_breakdown(method, poa, temp_cell, expected, tol):
+    """test PVSst recombination loss"""
+    pvsyst_fs_495 = get_pvsyst_fs_495()
+    # first evaluate PVSyst model with thin-film recombination loss current
+    # at reference conditions
+    x = pvsystem.calcparams_pvsyst(
+        effective_irradiance=poa, temp_cell=temp_cell,
+        alpha_sc=pvsyst_fs_495['alpha_sc'],
+        gamma_ref=pvsyst_fs_495['gamma_ref'],
+        mu_gamma=pvsyst_fs_495['mu_gamma'], I_L_ref=pvsyst_fs_495['I_L_ref'],
+        I_o_ref=pvsyst_fs_495['I_o_ref'], R_sh_ref=pvsyst_fs_495['R_sh_ref'],
+        R_sh_0=pvsyst_fs_495['R_sh_0'], R_sh_exp=pvsyst_fs_495['R_sh_exp'],
+        R_s=pvsyst_fs_495['R_s'],
+        cells_in_series=pvsyst_fs_495['cells_in_series'],
+        EgRef=pvsyst_fs_495['EgRef']
+    )
+    il_pvsyst, io_pvsyst, rs_pvsyst, rsh_pvsyst, nnsvt_pvsyst = x
+    voc_est_pvsyst = estimate_voc(photocurrent=il_pvsyst,
+                                  saturation_current=io_pvsyst,
+                                  nNsVth=nnsvt_pvsyst)
+    vd_pvsyst = np.linspace(0, voc_est_pvsyst, 1000)
+    pvsyst = bishop88(
+        diode_voltage=vd_pvsyst, photocurrent=il_pvsyst,
+        saturation_current=io_pvsyst, resistance_series=rs_pvsyst,
+        resistance_shunt=rsh_pvsyst, nNsVth=nnsvt_pvsyst,
+        breakdown_factor=1.-4, breakdown_voltage=-5.5,
+        breakdown_exp=3.28,
+    )
+    # test max power
+    assert np.isclose(max(pvsyst[2]), expected['pmp'], *tol)
+
+    # test short circuit current
+    isc_pvsyst = np.interp(0, pvsyst[1], pvsyst[0])
+    assert np.isclose(isc_pvsyst, expected['isc'], *tol)
+
+    # test open circuit voltage
+    voc_pvsyst = np.interp(0, pvsyst[0][::-1], pvsyst[1][::-1])
+    assert np.isclose(voc_pvsyst, expected['voc'], *tol)
+
+    # repeat tests as above with specialized bishop88 functions
+    y = dict(breakdown_factor=1.-4, breakdown_voltage=-5.5,
+             breakdown_exp=3.28)
+
+    mpp_88 = bishop88_mpp(*x, **y, method=method)
+    assert np.isclose(mpp_88[2], expected['pmp'], *tol)
+
+    isc_88 = bishop88_i_from_v(0, *x, **y, method=method)
+    assert np.isclose(isc_88, expected['isc'], *tol)
+
+    voc_88 = bishop88_v_from_i(0, *x, **y, method=method)
+    assert np.isclose(voc_88, expected['voc'], *tol)
+
+    ioc_88 = bishop88_i_from_v(voc_88, *x, **y, method=method)
+    assert np.isclose(ioc_88, 0.0, *tol)
+
+    vsc_88 = bishop88_v_from_i(isc_88, *x, **y, method=method)
+    assert np.isclose(vsc_88, 0.0, *tol)
\ No newline at end of file

From 4f07c0426a285dcf61f99017d160ee7d5c7c097b Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Tue, 31 Mar 2020 16:59:29 -0600
Subject: [PATCH 02/12] lint, correct test_pvsyst_breakdown

---
 pvlib/singlediode.py            | 18 +++++++++---------
 pvlib/tests/test_singlediode.py | 21 ++++++++++++---------
 2 files changed, 21 insertions(+), 18 deletions(-)

diff --git a/pvlib/singlediode.py b/pvlib/singlediode.py
index e489664f64..2bdd165a1b 100644
--- a/pvlib/singlediode.py
+++ b/pvlib/singlediode.py
@@ -179,8 +179,8 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
         i_breakdown = breakdown_factor * diode_voltage * g_sh * brk_pwr
     else:
         i_breakdown = 0.
-    i = (photocurrent - saturation_current * np.expm1(v_star)
-         - diode_voltage * g_sh - i_recomb - i_breakdown)
+    i = (photocurrent - saturation_current * np.expm1(v_star) -
+         diode_voltage * g_sh - i_recomb - i_breakdown)
     v = diode_voltage - i * resistance_series
     retval = (i, v, i*v)
     if gradients:
@@ -194,8 +194,9 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
             brk_fctr = breakdown_factor * g_sh
             grad_i_brk = brk_fctr * (brk_pwr + diode_voltage *
                                      -breakdown_exp * brk_pwr_1)
-            grad2i_brk = brk_fctr * -breakdown_exp * (2 * brk_pwr_1 + \
-                diode_voltage * (-breakdown_exp - 1) * brk_pwr_2)
+            grad2i_brk = brk_fctr * -breakdown_exp * (
+                2 * brk_pwr_1 + diode_voltage * (-breakdown_exp - 1) *
+                brk_pwr_2)
         else:
             grad_i_brk = 0.
             grad2i_brk = 0.
@@ -451,11 +452,10 @@ def fmpp(x, *a):
     if method.lower() == 'brentq':
         # break out arguments for numpy.vectorize to handle broadcasting
         vec_fun = np.vectorize(
-            lambda voc, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
-                vbr_a, vbr, vbr_exp:
-                brentq(fmpp, 0.0, voc,
-                       args=(iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
-                             vbr_a, vbr, vbr_exp))
+            lambda voc, iph, isat, rs, rsh, gamma, d2mutau, NsVbi, vbr_a, vbr,
+            vbr_exp: brentq(fmpp, 0.0, voc,
+                            args=(iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
+                            vbr_a, vbr, vbr_exp))
         )
         vd = vec_fun(voc_est, *args)
     elif method.lower() == 'newton':
diff --git a/pvlib/tests/test_singlediode.py b/pvlib/tests/test_singlediode.py
index 41470a3762..c538a99c1c 100644
--- a/pvlib/tests/test_singlediode.py
+++ b/pvlib/tests/test_singlediode.py
@@ -189,21 +189,23 @@ def test_pvsyst_recombination_loss(method, poa, temp_cell, expected, tol):
 
 @requires_scipy
 @pytest.mark.parametrize(
-    'poa, temp_cell, expected, tol', [
-        # reference conditions
+    'brk_params, poa, temp_cell, expected, tol', [
+        # reference conditions without breakdown model
         (
+            (0., -5.5, 3.28),
             get_pvsyst_fs_495()['irrad_ref'],
             get_pvsyst_fs_495()['temp_ref'],
             {
-                'pmp': (get_pvsyst_fs_495()['I_mp_ref'] *
+                'pmp': (get_pvsyst_fs_495()['I_mp_ref'] *  # noqa: W504
                         get_pvsyst_fs_495()['V_mp_ref']),
                 'isc': get_pvsyst_fs_495()['I_sc_ref'],
                 'voc': get_pvsyst_fs_495()['V_oc_ref']
             },
             (5e-4, 0.04)
         ),
-        # other conditions
+        # other conditions with breakdown model on
         (
+            (1.e-4, -5.5, 3.28),
             POA,
             TCELL,
             {
@@ -216,7 +218,7 @@ def test_pvsyst_recombination_loss(method, poa, temp_cell, expected, tol):
     ]
 )
 @pytest.mark.parametrize('method', ['newton', 'brentq'])
-def test_pvsyst_breakdown(method, poa, temp_cell, expected, tol):
+def test_pvsyst_breakdown(method, brk_params, poa, temp_cell, expected, tol):
     """test PVSst recombination loss"""
     pvsyst_fs_495 = get_pvsyst_fs_495()
     # first evaluate PVSyst model with thin-film recombination loss current
@@ -233,6 +235,9 @@ def test_pvsyst_breakdown(method, poa, temp_cell, expected, tol):
         EgRef=pvsyst_fs_495['EgRef']
     )
     il_pvsyst, io_pvsyst, rs_pvsyst, rsh_pvsyst, nnsvt_pvsyst = x
+
+    breakdown_factor, breakdown_voltage, breakdown_exp = brk_params
+
     voc_est_pvsyst = estimate_voc(photocurrent=il_pvsyst,
                                   saturation_current=io_pvsyst,
                                   nNsVth=nnsvt_pvsyst)
@@ -242,7 +247,7 @@ def test_pvsyst_breakdown(method, poa, temp_cell, expected, tol):
         saturation_current=io_pvsyst, resistance_series=rs_pvsyst,
         resistance_shunt=rsh_pvsyst, nNsVth=nnsvt_pvsyst,
         breakdown_factor=1.-4, breakdown_voltage=-5.5,
-        breakdown_exp=3.28,
+        breakdown_exp=3.28
     )
     # test max power
     assert np.isclose(max(pvsyst[2]), expected['pmp'], *tol)
@@ -256,8 +261,6 @@ def test_pvsyst_breakdown(method, poa, temp_cell, expected, tol):
     assert np.isclose(voc_pvsyst, expected['voc'], *tol)
 
     # repeat tests as above with specialized bishop88 functions
-    y = dict(breakdown_factor=1.-4, breakdown_voltage=-5.5,
-             breakdown_exp=3.28)
 
     mpp_88 = bishop88_mpp(*x, **y, method=method)
     assert np.isclose(mpp_88[2], expected['pmp'], *tol)
@@ -272,4 +275,4 @@ def test_pvsyst_breakdown(method, poa, temp_cell, expected, tol):
     assert np.isclose(ioc_88, 0.0, *tol)
 
     vsc_88 = bishop88_v_from_i(isc_88, *x, **y, method=method)
-    assert np.isclose(vsc_88, 0.0, *tol)
\ No newline at end of file
+    assert np.isclose(vsc_88, 0.0, *tol)

From 90692a7f05bbf5b0a3e79a818ba419b3714f6387 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Wed, 1 Apr 2020 10:16:45 -0600
Subject: [PATCH 03/12] default value for breakdown_factor

---
 pvlib/singlediode.py            | 26 +++++++++++++-------------
 pvlib/tests/test_singlediode.py | 10 +++++-----
 2 files changed, 18 insertions(+), 18 deletions(-)

diff --git a/pvlib/singlediode.py b/pvlib/singlediode.py
index 2bdd165a1b..c7529905eb 100644
--- a/pvlib/singlediode.py
+++ b/pvlib/singlediode.py
@@ -122,7 +122,7 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
         intrinsic layer. [V].
     breakdown_factor : numeric, default 0
         fraction of ohmic current involved in avalanche breakdown :math:`a`.
-        Default of 0 excludes the reverse bias terms from the model. [unitless]
+        Default of 0 excludes the reverse bias term from the model. [unitless]
     breakdown_voltage : numeric, default -5.5
         reverse breakdown voltage of the photovoltaic junction :math:`V_{br}`
         [V]
@@ -217,7 +217,7 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
 
 def bishop88_i_from_v(voltage, photocurrent, saturation_current,
                       resistance_series, resistance_shunt, nNsVth,
-                      d2mutau=0, NsVbi=np.Inf, breakdown_factor=1.e-4,
+                      d2mutau=0, NsVbi=np.Inf, breakdown_factor=0.,
                       breakdown_voltage=-5.5, breakdown_exp=3.28,
                       method='newton'):
     """
@@ -249,9 +249,9 @@ def bishop88_i_from_v(voltage, photocurrent, saturation_current,
         (a-Si) modules that is the product of the PV module number of series
         cells ``Ns`` and the builtin voltage ``Vbi`` of the intrinsic layer.
         [V].
-    breakdown_factor : numeric, default 1.0e-4
-        fraction of ohmic current involved in avalanche breakdown :math:`a`
-        [unitless]
+    breakdown_factor : numeric, default 0
+        fraction of ohmic current involved in avalanche breakdown :math:`a`.
+        Default of 0 excludes the reverse bias term from the model. [unitless]
     breakdown_voltage : numeric, default -5.5
         reverse breakdown voltage of the photovoltaic junction :math:`V_{br}`
         [V]
@@ -305,7 +305,7 @@ def vd_from_brent(voc, v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
 
 def bishop88_v_from_i(current, photocurrent, saturation_current,
                       resistance_series, resistance_shunt, nNsVth,
-                      d2mutau=0, NsVbi=np.Inf, breakdown_factor=1.e-4,
+                      d2mutau=0, NsVbi=np.Inf, breakdown_factor=0.,
                       breakdown_voltage=-5.5, breakdown_exp=3.28,
                       method='newton'):
     """
@@ -337,9 +337,9 @@ def bishop88_v_from_i(current, photocurrent, saturation_current,
         (a-Si) modules that is the product of the PV module number of series
         cells ``Ns`` and the builtin voltage ``Vbi`` of the intrinsic layer.
         [V].
-    breakdown_factor : numeric, default 1.0e-4
-        fraction of ohmic current involved in avalanche breakdown :math:`a`
-        [unitless]
+    breakdown_factor : numeric, default 0
+        fraction of ohmic current involved in avalanche breakdown :math:`a`.
+        Default of 0 excludes the reverse bias term from the model. [unitless]
     breakdown_voltage : numeric, default -5.5
         reverse breakdown voltage of the photovoltaic junction :math:`V_{br}`
         [V]
@@ -392,7 +392,7 @@ def vd_from_brent(voc, i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
 
 def bishop88_mpp(photocurrent, saturation_current, resistance_series,
                  resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.Inf,
-                 breakdown_factor=1.e-4, breakdown_voltage=-5.5,
+                 breakdown_factor=0., breakdown_voltage=-5.5,
                  breakdown_exp=3.28, method='newton'):
     """
     Find max power point.
@@ -421,9 +421,9 @@ def bishop88_mpp(photocurrent, saturation_current, resistance_series,
         (a-Si) modules that is the product of the PV module number of series
         cells ``Ns`` and the builtin voltage ``Vbi`` of the intrinsic layer.
         [V].
-    breakdown_factor : numeric, default 1.0e-4
-        fraction of ohmic current involved in avalanche breakdown :math:`a`
-        [unitless]
+    breakdown_factor : numeric, default 0
+        fraction of ohmic current involved in avalanche breakdown :math:`a`.
+        Default of 0 excludes the reverse bias term from the model. [unitless]
     breakdown_voltage : numeric, default -5.5
         reverse breakdown voltage of the photovoltaic junction :math:`V_{br}`
         [V]
diff --git a/pvlib/tests/test_singlediode.py b/pvlib/tests/test_singlediode.py
index c538a99c1c..1768e5f3e5 100644
--- a/pvlib/tests/test_singlediode.py
+++ b/pvlib/tests/test_singlediode.py
@@ -262,17 +262,17 @@ def test_pvsyst_breakdown(method, brk_params, poa, temp_cell, expected, tol):
 
     # repeat tests as above with specialized bishop88 functions
 
-    mpp_88 = bishop88_mpp(*x, **y, method=method)
+    mpp_88 = bishop88_mpp(*x, **brk_params, method=method)
     assert np.isclose(mpp_88[2], expected['pmp'], *tol)
 
-    isc_88 = bishop88_i_from_v(0, *x, **y, method=method)
+    isc_88 = bishop88_i_from_v(0, *x, **brk_params, method=method)
     assert np.isclose(isc_88, expected['isc'], *tol)
 
-    voc_88 = bishop88_v_from_i(0, *x, **y, method=method)
+    voc_88 = bishop88_v_from_i(0, *x, **brk_params, method=method)
     assert np.isclose(voc_88, expected['voc'], *tol)
 
-    ioc_88 = bishop88_i_from_v(voc_88, *x, **y, method=method)
+    ioc_88 = bishop88_i_from_v(voc_88, *x, **brk_params, method=method)
     assert np.isclose(ioc_88, 0.0, *tol)
 
-    vsc_88 = bishop88_v_from_i(isc_88, *x, **y, method=method)
+    vsc_88 = bishop88_v_from_i(isc_88, *x, **brk_params, method=method)
     assert np.isclose(vsc_88, 0.0, *tol)

From 70553fa068d30fc9b5f3568f714f58059522e9d6 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Wed, 1 Apr 2020 11:32:15 -0600
Subject: [PATCH 04/12] handle recombination in breakdown test, lint

---
 pvlib/singlediode.py            |  2 +-
 pvlib/tests/test_singlediode.py | 16 +++++++++++-----
 2 files changed, 12 insertions(+), 6 deletions(-)

diff --git a/pvlib/singlediode.py b/pvlib/singlediode.py
index c7529905eb..b9d18bf4ef 100644
--- a/pvlib/singlediode.py
+++ b/pvlib/singlediode.py
@@ -455,7 +455,7 @@ def fmpp(x, *a):
             lambda voc, iph, isat, rs, rsh, gamma, d2mutau, NsVbi, vbr_a, vbr,
             vbr_exp: brentq(fmpp, 0.0, voc,
                             args=(iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
-                            vbr_a, vbr, vbr_exp))
+                                  vbr_a, vbr, vbr_exp))
         )
         vd = vec_fun(voc_est, *args)
     elif method.lower() == 'newton':
diff --git a/pvlib/tests/test_singlediode.py b/pvlib/tests/test_singlediode.py
index 1768e5f3e5..56bc3f1151 100644
--- a/pvlib/tests/test_singlediode.py
+++ b/pvlib/tests/test_singlediode.py
@@ -189,10 +189,12 @@ def test_pvsyst_recombination_loss(method, poa, temp_cell, expected, tol):
 
 @requires_scipy
 @pytest.mark.parametrize(
-    'brk_params, poa, temp_cell, expected, tol', [
+    'brk_params, recomb_params, poa, temp_cell, expected, tol', [
         # reference conditions without breakdown model
         (
             (0., -5.5, 3.28),
+            (get_pvsyst_fs_495()['d2mutau'],
+             VOLTAGE_BUILTIN * get_pvsyst_fs_495()['cells_in_series']),
             get_pvsyst_fs_495()['irrad_ref'],
             get_pvsyst_fs_495()['temp_ref'],
             {
@@ -203,9 +205,10 @@ def test_pvsyst_recombination_loss(method, poa, temp_cell, expected, tol):
             },
             (5e-4, 0.04)
         ),
-        # other conditions with breakdown model on
+        # other conditions with breakdown model on and recombination model off
         (
             (1.e-4, -5.5, 3.28),
+            (0., np.Inf),
             POA,
             TCELL,
             {
@@ -218,7 +221,8 @@ def test_pvsyst_recombination_loss(method, poa, temp_cell, expected, tol):
     ]
 )
 @pytest.mark.parametrize('method', ['newton', 'brentq'])
-def test_pvsyst_breakdown(method, brk_params, poa, temp_cell, expected, tol):
+def test_pvsyst_breakdown(method, brk_params, recomb_params, poa, temp_cell,
+                          expected, tol):
     """test PVSst recombination loss"""
     pvsyst_fs_495 = get_pvsyst_fs_495()
     # first evaluate PVSyst model with thin-film recombination loss current
@@ -236,6 +240,7 @@ def test_pvsyst_breakdown(method, brk_params, poa, temp_cell, expected, tol):
     )
     il_pvsyst, io_pvsyst, rs_pvsyst, rsh_pvsyst, nnsvt_pvsyst = x
 
+    d2mutau, NsVbi = recomb_params
     breakdown_factor, breakdown_voltage, breakdown_exp = brk_params
 
     voc_est_pvsyst = estimate_voc(photocurrent=il_pvsyst,
@@ -246,8 +251,9 @@ def test_pvsyst_breakdown(method, brk_params, poa, temp_cell, expected, tol):
         diode_voltage=vd_pvsyst, photocurrent=il_pvsyst,
         saturation_current=io_pvsyst, resistance_series=rs_pvsyst,
         resistance_shunt=rsh_pvsyst, nNsVth=nnsvt_pvsyst,
-        breakdown_factor=1.-4, breakdown_voltage=-5.5,
-        breakdown_exp=3.28
+        d2mutau=d2mutau, NsVbi=NsVbi,
+        breakdown_factor=breakdown_factor, breakdown_voltage=breakdown_voltage,
+        breakdown_exp=breakdown_exp
     )
     # test max power
     assert np.isclose(max(pvsyst[2]), expected['pmp'], *tol)

From 4bb530c591364e7a588a9f8910298ef287d0670d Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Wed, 1 Apr 2020 12:29:16 -0600
Subject: [PATCH 05/12] correct parameter parsing

---
 pvlib/tests/test_singlediode.py | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/pvlib/tests/test_singlediode.py b/pvlib/tests/test_singlediode.py
index 56bc3f1151..b7acaa4658 100644
--- a/pvlib/tests/test_singlediode.py
+++ b/pvlib/tests/test_singlediode.py
@@ -267,18 +267,21 @@ def test_pvsyst_breakdown(method, brk_params, recomb_params, poa, temp_cell,
     assert np.isclose(voc_pvsyst, expected['voc'], *tol)
 
     # repeat tests as above with specialized bishop88 functions
+    y = {'d2mutau': recomb_params[0], 'NsVbi': recomb_params[1],
+         'breakdown_factor': brk_params[0], 'breakdown_voltage': brk_params[1],
+         'breakdown_exp': brk_params[2]}
 
-    mpp_88 = bishop88_mpp(*x, **brk_params, method=method)
+    mpp_88 = bishop88_mpp(*x, **y, method=method)
     assert np.isclose(mpp_88[2], expected['pmp'], *tol)
 
-    isc_88 = bishop88_i_from_v(0, *x, **brk_params, method=method)
+    isc_88 = bishop88_i_from_v(0, *x, **y, method=method)
     assert np.isclose(isc_88, expected['isc'], *tol)
 
-    voc_88 = bishop88_v_from_i(0, *x, **brk_params, method=method)
+    voc_88 = bishop88_v_from_i(0, *x, **y, method=method)
     assert np.isclose(voc_88, expected['voc'], *tol)
 
-    ioc_88 = bishop88_i_from_v(voc_88, *x, **brk_params, method=method)
+    ioc_88 = bishop88_i_from_v(voc_88, *x, **y, method=method)
     assert np.isclose(ioc_88, 0.0, *tol)
 
-    vsc_88 = bishop88_v_from_i(isc_88, *x, **brk_params, method=method)
+    vsc_88 = bishop88_v_from_i(isc_88, *x, **y, method=method)
     assert np.isclose(vsc_88, 0.0, *tol)

From 4bcb4683937ceaa66c8cefb583b3f97b9beddea4 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Wed, 1 Apr 2020 13:04:12 -0600
Subject: [PATCH 06/12] whatsnew

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

diff --git a/docs/sphinx/source/whatsnew/v0.7.2.rst b/docs/sphinx/source/whatsnew/v0.7.2.rst
index 41bec86778..cfd77eec69 100644
--- a/docs/sphinx/source/whatsnew/v0.7.2.rst
+++ b/docs/sphinx/source/whatsnew/v0.7.2.rst
@@ -25,6 +25,10 @@ Enhancements
 * Add new module :py:mod:`pvlib.snow` to contain models related to snow coverage and effects on a PV system. (:pull:`764`)
 * Add snow coverage model :py:func:`pvlib.snow.coverage_nrel` and function to identify when modules are fully covered by snow :py:func:`pvlib.snow.fully_covered_nrel`. (:issue:`577`)
 * Add function :py:func:`pvlib.snow.dc_loss_nrel` for effect of snow coverage on DC output. (:pull:`764`)
+* Add capability to calculate current at reverse bias using an avalanche
+  breakdown model, affects :py:func:`pvlib.singlediode.bishop88`,
+  :py:func:`pvlib.singlediode.bishop88_i_from_v`, :py:func:`pvlib.singlediode.bishop88_v_from_i`,
+  :py:func:`pvlib.singlediode.bishop88_mpp`. (:pull:`948`)
 
 Bug fixes
 ~~~~~~~~~

From bcb93bf0b5edf104bcf82fee909ba6d93bdc1232 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Mon, 6 Apr 2020 09:19:26 -0600
Subject: [PATCH 07/12] Update pvlib/tests/test_singlediode.py

Co-Authored-By: Cameron Stark <CameronTStark@users.noreply.github.com>
---
 pvlib/tests/test_singlediode.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pvlib/tests/test_singlediode.py b/pvlib/tests/test_singlediode.py
index b7acaa4658..ca2a2e1270 100644
--- a/pvlib/tests/test_singlediode.py
+++ b/pvlib/tests/test_singlediode.py
@@ -223,7 +223,7 @@ def test_pvsyst_recombination_loss(method, poa, temp_cell, expected, tol):
 @pytest.mark.parametrize('method', ['newton', 'brentq'])
 def test_pvsyst_breakdown(method, brk_params, recomb_params, poa, temp_cell,
                           expected, tol):
-    """test PVSst recombination loss"""
+    """test PVSyst recombination loss"""
     pvsyst_fs_495 = get_pvsyst_fs_495()
     # first evaluate PVSyst model with thin-film recombination loss current
     # at reference conditions

From 32f80171e9aa95e8e5ee0ec18709771fbe7ca9f9 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Mon, 6 Apr 2020 09:22:03 -0600
Subject: [PATCH 08/12] Update pvlib/singlediode.py

Co-Authored-By: Cameron Stark <CameronTStark@users.noreply.github.com>
---
 pvlib/singlediode.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/pvlib/singlediode.py b/pvlib/singlediode.py
index b9d18bf4ef..fd14350237 100644
--- a/pvlib/singlediode.py
+++ b/pvlib/singlediode.py
@@ -195,7 +195,12 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
             grad_i_brk = brk_fctr * (brk_pwr + diode_voltage *
                                      -breakdown_exp * brk_pwr_1)
             grad2i_brk = brk_fctr * -breakdown_exp * (
-                2 * brk_pwr_1 + diode_voltage * (-breakdown_exp - 1) *
+            grad2i_brk = (
+                brk_fctr
+                * -breakdown_exp
+                * (2 * brk_pwr_1 + diode_voltage * (-breakdown_exp - 1) * brk_pwr_2)
+            )
+
                 brk_pwr_2)
         else:
             grad_i_brk = 0.

From 9dbfdf20ac5f3859f8e4588effcdd7029666cab4 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Mon, 6 Apr 2020 09:22:22 -0600
Subject: [PATCH 09/12] Update pvlib/singlediode.py

Co-Authored-By: Cameron Stark <CameronTStark@users.noreply.github.com>
---
 pvlib/singlediode.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/pvlib/singlediode.py b/pvlib/singlediode.py
index fd14350237..c24854a504 100644
--- a/pvlib/singlediode.py
+++ b/pvlib/singlediode.py
@@ -179,7 +179,12 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
         i_breakdown = breakdown_factor * diode_voltage * g_sh * brk_pwr
     else:
         i_breakdown = 0.
-    i = (photocurrent - saturation_current * np.expm1(v_star) -
+    i = (
+        photocurrent
+        - saturation_current * np.expm1(v_star)
+        - diode_voltage * g_sh
+        - i_recomb
+    )
          diode_voltage * g_sh - i_recomb - i_breakdown)
     v = diode_voltage - i * resistance_series
     retval = (i, v, i*v)

From 0ddde4c7cfb096fc8329dc3840773175065dcba3 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Mon, 6 Apr 2020 09:36:11 -0600
Subject: [PATCH 10/12] fix formatting

---
 pvlib/singlediode.py | 12 ++++--------
 1 file changed, 4 insertions(+), 8 deletions(-)

diff --git a/pvlib/singlediode.py b/pvlib/singlediode.py
index c24854a504..f2d2bc54fa 100644
--- a/pvlib/singlediode.py
+++ b/pvlib/singlediode.py
@@ -181,11 +181,9 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
         i_breakdown = 0.
     i = (
         photocurrent
-        - saturation_current * np.expm1(v_star)
-        - diode_voltage * g_sh
-        - i_recomb
+        - saturation_current * np.expm1(v_star) - diode_voltage * g_sh
+        - i_recomb - i_breakdown
     )
-         diode_voltage * g_sh - i_recomb - i_breakdown)
     v = diode_voltage - i * resistance_series
     retval = (i, v, i*v)
     if gradients:
@@ -199,14 +197,12 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
             brk_fctr = breakdown_factor * g_sh
             grad_i_brk = brk_fctr * (brk_pwr + diode_voltage *
                                      -breakdown_exp * brk_pwr_1)
-            grad2i_brk = brk_fctr * -breakdown_exp * (
             grad2i_brk = (
                 brk_fctr
                 * -breakdown_exp
-                * (2 * brk_pwr_1 + diode_voltage * (-breakdown_exp - 1) * brk_pwr_2)
+                * (2 * brk_pwr_1 + diode_voltage * (-breakdown_exp - 1)
+                * brk_pwr_2)
             )
-
-                brk_pwr_2)
         else:
             grad_i_brk = 0.
             grad2i_brk = 0.

From 21a7edc5f0b18770479b1ec8e808f60d50bee525 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Mon, 6 Apr 2020 10:19:18 -0600
Subject: [PATCH 11/12] improve docstring

---
 pvlib/singlediode.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/pvlib/singlediode.py b/pvlib/singlediode.py
index f2d2bc54fa..97569b71a8 100644
--- a/pvlib/singlediode.py
+++ b/pvlib/singlediode.py
@@ -264,8 +264,8 @@ def bishop88_i_from_v(voltage, photocurrent, saturation_current,
     breakdown_exp : numeric, default 3.28
         avalanche breakdown exponent :math:`m` [unitless]
     method : str, default 'newton'
-        must be either ``'newton'`` or ``'brentq'``. If ``'breakdown_factor'``
-        is not 0, must be ``'newton'``.
+       Either ``'newton'`` or ``'brentq'``. ''method'' must be ``'newton'``
+       if ``breakdown_factor`` is not 0.
 
     Returns
     -------
@@ -352,8 +352,8 @@ def bishop88_v_from_i(current, photocurrent, saturation_current,
     breakdown_exp : numeric, default 3.28
         avalanche breakdown exponent :math:`m` [unitless]
     method : str, default 'newton'
-        must be either ``'newton'`` or ``'brentq'``. If ``'breakdown_factor'``
-        is not 0, must be ``'newton'``.
+       Either ``'newton'`` or ``'brentq'``. ''method'' must be ``'newton'``
+       if ``breakdown_factor`` is not 0.
 
     Returns
     -------
@@ -436,8 +436,8 @@ def bishop88_mpp(photocurrent, saturation_current, resistance_series,
     breakdown_exp : numeric, default 3.28
         avalanche breakdown exponent :math:`m` [unitless]
     method : str, default 'newton'
-        must be either ``'newton'`` or ``'brentq'``. If ``'breakdown_factor'``
-        is not 0, must be ``'newton'``.
+       Either ``'newton'`` or ``'brentq'``. ''method'' must be ``'newton'``
+       if ``breakdown_factor`` is not 0.
 
     Returns
     -------

From 9445d32541c2c86f6d12e925daf1a680ec7e139e Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Mon, 6 Apr 2020 10:45:52 -0600
Subject: [PATCH 12/12] lint

---
 pvlib/singlediode.py | 16 +++++-----------
 1 file changed, 5 insertions(+), 11 deletions(-)

diff --git a/pvlib/singlediode.py b/pvlib/singlediode.py
index 97569b71a8..0576ee6cbd 100644
--- a/pvlib/singlediode.py
+++ b/pvlib/singlediode.py
@@ -179,11 +179,8 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
         i_breakdown = breakdown_factor * diode_voltage * g_sh * brk_pwr
     else:
         i_breakdown = 0.
-    i = (
-        photocurrent
-        - saturation_current * np.expm1(v_star) - diode_voltage * g_sh
-        - i_recomb - i_breakdown
-    )
+    i = (photocurrent - saturation_current * np.expm1(v_star)  # noqa: W503
+         - diode_voltage * g_sh - i_recomb - i_breakdown)   # noqa: W503
     v = diode_voltage - i * resistance_series
     retval = (i, v, i*v)
     if gradients:
@@ -197,12 +194,9 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
             brk_fctr = breakdown_factor * g_sh
             grad_i_brk = brk_fctr * (brk_pwr + diode_voltage *
                                      -breakdown_exp * brk_pwr_1)
-            grad2i_brk = (
-                brk_fctr
-                * -breakdown_exp
-                * (2 * brk_pwr_1 + diode_voltage * (-breakdown_exp - 1)
-                * brk_pwr_2)
-            )
+            grad2i_brk = (brk_fctr * -breakdown_exp        # noqa: W503
+                          * (2 * brk_pwr_1 + diode_voltage   # noqa: W503
+                             * (-breakdown_exp - 1) * brk_pwr_2))  # noqa: W503
         else:
             grad_i_brk = 0.
             grad2i_brk = 0.