singlediode.py

"""
Low-level functions for solving the single diode equation.
"""

import numpy as np
import pandas as pd
from pvlib.tools import _golden_sect_DataFrame

from scipy.optimize import brentq, newton
from scipy.special import lambertw

# newton method default parameters for this module
NEWTON_DEFAULT_PARAMS = {
    'tol': 1e-6,
    'maxiter': 100
}

# intrinsic voltage per cell junction for a:Si, CdTe, Mertens et al.
VOLTAGE_BUILTIN = 0.9  # [V]


def estimate_voc(photocurrent, saturation_current, nNsVth):
    """
    Rough estimate of open circuit voltage useful for bounding searches for
    ``i`` of ``v`` when using :func:`~pvlib.pvsystem.singlediode`.

    Parameters
    ----------
    photocurrent : numeric
        photo-generated current [A]
    saturation_current : numeric
        diode reverse saturation current [A]
    nNsVth : numeric
        product of thermal voltage ``Vth`` [V], diode ideality factor ``n``,
        and number of series cells ``Ns``

    Returns
    -------
    numeric
        rough estimate of open circuit voltage [V]

    Notes
    -----
    Calculating the open circuit voltage, :math:`V_{oc}`, of an ideal device
    with infinite shunt resistance, :math:`R_{sh} \\to \\infty`, and zero
    series resistance, :math:`R_s = 0`, yields the following equation [1]. As
    an estimate of :math:`V_{oc}` it is useful as an upper bound for the
    bisection method.

    .. math::

        V_{oc, est}=n Ns V_{th} \\log \\left( \\frac{I_L}{I_0} + 1 \\right)

    .. [1] http://www.pveducation.org/pvcdrom/open-circuit-voltage
    """

    return nNsVth * np.log(np.asarray(photocurrent) / saturation_current + 1.0)


def bishop88(diode_voltage, photocurrent, saturation_current,
             resistance_series, resistance_shunt, nNsVth, d2mutau=0,
             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} \left (\exp \frac{V_{d}}{nN_{s}V_{th}} - 1 \right )
        - \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}} \left (1 - \frac{V_{d}}{V_{br}} \right )^{-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.

    Parameters
    ----------
    diode_voltage : numeric
        diode voltage :math:`V_d` [V]
    photocurrent : numeric
        photo-generated current :math:`I_{L}` [A]
    saturation_current : numeric
        diode reverse saturation current :math:`I_{0}` [A]
    resistance_series : numeric
        series resistance :math:`R_{s}` [ohms]
    resistance_shunt: numeric
        shunt resistance :math:`R_{sh}` [ohms]
    nNsVth : numeric
        product of thermal voltage :math:`V_{th}` [V], diode ideality factor
        :math:`n`, and number of series cells :math:`N_{s}` [V]
    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]
    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 :math:`N_{s}` and the builtin voltage :math:`V_{bi}` of the
        intrinsic layer. [V].
    breakdown_factor : float, 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 : float, default -5.5
        reverse breakdown voltage of the photovoltaic junction :math:`V_{br}`
        [V]
    breakdown_exp : float, default 3.28
        avalanche breakdown exponent :math:`m` [unitless]
    gradients : bool
        False returns only I, V, and P. True also returns gradients

    Returns
    -------
    tuple
        currents [A], voltages [V], power [W], and optionally
        :math:`\frac{dI}{dV_d}`, :math:`\frac{dV}{dV_d}`,
        :math:`\frac{dI}{dV}`, :math:`\frac{dP}{dV}`, and
        :math:`\frac{d^2 P}{dV dV_d}`

    Notes
    -----
    The PVSyst thin-film recombination losses parameters ``d2mutau`` and
    ``NsVbi`` should only be applied to cadmium-telluride (CdTe) and amorphous-
    silicon (a-Si) PV modules, [2]_, [3]_. The builtin voltage :math:`V_{bi}`
    should account for all junctions. For example: tandem and triple junction
    cells would have builtin voltages of 1.8[V] and 2.7[V] respectively, based
    on the default of 0.9[V] for a single junction. The parameter ``NsVbi``
    should only account for the number of series cells in a single parallel
    sub-string if the module has cells in parallel greater than 1.

    References
    ----------
    .. [1] "Computer simulation of the effects of electrical mismatches in
       photovoltaic cell interconnection circuits" JW Bishop, Solar Cell (1988)
       :doi:`10.1016/0379-6787(88)90059-2`

    .. [2] "Improved equivalent circuit and Analytical Model for Amorphous
       Silicon Solar Cells and Modules." J. Mertens, et al., IEEE Transactions
       on Electron Devices, Vol 45, No 2, Feb 1998.
       :doi:`10.1109/16.658676`

    .. [3] "Performance assessment of a simulation model for PV modules of any
       available technology", André Mermoud and Thibault Lejeune, 25th EUPVSEC,
       2010
       :doi:`10.4229/25thEUPVSEC2010-4BV.1.114`
    """
    # calculate recombination loss current where d2mutau > 0
    is_recomb = d2mutau > 0  # True where there is thin-film recombination loss
    v_recomb = np.where(is_recomb, NsVbi - diode_voltage, np.inf)
    i_recomb = np.where(is_recomb, photocurrent * d2mutau / v_recomb, 0)
    # 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)  # 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:
        # calculate recombination loss current gradients where d2mutau > 0
        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
        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        # 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.
        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 - grad2i_brk  # d2i/dvd
        grad2v = -grad2i * resistance_series  # d2v/dvd
        grad2p = (
            grad_v * grad + v * (grad2i/grad_v - grad_i*grad2v/grad_v**2)
            + grad_i
        )  # d2p/dv/dvd
        retval += (grad_i, grad_v, grad, grad_p, grad2p)
    return retval


def bishop88_i_from_v(voltage, photocurrent, saturation_current,
                      resistance_series, resistance_shunt, nNsVth,
                      d2mutau=0, NsVbi=np.inf, breakdown_factor=0.,
                      breakdown_voltage=-5.5, breakdown_exp=3.28,
                      method='newton', method_kwargs=None):
    """
    Find current given any voltage.

    Parameters
    ----------
    voltage : numeric
        voltage (V) in volts [V]
    photocurrent : numeric
        photogenerated current (Iph or IL) [A]
    saturation_current : numeric
        diode dark or saturation current (Io or Isat) [A]
    resistance_series : numeric
        series resistance (Rs) in [Ohm]
    resistance_shunt : numeric
        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]
    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]
    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].
    breakdown_factor : float, 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 : float, default -5.5
        reverse breakdown voltage of the photovoltaic junction :math:`V_{br}`
        [V]
    breakdown_exp : float, default 3.28
        avalanche breakdown exponent :math:`m` [unitless]
    method : str, default 'newton'
       Either ``'newton'`` or ``'brentq'``. ''method'' must be ``'newton'``
       if ``breakdown_factor`` is not 0.
    method_kwargs : dict, optional
        Keyword arguments passed to root finder method. See
        :py:func:`scipy:scipy.optimize.brentq` and
        :py:func:`scipy:scipy.optimize.newton` parameters.
        ``'full_output': True`` is allowed, and ``optimizer_output`` would be
        returned. See examples section.

    Returns
    -------
    current : numeric
        current (I) at the specified voltage (V). [A]
    optimizer_output : tuple, optional, if specified in ``method_kwargs``
        see root finder documentation for selected method.
        Found root is diode voltage in [1]_.

    Examples
    --------
    Using the following arguments that may come from any
    `calcparams_.*` function in :py:mod:`pvlib.pvsystem`:

    >>> args = {'photocurrent': 1., 'saturation_current': 9e-10, 'nNsVth': 4.,
    ...         'resistance_series': 4., 'resistance_shunt': 5000.0}

    Use default values:

    >>> i = bishop88_i_from_v(0.0, **args)

    Specify tolerances and maximum number of iterations:

    >>> i = bishop88_i_from_v(0.0, **args, method='newton',
    ...     method_kwargs={'tol': 1e-3, 'rtol': 1e-3, 'maxiter': 20})

    Retrieve full output from the root finder:

    >>> i, method_output = bishop88_i_from_v(0.0, **args, method='newton',
    ...     method_kwargs={'full_output': True})

    References
    ----------
    .. [1] "Computer simulation of the effects of electrical mismatches in
       photovoltaic cell interconnection circuits" JW Bishop, Solar Cell (1988)
       :doi:`10.1016/0379-6787(88)90059-2`
    """
    # collect args
    args = (photocurrent, saturation_current,
            resistance_series, resistance_shunt, nNsVth, d2mutau, NsVbi,
            breakdown_factor, breakdown_voltage, breakdown_exp)
    method = method.lower()

    # method_kwargs create dict if not provided
    # this pattern avoids bugs with Mutable Default Parameters
    if not method_kwargs:
        method_kwargs = {}

    def fv(x, v, *a):
        # calculate voltage residual given diode voltage "x"
        return bishop88(x, *a)[1] - v

    if method == 'brentq':
        # first bound the search using voc
        voc_est = estimate_voc(photocurrent, saturation_current, nNsVth)
        # start iteration slightly less than NsVbi when voc_est > NsVbi, to
        # avoid the asymptote at NsVbi
        xp = np.where(voc_est < NsVbi, voc_est, 0.9999*NsVbi)

        # 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,
                          breakdown_factor, breakdown_voltage, breakdown_exp):
            return brentq(fv, 0.0, voc,
                          args=(v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
                                breakdown_factor, breakdown_voltage,
                                breakdown_exp),
                          **method_kwargs)

        vd_from_brent_vectorized = np.vectorize(vd_from_brent)
        vd = vd_from_brent_vectorized(xp, voltage, *args)
    elif method == 'newton':
        x0, (voltage, *args), method_kwargs = \
            _prepare_newton_inputs(voltage, (voltage, *args), method_kwargs)
        vd = newton(func=lambda x, *a: fv(x, voltage, *a), x0=x0,
                    fprime=lambda x, *a: bishop88(x, *a, gradients=True)[4],
                    args=args, **method_kwargs)
    else:
        raise NotImplementedError("Method '%s' isn't implemented" % method)

    # When 'full_output' parameter is specified, returned 'vd' is a tuple with
    # many elements, where the root is the first one. So we use it to output
    # the bishop88 result and return tuple(scalar, tuple with method results)
    if method_kwargs.get('full_output') is True:
        return (bishop88(vd[0], *args)[0], vd)
    else:
        return bishop88(vd, *args)[0]


def bishop88_v_from_i(current, photocurrent, saturation_current,
                      resistance_series, resistance_shunt, nNsVth,
                      d2mutau=0, NsVbi=np.inf, breakdown_factor=0.,
                      breakdown_voltage=-5.5, breakdown_exp=3.28,
                      method='newton', method_kwargs=None):
    """
    Find voltage given any current.

    Parameters
    ----------
    current : numeric
        current (I) in amperes [A]
    photocurrent : numeric
        photogenerated current (Iph or IL) [A]
    saturation_current : numeric
        diode dark or saturation current (Io or Isat) [A]
    resistance_series : numeric
        series resistance (Rs) in [Ohm]
    resistance_shunt : numeric
        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]
    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]
    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].
    breakdown_factor : float, 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 : float, default -5.5
        reverse breakdown voltage of the photovoltaic junction :math:`V_{br}`
        [V]
    breakdown_exp : float, default 3.28
        avalanche breakdown exponent :math:`m` [unitless]
    method : str, default 'newton'
       Either ``'newton'`` or ``'brentq'``. ''method'' must be ``'newton'``
       if ``breakdown_factor`` is not 0.
    method_kwargs : dict, optional
        Keyword arguments passed to root finder method. See
        :py:func:`scipy:scipy.optimize.brentq` and
        :py:func:`scipy:scipy.optimize.newton` parameters.
        ``'full_output': True`` is allowed, and ``optimizer_output`` would be
        returned. See examples section.

    Returns
    -------
    voltage : numeric
        voltage (V) at the specified current (I) in volts [V]
    optimizer_output : tuple, optional, if specified in ``method_kwargs``
        see root finder documentation for selected method.
        Found root is diode voltage in [1]_.

    Examples
    --------
    Using the following arguments that may come from any
    `calcparams_.*` function in :py:mod:`pvlib.pvsystem`:

    >>> args = {'photocurrent': 1., 'saturation_current': 9e-10, 'nNsVth': 4.,
    ...         'resistance_series': 4., 'resistance_shunt': 5000.0}

    Use default values:

    >>> v = bishop88_v_from_i(0.0, **args)

    Specify tolerances and maximum number of iterations:

    >>> v = bishop88_v_from_i(0.0, **args, method='newton',
    ...     method_kwargs={'tol': 1e-3, 'rtol': 1e-3, 'maxiter': 20})

    Retrieve full output from the root finder:

    >>> v, method_output = bishop88_v_from_i(0.0, **args, method='newton',
    ...     method_kwargs={'full_output': True})

    References
    ----------
    .. [1] "Computer simulation of the effects of electrical mismatches in
       photovoltaic cell interconnection circuits" JW Bishop, Solar Cell (1988)
       :doi:`10.1016/0379-6787(88)90059-2`
    """
    # collect args
    args = (photocurrent, saturation_current,
            resistance_series, resistance_shunt, nNsVth, d2mutau, NsVbi,
            breakdown_factor, breakdown_voltage, breakdown_exp)
    method = method.lower()

    # method_kwargs create dict if not provided
    # this pattern avoids bugs with Mutable Default Parameters
    if not method_kwargs:
        method_kwargs = {}

    # first bound the search using voc
    voc_est = estimate_voc(photocurrent, saturation_current, nNsVth)
    # start iteration slightly less than NsVbi when voc_est > NsVbi, to avoid
    # the asymptote at NsVbi
    xp = np.where(voc_est < NsVbi, voc_est, 0.9999*NsVbi)

    def fi(x, i, *a):
        # calculate current residual given diode voltage "x"
        return bishop88(x, *a)[0] - i

    if method == 'brentq':
        # 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,
                          breakdown_factor, breakdown_voltage, breakdown_exp):
            return brentq(fi, 0.0, voc,
                          args=(i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
                                breakdown_factor, breakdown_voltage,
                                breakdown_exp),
                          **method_kwargs)

        vd_from_brent_vectorized = np.vectorize(vd_from_brent)
        vd = vd_from_brent_vectorized(xp, current, *args)
    elif method == 'newton':
        x0, (current, *args), method_kwargs = \
            _prepare_newton_inputs(xp, (current, *args), method_kwargs)
        vd = newton(func=lambda x, *a: fi(x, current, *a), x0=x0,
                    fprime=lambda x, *a: bishop88(x, *a, gradients=True)[3],
                    args=args, **method_kwargs)
    else:
        raise NotImplementedError("Method '%s' isn't implemented" % method)

    # When 'full_output' parameter is specified, returned 'vd' is a tuple with
    # many elements, where the root is the first one. So we use it to output
    # the bishop88 result and return tuple(scalar, tuple with method results)
    if method_kwargs.get('full_output') is True:
        return (bishop88(vd[0], *args)[1], vd)
    else:
        return bishop88(vd, *args)[1]


def bishop88_mpp(photocurrent, saturation_current, resistance_series,
                 resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.inf,
                 breakdown_factor=0., breakdown_voltage=-5.5,
                 breakdown_exp=3.28, method='newton', method_kwargs=None):
    """
    Find max power point.

    Parameters
    ----------
    photocurrent : numeric
        photogenerated current (Iph or IL) [A]
    saturation_current : numeric
        diode dark or saturation current (Io or Isat) [A]
    resistance_series : numeric
        series resistance (Rs) in [Ohm]
    resistance_shunt : numeric
        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]
    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]
    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].
    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]
    breakdown_exp : numeric, default 3.28
        avalanche breakdown exponent :math:`m` [unitless]
    method : str, default 'newton'
       Either ``'newton'`` or ``'brentq'``. ''method'' must be ``'newton'``
       if ``breakdown_factor`` is not 0.
    method_kwargs : dict, optional
        Keyword arguments passed to root finder method. See
        :py:func:`scipy:scipy.optimize.brentq` and
        :py:func:`scipy:scipy.optimize.newton` parameters.
        ``'full_output': True`` is allowed, and ``optimizer_output`` would be
        returned. See examples section.

    Returns
    -------
    tuple
        max power current ``i_mp`` [A], max power voltage ``v_mp`` [V], and
        max power ``p_mp`` [W]
    optimizer_output : tuple, optional, if specified in ``method_kwargs``
        see root finder documentation for selected method.
        Found root is diode voltage in [1]_.

    Examples
    --------
    Using the following arguments that may come from any
    `calcparams_.*` function in :py:mod:`pvlib.pvsystem`:

    >>> args = {'photocurrent': 1., 'saturation_current': 9e-10, 'nNsVth': 4.,
    ...         'resistance_series': 4., 'resistance_shunt': 5000.0}

    Use default values:

    >>> i_mp, v_mp, p_mp = bishop88_mpp(**args)

    Specify tolerances and maximum number of iterations:

    >>> i_mp, v_mp, p_mp = bishop88_mpp(**args, method='newton',
    ...     method_kwargs={'tol': 1e-3, 'rtol': 1e-3, 'maxiter': 20})

    Retrieve full output from the root finder:

    >>> (i_mp, v_mp, p_mp), method_output = bishop88_mpp(**args,
    ...     method='newton', method_kwargs={'full_output': True})

    References
    ----------
    .. [1] "Computer simulation of the effects of electrical mismatches in
       photovoltaic cell interconnection circuits" JW Bishop, Solar Cell (1988)
       :doi:`10.1016/0379-6787(88)90059-2`
    """
    # collect args
    args = (photocurrent, saturation_current,
            resistance_series, resistance_shunt, nNsVth, d2mutau, NsVbi,
            breakdown_factor, breakdown_voltage, breakdown_exp)
    method = method.lower()

    # method_kwargs create dict if not provided
    # this pattern avoids bugs with Mutable Default Parameters
    if not method_kwargs:
        method_kwargs = {}

    # first bound the search using voc
    voc_est = estimate_voc(photocurrent, saturation_current, nNsVth)
    # start iteration slightly less than NsVbi when voc_est > NsVbi, to avoid
    # the asymptote at NsVbi
    xp = np.where(voc_est < NsVbi, voc_est, 0.9999*NsVbi)

    def fmpp(x, *a):
        return bishop88(x, *a, gradients=True)[6]

    if method == '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),
                            **method_kwargs)
        )
        vd = vec_fun(xp, *args)
    elif method == 'newton':
        # make sure all args are numpy arrays if max size > 1
        # if voc_est is an array, then make a copy to use for initial guess, v0

        x0, args, method_kwargs = \
            _prepare_newton_inputs(xp, args, method_kwargs)
        vd = newton(func=fmpp, x0=x0,
                    fprime=lambda x, *a: bishop88(x, *a, gradients=True)[7],
                    args=args, **method_kwargs)
    else:
        raise NotImplementedError("Method '%s' isn't implemented" % method)

    # When 'full_output' parameter is specified, returned 'vd' is a tuple with
    # many elements, where the root is the first one. So we use it to output
    # the bishop88 result and return
    # tuple(tuple with bishop88 solution, tuple with method results)
    if method_kwargs.get('full_output') is True:
        return (bishop88(vd[0], *args), vd)
    else:
        return bishop88(vd, *args)


def _shape_of_max_size(*args):
    return max(((np.size(a), np.shape(a)) for a in args),
               key=lambda t: t[0])[1]


def _prepare_newton_inputs(x0, args, method_kwargs):
    """
    Make inputs compatible with Scipy's newton by:
    - converting all arguments (`x0` and `args`) into numpy.ndarrays if any
      argument is not a scalar.
    - broadcasting the initial guess `x0` to the shape of the argument with
      the greatest size.

    Parameters
    ----------
    x0: numeric
        Initial guess for newton.
    args: Iterable(numeric)
        Iterable of additional arguments to use in SciPy's newton.
    method_kwargs: dict
        Options to pass to newton.

    Returns
    -------
    tuple
        The updated initial guess, arguments, and options for newton.
    """
    if not (np.isscalar(x0) and all(map(np.isscalar, args))):
        args = tuple(map(np.asarray, args))
        x0 = np.broadcast_to(x0, _shape_of_max_size(x0, *args))

    # set abs tolerance and maxiter from method_kwargs if not provided
    # apply defaults, but giving priority to user-specified values
    method_kwargs = {**NEWTON_DEFAULT_PARAMS, **method_kwargs}

    return x0, args, method_kwargs


def _lambertw_v_from_i(current, photocurrent, saturation_current,
                       resistance_series, resistance_shunt, nNsVth):
    # Record if inputs were all scalar
    output_is_scalar = all(map(np.isscalar,
                               (current, photocurrent, saturation_current,
                                resistance_series, resistance_shunt, nNsVth)))

    # This transforms Gsh=1/Rsh, including ideal Rsh=np.inf into Gsh=0., which
    #  is generally more numerically stable
    conductance_shunt = 1. / resistance_shunt

    # Ensure that we are working with read-only views of numpy arrays
    # Turns Series into arrays so that we don't have to worry about
    #  multidimensional broadcasting failing
    I, IL, I0, Rs, Gsh, a = \
        np.broadcast_arrays(current, photocurrent, saturation_current,
                            resistance_series, conductance_shunt, nNsVth)

    # Intitalize output V (I might not be float64)
    V = np.full_like(I, np.nan, dtype=np.float64)

    # Determine indices where 0 < Gsh requires implicit model solution
    idx_p = 0. < Gsh

    # Determine indices where 0 = Gsh allows explicit model solution
    idx_z = 0. == Gsh

    # Explicit solutions where Gsh=0
    if np.any(idx_z):
        V[idx_z] = a[idx_z] * np.log1p((IL[idx_z] - I[idx_z]) / I0[idx_z]) - \
                   I[idx_z] * Rs[idx_z]

    # Only compute using LambertW if there are cases with Gsh>0
    if np.any(idx_p):

        # use only the relevant subset for what follows
        I = I[idx_p]
        IL = IL[idx_p]
        I0 = I0[idx_p]
        Rs = Rs[idx_p]
        Gsh = Gsh[idx_p]
        a = a[idx_p]

        # LambertW argument, cannot be float128, may overflow to np.inf
        # overflow is explicitly handled below, so ignore warnings here
        with np.errstate(over='ignore'):
            argW = I0 / (Gsh * a) * np.exp((-I + IL + I0) / (Gsh * a))

        # lambertw typically returns complex value with zero imaginary part
        # may overflow to np.inf
        lambertwterm = lambertw(argW).real

        # Record indices where lambertw input overflowed output
        idx_inf = np.isinf(lambertwterm)

        # Only re-compute LambertW if it overflowed
        if np.any(idx_inf):
            # Calculate using log(argW) in case argW is really big
            logargW = (np.log(I0[idx_inf]) - np.log(Gsh[idx_inf]) -
                       np.log(a[idx_inf]) +
                       (-I[idx_inf] + IL[idx_inf] + I0[idx_inf]) /
                       (Gsh[idx_inf] * a[idx_inf]))

            # Three iterations of Newton-Raphson method to solve
            #  w+log(w)=logargW. The initial guess is w=logargW. Where direct
            #  evaluation (above) results in NaN from overflow, 3 iterations
            #  of Newton's method gives approximately 8 digits of precision.
            w = logargW
            for _ in range(0, 3):
                w = w * (1. - np.log(w) + logargW) / (1. + w)
            lambertwterm[idx_inf] = w

        # Eqn. 3 in Jain and Kapoor, 2004
        #  V = -I*(Rs + Rsh) + IL*Rsh - a*lambertwterm + I0*Rsh
        # Recast in terms of Gsh=1/Rsh for better numerical stability.
        V[idx_p] = (IL + I0 - I) / Gsh - I * Rs - a * lambertwterm

    if output_is_scalar:
        return V.item()
    else:
        return V


def _lambertw_i_from_v(voltage, photocurrent, saturation_current,
                       resistance_series, resistance_shunt, nNsVth):
    # Record if inputs were all scalar
    output_is_scalar = all(map(np.isscalar,
                               (voltage, photocurrent, saturation_current,
                                resistance_series, resistance_shunt, nNsVth)))

    # This transforms Gsh=1/Rsh, including ideal Rsh=np.inf into Gsh=0., which
    #  is generally more numerically stable
    conductance_shunt = 1. / resistance_shunt

    # Ensure that we are working with read-only views of numpy arrays
    # Turns Series into arrays so that we don't have to worry about
    #  multidimensional broadcasting failing
    V, IL, I0, Rs, Gsh, a = \
        np.broadcast_arrays(voltage, photocurrent, saturation_current,
                            resistance_series, conductance_shunt, nNsVth)

    # Intitalize output I (V might not be float64)
    I = np.full_like(V, np.nan, dtype=np.float64)           # noqa: E741, N806

    # Determine indices where 0 < Rs requires implicit model solution
    idx_p = 0. < Rs

    # Determine indices where 0 = Rs allows explicit model solution
    idx_z = 0. == Rs

    # Explicit solutions where Rs=0
    if np.any(idx_z):
        I[idx_z] = IL[idx_z] - I0[idx_z] * np.expm1(V[idx_z] / a[idx_z]) - \
                   Gsh[idx_z] * V[idx_z]

    # Only compute using LambertW if there are cases with Rs>0
    # Does NOT handle possibility of overflow, github issue 298
    if np.any(idx_p):
        # LambertW argument, cannot be float128, may overflow to np.inf
        argW = Rs[idx_p] * I0[idx_p] / (
                    a[idx_p] * (Rs[idx_p] * Gsh[idx_p] + 1.)) * \
               np.exp((Rs[idx_p] * (IL[idx_p] + I0[idx_p]) + V[idx_p]) /
                      (a[idx_p] * (Rs[idx_p] * Gsh[idx_p] + 1.)))

        # lambertw typically returns complex value with zero imaginary part
        # may overflow to np.inf
        lambertwterm = lambertw(argW).real

        # Eqn. 2 in Jain and Kapoor, 2004
        #  I = -V/(Rs + Rsh) - (a/Rs)*lambertwterm + Rsh*(IL + I0)/(Rs + Rsh)
        # Recast in terms of Gsh=1/Rsh for better numerical stability.
        I[idx_p] = (IL[idx_p] + I0[idx_p] - V[idx_p] * Gsh[idx_p]) / \
                   (Rs[idx_p] * Gsh[idx_p] + 1.) - (
                               a[idx_p] / Rs[idx_p]) * lambertwterm

    if output_is_scalar:
        return I.item()
    else:
        return I


def _lambertw(photocurrent, saturation_current, resistance_series,
              resistance_shunt, nNsVth, ivcurve_pnts=None):
    # collect args
    params = {'photocurrent': photocurrent,
              'saturation_current': saturation_current,
              'resistance_series': resistance_series,
              'resistance_shunt': resistance_shunt, 'nNsVth': nNsVth}

    # Compute short circuit current
    i_sc = _lambertw_i_from_v(0., **params)

    # Compute open circuit voltage
    v_oc = _lambertw_v_from_i(0., **params)

    # Set small elements <0 in v_oc to 0
    if isinstance(v_oc, np.ndarray):
        v_oc[(v_oc < 0) & (v_oc > -1e-12)] = 0.
    elif isinstance(v_oc, (float, int)):
        if v_oc < 0 and v_oc > -1e-12:
            v_oc = 0.

    # Find the voltage, v_mp, where the power is maximized.
    # Start the golden section search at v_oc * 1.14
    p_mp, v_mp = _golden_sect_DataFrame(params, 0., v_oc * 1.14, _pwr_optfcn)

    # Find Imp using Lambert W
    i_mp = _lambertw_i_from_v(v_mp, **params)

    # Find Ix and Ixx using Lambert W
    i_x = _lambertw_i_from_v(0.5 * v_oc, **params)

    i_xx = _lambertw_i_from_v(0.5 * (v_oc + v_mp), **params)

    out = (i_sc, v_oc, i_mp, v_mp, p_mp, i_x, i_xx)

    # create ivcurve
    if ivcurve_pnts:
        ivcurve_v = (np.asarray(v_oc)[..., np.newaxis] *
                     np.linspace(0, 1, ivcurve_pnts))

        ivcurve_i = _lambertw_i_from_v(ivcurve_v.T, **params).T

        out += (ivcurve_i, ivcurve_v)

    return out


def _pwr_optfcn(df, loc):
    '''
    Function to find power from ``i_from_v``.
    '''

    current = _lambertw_i_from_v(df[loc], df['photocurrent'],
                                 df['saturation_current'],
                                 df['resistance_series'],
                                 df['resistance_shunt'], df['nNsVth'])

    return current * df[loc]


def batzelis_keypoints(photocurrent, saturation_current, resistance_series,
                       resistance_shunt, nNsVth):
    """
    Estimate maximum power, open-circuit, and short-circuit points from
    single-diode equation parameters using Batzelis's method.

    This method is described in Section II.B of [1]_.

    Parameters
    ----------
    photocurrent : numeric
        Light-generated current. [A]
    saturation_current : numeric
        Diode saturation current. [A]
    resistance_series : numeric
        Series resistance. [Ohm]
    resistance_shunt : numeric
        Shunt resistance. [Ohm]
    nNsVth : numeric
        The product of the usual diode ideality factor (n, unitless),
        number of cells in series (Ns), and cell thermal voltage at
        specified effective irradiance and cell temperature. [V]

    Returns
    -------
    dict
        The returned dict-like object contains the keys/columns:

        * ``p_mp`` - power at maximum power point. [W]
        * ``i_mp`` - current at maximum power point. [A]
        * ``v_mp`` - voltage at maximum power point. [V]
        * ``i_sc`` - short circuit current. [A]
        * ``v_oc`` - open circuit voltage. [V]

    References
    ----------
    .. [1] E. I. Batzelis, "Simple PV Performance Equations Theoretically Well
       Founded on the Single-Diode Model," Journal of Photovoltaics vol. 7,
       no. 5, pp. 1400-1409, Sep 2017, :doi:`10.1109/JPHOTOV.2017.2711431`
    """
    # convenience variables
    Iph = photocurrent
    Is = saturation_current
    Rsh = resistance_shunt
    Rs = resistance_series
    a = nNsVth

    # Eqs 3-4
    isc = Iph * Rsh / (Rs + Rsh)
    with np.errstate(divide='ignore'):  # zero Iph
        voc = a * np.log(Iph / Is)

    # Eqs 5-8
    w = np.real(lambertw(np.e * Iph / Is))
    # vmp = (1 + Rs/Rsh) * a * (w - 1) - Rs * Iph * (1 - 1/w)  # not needed
    with np.errstate(divide='ignore', invalid='ignore'):  # zero Iph -> zero w
        imp = Iph * (1 - 1/w) - a * (w - 1) / Rsh

    vmp = a * (w - 1) - Rs * imp

    vmp = np.where(Iph > 0, vmp, 0)
    voc = np.where(Iph > 0, voc, 0)
    imp = np.where(Iph > 0, imp, 0)
    isc = np.where(Iph > 0, isc, 0)

    out = {
        'p_mp': imp * vmp,
        'i_mp': imp,
        'v_mp': vmp,
        'i_sc': isc,
        'v_oc': voc,
    }

    # if pandas in, ensure pandas out
    pandas_inputs = [
        x for x in [Iph, Is, Rsh, Rs, a] if isinstance(x, pd.Series)]
    if pandas_inputs:
        out = pd.DataFrame(out, index=pandas_inputs[0].index)

    return out