From 9ab2e951f1ab4bbad4565c594422b29efdb559b7 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Wed, 14 Jun 2017 22:28:22 -0600
Subject: [PATCH 01/22] Add draft i_from_v_alt function with tests

---
 pvlib/pvsystem.py           | 115 ++++++++++++++++++++++++++++++++---
 pvlib/test/test_pvsystem.py | 118 +++++++++++++++++++++++++++++++++++-
 2 files changed, 224 insertions(+), 9 deletions(-)

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 823fd8b28e..0105f3e6e9 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1913,7 +1913,7 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
 
     Returns
     -------
-    current : np.array
+    current : np.ndarray or np.float64
 
     References
     ----------
@@ -1928,12 +1928,14 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
 
     # asarray turns Series into arrays so that we don't have to worry
     # about multidimensional broadcasting failing
-    Rsh = np.asarray(resistance_shunt)
-    Rs = np.asarray(resistance_series)
-    I0 = np.asarray(saturation_current)
-    IL = np.asarray(photocurrent)
-    V = np.asarray(voltage)
-
+    Rsh = np.asarray(resistance_shunt, np.float64)
+    Rs = np.asarray(resistance_series, np.float64)
+    nNsVth = np.asarray(nNsVth, np.float64)
+    V = np.asarray(voltage, np.float64)
+    I0 = np.asarray(saturation_current, np.float64)
+    IL = np.asarray(photocurrent, np.float64)
+    
+    # argW cannot be float128
     argW = (Rs*I0*Rsh *
             np.exp(Rsh*(Rs*(IL+I0)+V) / (nNsVth*(Rs+Rsh))) /
             (nNsVth*(Rs + Rsh)))
@@ -1941,10 +1943,107 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
 
     # Eqn. 4 in Jain and Kapoor, 2004
     I = -V/(Rs + Rsh) - (nNsVth/Rs)*lambertwterm + Rsh*(IL + I0)/(Rs + Rsh)
-
+    
     return I
 
 
+def current_sum_at_diode_node(V=None, I=None, IL=None, I0=None, nNsVth=None, Rs=None, Rsh=None):
+    '''
+    Parameters
+    ----------
+    V : numpy.ndarray or numpy.float64 or python scalar
+        Device voltage [V]
+    I : numpy.ndarray or numpy.float64 or python scalar
+        Device current [A]
+    IL : numpy.ndarray or numpy.float64 or python scalar
+        Device photocurrent [A]
+    I0 : numpy.ndarray or numpy.float64 or python scalar
+        Device saturation current [A]
+    nNsVth : numpy.ndarray or numpy.float64 or python scalar
+        Device thermal voltage [V]
+    Rs : numpy.ndarray or numpy.float64 or python scalar
+        Device series resistance [Ohm]
+    Rsh : numpy.ndarray or numpy.float64 or python scalar
+        Device shunt resistance [Ohm]
+    
+    Returns
+    -------
+    current_sum_at_diode_node : numpy.ndarray or numpy.float64
+        Sum of currents at the diode node in equivalent circuit model [A]
+    '''
+    
+    # TODO Ensure that this qualifies as a numpy ufunc
+    
+    return IL - I0*np.expm1((V + I*Rs)/nNsVth) - (V + I*Rs)/Rsh - I # current_sum_at_diode_node
+
+
+def i_from_v_alt(V=None, IL=None, I0=None, nNsVth=None, Rs=None, Rsh=None, return_meta_dict=False):
+    '''
+    Parameters
+    ----------
+    V : numpy.ndarray or numpy.float64 or python scalar
+        Device voltage [V]
+    IL : numpy.ndarray or numpy.float64 or python scalar
+        Device photocurrent [A]
+    I0 : numpy.ndarray or numpy.float64 or python scalar
+        Device saturation current [A]
+    nNsVth : numpy.ndarray or numpy.float64 or python scalar
+        Device thermal voltage [V]
+    Rs : numpy.ndarray or numpy.float64 or python scalar
+        Device series resistance [Ohm]
+    Rsh : numpy.ndarray or numpy.float64 or python scalar
+        Device shunt resistance [Ohm]
+    return_meta_dict : boolean
+        Return additional computation metadata dictionary
+    
+    Returns
+    -------
+    I : numpy.ndarray or numpy.float64
+        Device current [A]
+    meta_dict : python dictionary (optional, returned when return_meta_dict=True)
+        Metadata for computation
+            meta_dict['current_sum_at_diode_node'] : numpy.ndarray or numpy.float64 like I
+                Sum of currents at the diode node in equivalent circuit model [A]
+            meta_dict['zero_Rs_idx'] : boolean like I
+                Indices where zero series resistance gives best solution
+    '''
+
+    # TODO Check if this qualifies as a numpy ufunc
+        
+    # Work with np.ndarray inputs
+    V = np.asarray(V)
+    IL = np.asarray(IL)
+    I0 = np.asarray(I0)
+    nNsVth = np.asarray(nNsVth)
+    Rs = np.asarray(Rs)
+    Rsh = np.asarray(Rsh)
+    
+    # Default computation of I using Rs=np.full_like(Rs, 0.), does not lose Rs shape or type info
+    I_times_Rs_zero = np.full_like(Rs, 0.)
+    I = IL - I0*np.expm1((V + I_times_Rs_zero)/nNsVth) - (V + I_times_Rs_zero)/Rsh
+    current_sum_at_diode_node_out = current_sum_at_diode_node(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    
+    # Computation of I using LambertW for provided Rs
+    I_LambertW = i_from_v(Rsh, Rs, nNsVth, V, I0, IL)
+    current_sum_at_diode_node_LambertW = current_sum_at_diode_node(V=V, I=I_LambertW, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    
+    # Compute selection indices (may be a scalar boolean)
+    zero_Rs_idx = np.logical_and(np.isfinite(current_sum_at_diode_node_LambertW), np.absolute(current_sum_at_diode_node_LambertW) <= np.absolute(current_sum_at_diode_node_out))
+    
+    if np.isscalar(I):
+        if zero_Rs_idx:
+            I = I_LambertW
+            current_sum_at_diode_node_out = current_sum_at_diode_node_LambertW
+    else:
+        I[zero_Rs_idx] = I_LambertW[zero_Rs_idx]
+        current_sum_at_diode_node_out[zero_Rs_idx] = current_sum_at_diode_node_LambertW[zero_Rs_idx]
+    
+    if return_meta_dict:
+        return I, {'current_sum_at_diode_node' : current_sum_at_diode_node_out, 'zero_Rs_idx' : zero_Rs_idx}
+    else:
+        return I
+
+
 def snlinverter(v_dc, p_dc, inverter):
     r'''
     Converts DC power and voltage to AC power using Sandia's
diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index b2d22084cd..7740e2b386 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -9,7 +9,7 @@
 
 import pytest
 from pandas.util.testing import assert_series_equal, assert_frame_equal
-from numpy.testing import assert_allclose
+from numpy.testing import assert_allclose, assert_array_equal
 
 from pvlib import tmy
 from pvlib import pvsystem
@@ -353,6 +353,122 @@ def test_PVSystem_calcparams_desoto(cec_module_params):
     assert_allclose(nNsVth, 0.473)
 
 
+def test_current_sum_at_diode_node():
+    V = np.array([40., 0. , 0])
+    I = np.array([0., 0., 3.])
+    IL = 7.
+    I0 = 6.e-7
+    nNsVth = 0.5 
+    Rs = 0.1
+    Rsh = 20.
+    
+    results_1 = np.full_like(V, np.nan)
+    results_2 = np.full_like(V, np.nan)
+    
+    results_1[0] = pvsystem.current_sum_at_diode_node(V=V[0], I=I[0], IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    results_2[0] = IL - I0*np.expm1(V[0]/nNsVth) - V[0]/Rsh
+    
+    results_1[1] = pvsystem.current_sum_at_diode_node(V=V[1], I=I[1], IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    results_2[1] = IL
+    
+    results_1[2] = pvsystem.current_sum_at_diode_node(V=V[2], I=I[2], IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    results_2[2] = IL - I0*np.expm1(I[2]*Rs/nNsVth) - I[2]*Rs/Rsh - I[2]
+    
+    assert_array_equal(results_1, results_2)
+    
+    results_vec = pvsystem.current_sum_at_diode_node(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    assert_array_equal(results_vec, results_2)
+
+    V = 0.
+    I = 0.
+    nNsVth = np.asarray([0.45, 0.5, 0.55])
+    results_vec = pvsystem.current_sum_at_diode_node(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    assert_array_equal(results_vec, np.asarray([IL, IL, IL]))
+
+    nNsVth = pd.Series([0.45, 0.5, 0.55])
+    results_series = pvsystem.current_sum_at_diode_node(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    assert_series_equal(results_series, pd.Series([IL, IL, IL]))    
+
+    V = 40.
+    I = 3.
+    nNsVth = 0.5
+    Rsh = np.inf
+    results_inf_Rsh_1 = pvsystem.current_sum_at_diode_node(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    results_inf_Rsh_2 = IL - I0*np.expm1((V + I*Rs)/nNsVth) - I
+    assert_array_equal(results_inf_Rsh_1, results_inf_Rsh_2)
+
+
+@requires_scipy
+def test_i_from_v_alt():
+    # Solution set of Python scalars
+    Rsh = 20.
+    Rs = 0.1
+    nNsVth = 0.5
+    V = 40.
+    I0 = 6.e-7
+    IL = 7.
+    I = -299.746389916
+    
+    # Convergence criteria
+    atol = 1.e-12
+    
+    # Can handle all python scalar inputs
+    I_out = pvsystem.i_from_v_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL)
+    I_expected = np.float64(I)
+    assert(isinstance(I_out, type(I_expected)))
+    assert_allclose(I_out, I_expected)
+    _, meta_dict = pvsystem.i_from_v_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL, return_meta_dict=True)
+    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
+
+    # Can handle all rank-0 array inputs
+    I_out = pvsystem.i_from_v_alt(Rsh=np.array(Rsh), Rs=np.array(Rs), nNsVth=np.array(nNsVth), V=np.array(V), I0=np.array(I0), IL=np.array(IL))
+    I_expected = np.float64(I)
+    assert(isinstance(I_out, type(I_expected)))
+    assert_allclose(I_out, I_expected)
+    _, meta_dict = pvsystem.i_from_v_alt(Rsh=np.array(Rsh), Rs=np.array(Rs), nNsVth=np.array(nNsVth), V=np.array(V), I0=np.array(I0), IL=np.array(IL), return_meta_dict=True)
+    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
+
+    # Can handle all rank-1 singleton array inputs
+    I_out = pvsystem.i_from_v_alt(Rsh=np.array([Rsh]), Rs=np.array([Rs]), nNsVth=np.array([nNsVth]), V=np.array([V]), I0=np.array([I0]), IL=np.array([IL]))
+    I_expected = np.array([I])
+    assert(isinstance(I_out, type(I_expected)))
+    assert_allclose(I_out, I_expected)
+    _, meta_dict = pvsystem.i_from_v_alt(Rsh=np.array([Rsh]), Rs=np.array([Rs]), nNsVth=np.array([nNsVth]), V=np.array([V]), I0=np.array([I0]), IL=np.array([IL]), return_meta_dict=True)
+    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
+
+    # Can handle all rank-1 non-singleton array inputs with a zero series resistance (Rs=0 gives I=IL=Isc at V=0)
+    I_out = pvsystem.i_from_v_alt(Rsh=np.array([Rsh, Rsh]), Rs=np.array([0., Rs]), nNsVth=np.array([nNsVth, nNsVth]), V=np.array([0., V]), I0=np.array([I0, I0]), IL=np.array([IL, IL]))
+    I_expected = np.array([IL, I])
+    assert(isinstance(I_out, type(I_expected)))
+    assert_allclose(I_out, I_expected)
+    _, meta_dict = pvsystem.i_from_v_alt(Rsh=np.array([Rsh, Rsh]), Rs=np.array([0., Rs]), nNsVth=np.array([nNsVth, nNsVth]), V=np.array([0., V]), I0=np.array([I0, I0]), IL=np.array([IL, IL]), return_meta_dict=True)    
+    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
+    
+    # Can handle mixed inputs with a rank-2 array with zero series resistance (Rs=0 gives I=IL=Isc at V=0)
+    I_out = pvsystem.i_from_v_alt(Rsh=np.array([Rsh]), Rs=np.array([[0., 0.], [0., 0.]]), nNsVth=np.array(nNsVth), V=np.array(0.), I0=np.array([I0]), IL=np.array([IL]))
+    I_expected = np.array([[IL, IL], [IL, IL]])
+    assert(isinstance(I_out, type(I_expected)))
+    assert_allclose(I_out, I_expected)
+    _, meta_dict = pvsystem.i_from_v_alt(Rsh=np.array([Rsh]), Rs=np.array([[0., 0.], [0., 0.]]), nNsVth=np.array(nNsVth), V=np.array(0.), I0=np.array([I0]), IL=np.array([IL]), return_meta_dict=True)
+    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
+    
+    # Can handle ideal series and shunt
+    V_oc = nNsVth*np.log(IL/I0 + 1.)
+    I_out = pvsystem.i_from_v_alt(Rsh=np.inf, Rs=0., nNsVth=nNsVth, V=np.array([0., V_oc/2., V_oc]), I0=I0, IL=IL)
+    I_expected = np.array([IL, IL - I0*np.expm1(V_oc/2./nNsVth), IL - I0*np.expm1(V_oc/nNsVth)])
+    assert(isinstance(I_out, type(I_expected)))
+    assert_allclose(I_out, I_expected)
+    _, meta_dict = pvsystem.i_from_v_alt(Rsh=np.inf, Rs=0., nNsVth=nNsVth, V=np.array([0., V_oc/2., V_oc]), I0=I0, IL=IL, return_meta_dict=True)
+    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
+    
+    # Can handle only ideal shunt
+    I_out = pvsystem.i_from_v_alt(Rsh=np.inf, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL)
+    assert(isinstance(I_out, np.float64))
+    assert(np.isfinite(I_out)) # No exact expression to evaluate
+    _, meta_dict = pvsystem.i_from_v_alt(Rsh=np.inf, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL, return_meta_dict=True)
+    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
+
+
 @requires_scipy
 def test_v_from_i():
     output = pvsystem.v_from_i(20, .1, .5, 3, 6e-7, 7)

From a19447b3f36f609d658307eeda9d208062c9cb74 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Wed, 14 Jun 2017 22:43:44 -0600
Subject: [PATCH 02/22] Better comments and more explicit typing

---
 pvlib/pvsystem.py           | 32 ++++++++++++++++++--------------
 pvlib/test/test_pvsystem.py |  1 +
 2 files changed, 19 insertions(+), 14 deletions(-)

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index e1b87e1758..f266b11a9d 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1830,7 +1830,7 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
 
     Returns
     -------
-    current : np.array
+    current : np.ndarray or np.float64
 
     References
     ----------
@@ -1843,12 +1843,16 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
     except ImportError:
         raise ImportError('This function requires scipy')
 
-    Rsh = resistance_shunt
-    Rs = resistance_series
-    I0 = saturation_current
-    IL = photocurrent
-    I = current
-
+    # asarray turns Series into arrays so that we don't have to worry
+    # about multidimensional broadcasting failing
+    Rsh = np.asarray(resistance_shunt, np.float64)
+    Rs = np.asarray(resistance_series, np.float64)
+    nNsVth = np.asarray(nNsVth, np.float64)
+    I = np.asarray(current, np.float64)
+    I0 = np.asarray(saturation_current, np.float64)
+    IL = np.asarray(photocurrent, np.float64)
+    
+    # argW cannot be float128
     argW = I0 * Rsh / nNsVth * np.exp(Rsh * (-I + IL + I0) / nNsVth)
     lambertwterm = lambertw(argW).real
 
@@ -2008,13 +2012,13 @@ def i_from_v_alt(V=None, IL=None, I0=None, nNsVth=None, Rs=None, Rsh=None, retur
 
     # TODO Check if this qualifies as a numpy ufunc
         
-    # Work with np.ndarray inputs
-    V = np.asarray(V)
-    IL = np.asarray(IL)
-    I0 = np.asarray(I0)
-    nNsVth = np.asarray(nNsVth)
-    Rs = np.asarray(Rs)
-    Rsh = np.asarray(Rsh)
+    # Ensure inputs are all np.ndarray with np.float64 type
+    V = np.asarray(V, np.float64)
+    IL = np.asarray(IL, np.float64)
+    I0 = np.asarray(I0, np.float64)
+    nNsVth = np.asarray(nNsVth, np.float64)
+    Rs = np.asarray(Rs, np.float64)
+    Rsh = np.asarray(Rsh, np.float64)
     
     # Default computation of I using Rs=np.full_like(Rs, 0.), does not lose Rs shape or type info
     I_times_Rs_zero = np.full_like(Rs, 0.)
diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index 69d4b46c3c..0a199d8092 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -461,6 +461,7 @@ def test_i_from_v_alt():
     _, meta_dict = pvsystem.i_from_v_alt(Rsh=np.inf, Rs=0., nNsVth=nNsVth, V=np.array([0., V_oc/2., V_oc]), I0=I0, IL=IL, return_meta_dict=True)
     assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
     
+    # THIS FAILS: WE SHOULD PROBABLY USE SHUNT CONDUCTANCE INSTEAD OF SHUNT RESISTANCE
     # Can handle only ideal shunt
     I_out = pvsystem.i_from_v_alt(Rsh=np.inf, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL)
     assert(isinstance(I_out, np.float64))

From 76616d51f18de8232e85b04e4954c5df692294d7 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Thu, 15 Jun 2017 13:32:01 -0600
Subject: [PATCH 03/22] Use transform from shunt resistance to shunt
 conductance

---
 pvlib/pvsystem.py           | 9 +++++----
 pvlib/test/test_pvsystem.py | 8 ++++----
 2 files changed, 9 insertions(+), 8 deletions(-)

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index f266b11a9d..73d82f7b16 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1937,14 +1937,15 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
     I0 = np.asarray(saturation_current, np.float64)
     IL = np.asarray(photocurrent, np.float64)
     
+    # This transforms any ideal Rsh=np.inf into Gsh=0., which is generally more numerically stable
+    Gsh = 1./Rsh
+    
     # argW cannot be float128
-    argW = (Rs*I0*Rsh *
-            np.exp(Rsh*(Rs*(IL+I0)+V) / (nNsVth*(Rs+Rsh))) /
-            (nNsVth*(Rs + Rsh)))
+    argW = Rs*I0/(nNsVth*(Rs*Gsh + 1.))*np.exp((Rs*(IL + I0) + V)/(nNsVth*(Rs*Gsh + 1.)))
     lambertwterm = lambertw(argW).real
 
     # Eqn. 4 in Jain and Kapoor, 2004
-    I = -V/(Rs + Rsh) - (nNsVth/Rs)*lambertwterm + Rsh*(IL + I0)/(Rs + Rsh)
+    I = (IL + I0 - V*Gsh)/(Rs*Gsh + 1.) - (nNsVth/Rs)*lambertwterm
     
     return I
 
diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index 0a199d8092..c034f2bb60 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -410,7 +410,7 @@ def test_i_from_v_alt():
     I = -299.746389916
     
     # Convergence criteria
-    atol = 1.e-12
+    atol = 1.e-11
     
     # Can handle all python scalar inputs
     I_out = pvsystem.i_from_v_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL)
@@ -588,10 +588,10 @@ def test_singlediode_series_ivcurve(cec_module_params):
                                          module_parameters=cec_module_params,
                                          EgRef=1.121,
                                          dEgdT=-0.0002677)
-
+    
     out = pvsystem.singlediode(IL, I0, Rs, Rsh, nNsVth, ivcurve_pnts=3)
-
-    expected = OrderedDict([('i_sc', array([        nan,  3.01054475,  6.00675648])),
+    
+    expected = OrderedDict([('i_sc', array([        0.,  3.01054475,  6.00675648])),
              ('v_oc', array([         nan,   9.96886962,  10.29530483])),
              ('i_mp', array([        nan,  2.65191983,  5.28594672])),
              ('v_mp', array([        nan,  8.33392491,  8.4159707 ])),

From e35412f4fa1ce55d11f7e250db297b33b5a274f3 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Thu, 15 Jun 2017 22:41:06 -0600
Subject: [PATCH 04/22] Add v_from_i_alt() with initial tests and use np.where

---
 pvlib/pvsystem.py           | 125 +++++++++++++++++++++++++++---------
 pvlib/test/test_pvsystem.py |  88 +++++++++++++++++++++++--
 2 files changed, 175 insertions(+), 38 deletions(-)

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 73d82f7b16..c130511cca 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1852,13 +1852,15 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
     I0 = np.asarray(saturation_current, np.float64)
     IL = np.asarray(photocurrent, np.float64)
     
+    # This transforms any ideal Rsh=np.inf into Gsh=0., which is generally more numerically stable
+    Gsh = 1./Rsh
+    
     # argW cannot be float128
-    argW = I0 * Rsh / nNsVth * np.exp(Rsh * (-I + IL + I0) / nNsVth)
+    argW = I0/(Gsh*nNsVth)*np.exp((-I + IL + I0)/(Gsh*nNsVth))
     lambertwterm = lambertw(argW).real
 
     # Calculate using log(argW) in case argW is really big
-    logargW = (np.log(I0) + np.log(Rsh) - np.log(nNsVth) +
-               Rsh * (-I + IL + I0) / nNsVth)
+    logargW = (np.log(I0) - np.log(Gsh) - np.log(nNsVth) + (-I + IL + I0)/(Gsh*nNsVth))
 
     # Three iterations of Newton-Raphson method to solve
     # w+log(w)=logargW. The initial guess is w=logargW. Where direct
@@ -1869,11 +1871,10 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
         w = w * (1 - np.log(w) + logargW) / (1 + w)
     lambertwterm_log = w
 
-    lambertwterm = np.where(np.isfinite(lambertwterm), lambertwterm,
-                            lambertwterm_log)
+    lambertwterm = np.where(np.isfinite(lambertwterm), lambertwterm, lambertwterm_log)
 
     # Eqn. 3 in Jain and Kapoor, 2004
-    V = -I*(Rs + Rsh) + IL*Rsh - nNsVth*lambertwterm + I0*Rsh
+    V = (IL + I0 - I)/Gsh - I*Rs - nNsVth*lambertwterm
 
     return V
 
@@ -1950,7 +1951,7 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
     return I
 
 
-def current_sum_at_diode_node(V=None, I=None, IL=None, I0=None, nNsVth=None, Rs=None, Rsh=None):
+def current_sum_at_diode_node(V, I, IL, I0, nNsVth, Rs, Rsh):
     '''
     Parameters
     ----------
@@ -1977,37 +1978,101 @@ def current_sum_at_diode_node(V=None, I=None, IL=None, I0=None, nNsVth=None, Rs=
     
     # TODO Ensure that this qualifies as a numpy ufunc
     
-    return IL - I0*np.expm1((V + I*Rs)/nNsVth) - (V + I*Rs)/Rsh - I # current_sum_at_diode_node
+    # current_sum_at_diode_node
+    return IL - I0*np.expm1((V + I*Rs)/nNsVth) - (V + I*Rs)/Rsh - I
 
 
-def i_from_v_alt(V=None, IL=None, I0=None, nNsVth=None, Rs=None, Rsh=None, return_meta_dict=False):
+def v_from_i_alt(I, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
     '''
     Parameters
     ----------
-    V : numpy.ndarray or numpy.float64 or python scalar
+    I : numpy.ndarray or scalar
+        Device current [A]
+    IL : numpy.ndarray or scalar
+        Device photocurrent [A]
+    I0 : numpy.ndarray or scalar
+        Device saturation current [A]
+    nNsVth : numpy.ndarray or scalar
+        Device thermal voltage [V]
+    Rs : numpy.ndarray or scalar
+        Device series resistance [Ohm]
+    Rsh : numpy.ndarray or scalar
+        Device shunt resistance [Ohm]
+    return_meta_dict : boolean
+        Return additional computation metadata dictionary
+    
+    Returns
+    -------
+    V : numpy.ndarray
         Device voltage [V]
-    IL : numpy.ndarray or numpy.float64 or python scalar
+    meta_dict : python dictionary (optional, returned when return_meta_dict=True)
+        Metadata for computation
+            meta_dict['current_sum_at_diode_node'] : numpy.ndarray like V
+                Sum of currents at the diode node in equivalent circuit model [A]
+            meta_dict['inf_Rsh_idx'] : boolean numpy.ndarray like V
+                Indices where infinite shunt resistance gives best solution
+    '''
+
+    # TODO Check if this qualifies as a numpy ufunc
+        
+    # Ensure inputs are all np.ndarray with np.float64 type
+    I = np.asarray(I, np.float64)
+    IL = np.asarray(IL, np.float64)
+    I0 = np.asarray(I0, np.float64)
+    nNsVth = np.asarray(nNsVth, np.float64)
+    Rs = np.asarray(Rs, np.float64)
+    Rsh = np.asarray(Rsh, np.float64)
+  
+    # Default computation of V using Rsh=np.full_like(Rsh, np.inf), does not lose Rsh shape or type info
+    zero_term = np.zeros_like(I*Rs/Rsh)
+    V = nNsVth*(np.log(IL - I - zero_term + I0) - np.log(I0)) - I*Rs
+    current_sum_at_diode_node_out = current_sum_at_diode_node(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    
+    # Computation of V using LambertW for provided Rsh
+    V_LambertW = v_from_i(Rsh, Rs, nNsVth, I, I0, IL)
+    current_sum_at_diode_node_LambertW = current_sum_at_diode_node(V=V_LambertW, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    
+    # Compute selection indices (may be a scalar boolean)
+    finite_Rsh_idx = np.logical_and(np.isfinite(current_sum_at_diode_node_LambertW), np.absolute(current_sum_at_diode_node_LambertW) <= np.absolute(current_sum_at_diode_node_out))
+    
+    # These are always np.ndarray
+    V = np.where(finite_Rsh_idx, V_LambertW, V)
+    current_sum_at_diode_node_out = np.where(finite_Rsh_idx, current_sum_at_diode_node_LambertW, current_sum_at_diode_node_out)
+        
+    if return_meta_dict:
+        return V, {'current_sum_at_diode_node' : current_sum_at_diode_node_out, 'inf_Rsh_idx' : np.logical_not(finite_Rsh_idx)}
+    else:
+        return V
+
+
+def i_from_v_alt(V, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
+    '''
+    Parameters
+    ----------
+    V : numpy.ndarray or scalar
+        Device voltage [V]
+    IL : numpy.ndarray or scalar
         Device photocurrent [A]
-    I0 : numpy.ndarray or numpy.float64 or python scalar
+    I0 : numpy.ndarray or scalar
         Device saturation current [A]
-    nNsVth : numpy.ndarray or numpy.float64 or python scalar
+    nNsVth : numpy.ndarray or scalar
         Device thermal voltage [V]
-    Rs : numpy.ndarray or numpy.float64 or python scalar
+    Rs : numpy.ndarray or scalar
         Device series resistance [Ohm]
-    Rsh : numpy.ndarray or numpy.float64 or python scalar
+    Rsh : numpy.ndarray or scalar
         Device shunt resistance [Ohm]
     return_meta_dict : boolean
         Return additional computation metadata dictionary
     
     Returns
     -------
-    I : numpy.ndarray or numpy.float64
+    I : numpy.ndarray
         Device current [A]
     meta_dict : python dictionary (optional, returned when return_meta_dict=True)
         Metadata for computation
-            meta_dict['current_sum_at_diode_node'] : numpy.ndarray or numpy.float64 like I
+            meta_dict['current_sum_at_diode_node'] : numpy.ndarray like I
                 Sum of currents at the diode node in equivalent circuit model [A]
-            meta_dict['zero_Rs_idx'] : boolean like I
+            meta_dict['zero_Rs_idx'] : boolean numpy.ndarray like I
                 Indices where zero series resistance gives best solution
     '''
 
@@ -2020,29 +2085,25 @@ def i_from_v_alt(V=None, IL=None, I0=None, nNsVth=None, Rs=None, Rsh=None, retur
     nNsVth = np.asarray(nNsVth, np.float64)
     Rs = np.asarray(Rs, np.float64)
     Rsh = np.asarray(Rsh, np.float64)
-    
+
     # Default computation of I using Rs=np.full_like(Rs, 0.), does not lose Rs shape or type info
-    I_times_Rs_zero = np.full_like(Rs, 0.)
-    I = IL - I0*np.expm1((V + I_times_Rs_zero)/nNsVth) - (V + I_times_Rs_zero)/Rsh
+    zero_term = np.zeros_like(Rs)
+    I = IL - I0*np.expm1((V + zero_term)/nNsVth) - (V + zero_term)/Rsh
     current_sum_at_diode_node_out = current_sum_at_diode_node(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
-    
+
     # Computation of I using LambertW for provided Rs
     I_LambertW = i_from_v(Rsh, Rs, nNsVth, V, I0, IL)
     current_sum_at_diode_node_LambertW = current_sum_at_diode_node(V=V, I=I_LambertW, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
-    
+
     # Compute selection indices (may be a scalar boolean)
-    zero_Rs_idx = np.logical_and(np.isfinite(current_sum_at_diode_node_LambertW), np.absolute(current_sum_at_diode_node_LambertW) <= np.absolute(current_sum_at_diode_node_out))
+    nonzero_Rs_idx = np.logical_and(np.isfinite(current_sum_at_diode_node_LambertW), np.absolute(current_sum_at_diode_node_LambertW) <= np.absolute(current_sum_at_diode_node_out))
     
-    if np.isscalar(I):
-        if zero_Rs_idx:
-            I = I_LambertW
-            current_sum_at_diode_node_out = current_sum_at_diode_node_LambertW
-    else:
-        I[zero_Rs_idx] = I_LambertW[zero_Rs_idx]
-        current_sum_at_diode_node_out[zero_Rs_idx] = current_sum_at_diode_node_LambertW[zero_Rs_idx]
+    # These are always np.ndarray
+    I = np.where(nonzero_Rs_idx, I_LambertW, I)
+    current_sum_at_diode_node_out = np.where(nonzero_Rs_idx, current_sum_at_diode_node_LambertW, current_sum_at_diode_node_out)
     
     if return_meta_dict:
-        return I, {'current_sum_at_diode_node' : current_sum_at_diode_node_out, 'zero_Rs_idx' : zero_Rs_idx}
+        return I, {'current_sum_at_diode_node' : current_sum_at_diode_node_out, 'zero_Rs_idx' : np.logical_not(nonzero_Rs_idx)}
     else:
         return I
 
diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index c034f2bb60..848c07ce22 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -414,7 +414,7 @@ def test_i_from_v_alt():
     
     # Can handle all python scalar inputs
     I_out = pvsystem.i_from_v_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL)
-    I_expected = np.float64(I)
+    I_expected = np.array(I)
     assert(isinstance(I_out, type(I_expected)))
     assert_allclose(I_out, I_expected)
     _, meta_dict = pvsystem.i_from_v_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL, return_meta_dict=True)
@@ -422,7 +422,7 @@ def test_i_from_v_alt():
 
     # Can handle all rank-0 array inputs
     I_out = pvsystem.i_from_v_alt(Rsh=np.array(Rsh), Rs=np.array(Rs), nNsVth=np.array(nNsVth), V=np.array(V), I0=np.array(I0), IL=np.array(IL))
-    I_expected = np.float64(I)
+    I_expected = np.array(I)
     assert(isinstance(I_out, type(I_expected)))
     assert_allclose(I_out, I_expected)
     _, meta_dict = pvsystem.i_from_v_alt(Rsh=np.array(Rsh), Rs=np.array(Rs), nNsVth=np.array(nNsVth), V=np.array(V), I0=np.array(I0), IL=np.array(IL), return_meta_dict=True)
@@ -453,7 +453,7 @@ def test_i_from_v_alt():
     assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
     
     # Can handle ideal series and shunt
-    V_oc = nNsVth*np.log(IL/I0 + 1.)
+    V_oc = nNsVth*(np.log(IL + I0) - np.log(I0))
     I_out = pvsystem.i_from_v_alt(Rsh=np.inf, Rs=0., nNsVth=nNsVth, V=np.array([0., V_oc/2., V_oc]), I0=I0, IL=IL)
     I_expected = np.array([IL, IL - I0*np.expm1(V_oc/2./nNsVth), IL - I0*np.expm1(V_oc/nNsVth)])
     assert(isinstance(I_out, type(I_expected)))
@@ -461,13 +461,89 @@ def test_i_from_v_alt():
     _, meta_dict = pvsystem.i_from_v_alt(Rsh=np.inf, Rs=0., nNsVth=nNsVth, V=np.array([0., V_oc/2., V_oc]), I0=I0, IL=IL, return_meta_dict=True)
     assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
     
-    # THIS FAILS: WE SHOULD PROBABLY USE SHUNT CONDUCTANCE INSTEAD OF SHUNT RESISTANCE
-    # Can handle only ideal shunt
+    # Can handle only ideal shunt resistance
     I_out = pvsystem.i_from_v_alt(Rsh=np.inf, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL)
-    assert(isinstance(I_out, np.float64))
+    assert(isinstance(I_out, np.ndarray))
     assert(np.isfinite(I_out)) # No exact expression to evaluate
     _, meta_dict = pvsystem.i_from_v_alt(Rsh=np.inf, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL, return_meta_dict=True)
     assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
+    
+    # TODO Stability as Rs->0^+ and/or Rsh->inf
+
+
+@requires_scipy
+def test_v_from_i_alt():
+    # Solution set of Python scalars
+    Rsh = 20.
+    Rs = 0.1
+    nNsVth = 0.5
+    I = 3.
+    I0 = 6.e-7
+    IL = 7.
+    V = 7.5049875193450521
+    
+    # Convergence criteria
+    atol = 1.e-11
+    
+    # Can handle all python scalar inputs
+    V_out = pvsystem.v_from_i_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, I=I, I0=I0, IL=IL)
+    V_expected = np.array(V)
+    assert(isinstance(V_out, type(V_expected)))
+    assert_allclose(V_out, V_expected)
+    _, meta_dict = pvsystem.v_from_i_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, I=I, I0=I0, IL=IL, return_meta_dict=True)
+    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
+
+    # Can handle all rank-0 array inputs
+    V_out = pvsystem.v_from_i_alt(Rsh=np.array(Rsh), Rs=np.array(Rs), nNsVth=np.array(nNsVth), I=np.array(I), I0=np.array(I0), IL=np.array(IL))
+    V_expected = np.array(V)
+    assert(isinstance(V_out, type(V_expected)))
+    assert_allclose(V_out, V_expected)
+    _, meta_dict = pvsystem.v_from_i_alt(Rsh=np.array(Rsh), Rs=np.array(Rs), nNsVth=np.array(nNsVth), I=np.array(I), I0=np.array(I0), IL=np.array(IL), return_meta_dict=True)
+    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
+
+    # Can handle all rank-1 singleton array inputs
+    V_out = pvsystem.v_from_i_alt(Rsh=np.array([Rsh]), Rs=np.array([Rs]), nNsVth=np.array([nNsVth]), I=np.array([I]), I0=np.array([I0]), IL=np.array([IL]))
+    V_expected = np.array([V])
+    assert(isinstance(V_out, type(V_expected)))
+    assert_allclose(V_out, V_expected)
+    _, meta_dict = pvsystem.v_from_i_alt(Rsh=np.array([Rsh]), Rs=np.array([Rs]), nNsVth=np.array([nNsVth]), I=np.array([I]), I0=np.array([I0]), IL=np.array([IL]), return_meta_dict=True)
+    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
+
+    # Can handle all rank-1 non-singleton array inputs with infinite shunt resistance (Rsh=inf gives V=Voc at I=0)
+    V_oc = nNsVth*(np.log(IL + I0) - np.log(I0))
+    V_out = pvsystem.v_from_i_alt(Rsh=np.array([np.inf, Rsh]), Rs=np.array([Rs, Rs]), nNsVth=np.array([nNsVth, nNsVth]), I=np.array([0., I]), I0=np.array([I0, I0]), IL=np.array([IL, IL]))
+    V_expected = np.array([V_oc, V])
+    assert(isinstance(V_out, type(V_expected)))
+    assert_allclose(V_out, V_expected)
+    _, meta_dict = pvsystem.v_from_i_alt(Rsh=np.array([np.inf, Rsh]), Rs=np.array([Rs, Rs]), nNsVth=np.array([nNsVth, nNsVth]), I=np.array([0., I]), I0=np.array([I0, I0]), IL=np.array([IL, IL]), return_meta_dict=True)    
+    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
+    
+    # Can handle mixed inputs with a rank-2 array with infinite shunt resistance (Rsh=inf gives V=Voc at I=0)
+    V_oc = nNsVth*(np.log(IL + I0) - np.log(I0))
+    V_out = pvsystem.v_from_i_alt(Rsh=np.array([[np.inf, np.inf], [np.inf, np.inf]]), Rs=np.array([Rs]), nNsVth=np.array(nNsVth), I=np.array(0.), I0=np.array([I0]), IL=np.array([IL]))
+    V_expected = np.array([[V_oc, V_oc], [V_oc, V_oc]])
+    assert(isinstance(V_out, type(V_expected)))
+    assert_allclose(V_out, V_expected)
+    _, meta_dict = pvsystem.v_from_i_alt(Rsh=np.array([[np.inf, np.inf], [np.inf, np.inf]]), Rs=np.array([Rs]), nNsVth=np.array(nNsVth), I=np.array(0.), I0=np.array([I0]), IL=np.array([IL]), return_meta_dict=True)
+    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
+    
+    # Can handle ideal series and shunt
+    V_oc = nNsVth*(np.log(IL + I0) - np.log(I0))
+    V_out = pvsystem.v_from_i_alt(Rsh=np.inf, Rs=0., nNsVth=nNsVth, I=np.array([IL, IL/2., 0.]), I0=I0, IL=IL)
+    V_expected = np.array([0., nNsVth*(np.log(IL - IL/2. + I0) - np.log(I0)), V_oc])
+    assert(isinstance(V_out, type(V_expected)))
+    assert_allclose(V_out, V_expected)
+    _, meta_dict = pvsystem.v_from_i_alt(Rsh=np.inf, Rs=0., nNsVth=nNsVth, I=np.array([IL, IL/2., 0.]), I0=I0, IL=IL, return_meta_dict=True)
+    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
+    
+    # Can handle only ideal series resistance
+    V_out = pvsystem.v_from_i_alt(Rsh=Rsh, Rs=0., nNsVth=nNsVth, I=I, I0=I0, IL=IL)
+    assert(isinstance(V_out, np.ndarray))
+    assert(np.isfinite(V_out)) # No exact expression to evaluate
+    _, meta_dict = pvsystem.v_from_i_alt(Rsh=Rsh, Rs=0., nNsVth=nNsVth, I=I, I0=I0, IL=IL, return_meta_dict=True)
+    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
+    
+    # TODO Stability as Rs->0^+ and/or Rsh->inf
 
 
 @requires_scipy

From 82217fd91f93248432cf0907c895c2568c13d9ec Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Tue, 20 Jun 2017 22:05:14 -0600
Subject: [PATCH 05/22] Use test fixtures

---
 pvlib/pvsystem.py           |  66 ++++---
 pvlib/test/test_pvsystem.py | 370 +++++++++++++++++++++++-------------
 2 files changed, 275 insertions(+), 161 deletions(-)

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index c130511cca..044fa80989 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1953,26 +1953,28 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
 
 def current_sum_at_diode_node(V, I, IL, I0, nNsVth, Rs, Rsh):
     '''
+    TODO Description
+
     Parameters
     ----------
-    V : numpy.ndarray or numpy.float64 or python scalar
+    V : float64 numpy.ndarray, scalar, or pandas.series
         Device voltage [V]
-    I : numpy.ndarray or numpy.float64 or python scalar
+    I : float64 numpy.ndarray, scalar, or pandas.series
         Device current [A]
-    IL : numpy.ndarray or numpy.float64 or python scalar
+    IL : float64 numpy.ndarray, scalar, or pandas.series
         Device photocurrent [A]
-    I0 : numpy.ndarray or numpy.float64 or python scalar
+    I0 : float64 numpy.ndarray, scalar, or pandas.series
         Device saturation current [A]
-    nNsVth : numpy.ndarray or numpy.float64 or python scalar
+    nNsVth : float64 numpy.ndarray, scalar, or pandas.series
         Device thermal voltage [V]
-    Rs : numpy.ndarray or numpy.float64 or python scalar
+    Rs : float64 numpy.ndarray, scalar, or pandas.series
         Device series resistance [Ohm]
-    Rsh : numpy.ndarray or numpy.float64 or python scalar
+    Rsh : float64 numpy.ndarray, scalar, or pandas.series
         Device shunt resistance [Ohm]
     
     Returns
     -------
-    current_sum_at_diode_node : numpy.ndarray or numpy.float64
+    current_sum_at_diode_node : float64 numpy.ndarray, scalar, or pandas.series
         Sum of currents at the diode node in equivalent circuit model [A]
     '''
     
@@ -1984,30 +1986,32 @@ def current_sum_at_diode_node(V, I, IL, I0, nNsVth, Rs, Rsh):
 
 def v_from_i_alt(I, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
     '''
+    TODO Description
+
     Parameters
     ----------
-    I : numpy.ndarray or scalar
+    I : float64 numpy.ndarray or scalar
         Device current [A]
-    IL : numpy.ndarray or scalar
+    IL : float64 numpy.ndarray or scalar
         Device photocurrent [A]
-    I0 : numpy.ndarray or scalar
+    I0 : float64 numpy.ndarray or scalar
         Device saturation current [A]
-    nNsVth : numpy.ndarray or scalar
+    nNsVth : float64 numpy.ndarray or scalar
         Device thermal voltage [V]
-    Rs : numpy.ndarray or scalar
+    Rs : float64 numpy.ndarray or scalar
         Device series resistance [Ohm]
-    Rsh : numpy.ndarray or scalar
+    Rsh : float64 numpy.ndarray or scalar
         Device shunt resistance [Ohm]
-    return_meta_dict : boolean
+    return_meta_dict : boolean scalar
         Return additional computation metadata dictionary
     
     Returns
     -------
-    V : numpy.ndarray
+    V : float64 numpy.ndarray
         Device voltage [V]
     meta_dict : python dictionary (optional, returned when return_meta_dict=True)
         Metadata for computation
-            meta_dict['current_sum_at_diode_node'] : numpy.ndarray like V
+            meta_dict['current_sum_at_diode_node'] : float64 numpy.ndarray like V
                 Sum of currents at the diode node in equivalent circuit model [A]
             meta_dict['inf_Rsh_idx'] : boolean numpy.ndarray like V
                 Indices where infinite shunt resistance gives best solution
@@ -2023,7 +2027,7 @@ def v_from_i_alt(I, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
     Rs = np.asarray(Rs, np.float64)
     Rsh = np.asarray(Rsh, np.float64)
   
-    # Default computation of V using Rsh=np.full_like(Rsh, np.inf), does not lose Rsh shape or type info
+    # Default computation of V using zero_term keeps Rsh shape info
     zero_term = np.zeros_like(I*Rs/Rsh)
     V = nNsVth*(np.log(IL - I - zero_term + I0) - np.log(I0)) - I*Rs
     current_sum_at_diode_node_out = current_sum_at_diode_node(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
@@ -2040,37 +2044,39 @@ def v_from_i_alt(I, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
     current_sum_at_diode_node_out = np.where(finite_Rsh_idx, current_sum_at_diode_node_LambertW, current_sum_at_diode_node_out)
         
     if return_meta_dict:
-        return V, {'current_sum_at_diode_node' : current_sum_at_diode_node_out, 'inf_Rsh_idx' : np.logical_not(finite_Rsh_idx)}
+        return V, {'current_sum_at_diode_node' : current_sum_at_diode_node_out, 'inf_Rsh_idx' : np.array(np.logical_not(finite_Rsh_idx))}
     else:
         return V
 
 
 def i_from_v_alt(V, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
     '''
+    TODO Description
+    
     Parameters
     ----------
-    V : numpy.ndarray or scalar
+    V : float64 numpy.ndarray or scalar
         Device voltage [V]
-    IL : numpy.ndarray or scalar
+    IL : float64 numpy.ndarray or scalar
         Device photocurrent [A]
-    I0 : numpy.ndarray or scalar
+    I0 : float64 numpy.ndarray or scalar
         Device saturation current [A]
-    nNsVth : numpy.ndarray or scalar
+    nNsVth : float64 numpy.ndarray or scalar
         Device thermal voltage [V]
-    Rs : numpy.ndarray or scalar
+    Rs : float64 numpy.ndarray or scalar
         Device series resistance [Ohm]
-    Rsh : numpy.ndarray or scalar
+    Rsh : float64 numpy.ndarray or scalar
         Device shunt resistance [Ohm]
-    return_meta_dict : boolean
+    return_meta_dict : boolean scalar
         Return additional computation metadata dictionary
     
     Returns
     -------
-    I : numpy.ndarray
+    I : float64 numpy.ndarray
         Device current [A]
     meta_dict : python dictionary (optional, returned when return_meta_dict=True)
         Metadata for computation
-            meta_dict['current_sum_at_diode_node'] : numpy.ndarray like I
+            meta_dict['current_sum_at_diode_node'] : float64 numpy.ndarray like I
                 Sum of currents at the diode node in equivalent circuit model [A]
             meta_dict['zero_Rs_idx'] : boolean numpy.ndarray like I
                 Indices where zero series resistance gives best solution
@@ -2086,7 +2092,7 @@ def i_from_v_alt(V, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
     Rs = np.asarray(Rs, np.float64)
     Rsh = np.asarray(Rsh, np.float64)
 
-    # Default computation of I using Rs=np.full_like(Rs, 0.), does not lose Rs shape or type info
+    # Default computation of I using zero_term keeps Rs shape info
     zero_term = np.zeros_like(Rs)
     I = IL - I0*np.expm1((V + zero_term)/nNsVth) - (V + zero_term)/Rsh
     current_sum_at_diode_node_out = current_sum_at_diode_node(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
@@ -2103,7 +2109,7 @@ def i_from_v_alt(V, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
     current_sum_at_diode_node_out = np.where(nonzero_Rs_idx, current_sum_at_diode_node_LambertW, current_sum_at_diode_node_out)
     
     if return_meta_dict:
-        return I, {'current_sum_at_diode_node' : current_sum_at_diode_node_out, 'zero_Rs_idx' : np.logical_not(nonzero_Rs_idx)}
+        return I, {'current_sum_at_diode_node' : current_sum_at_diode_node_out, 'zero_Rs_idx' : np.array(np.logical_not(nonzero_Rs_idx))}
     else:
         return I
 
diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index 848c07ce22..19c15884f3 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -398,150 +398,258 @@ def test_current_sum_at_diode_node():
     assert_array_equal(results_inf_Rsh_1, results_inf_Rsh_2)
 
 
+@pytest.fixture(params=[
+    { # Can handle all python scalar inputs
+     'Rsh' : 20.,
+     'Rs' : 0.1,
+     'nNsVth' : 0.5,
+     'I' : 3.,
+     'I0' : 6.e-7,
+     'IL' : 7.,
+     'V_expected' : np.array(7.5049875193450521),
+     'current_sum_at_diode_node_expected' : np.array(0.),
+     'inf_Rsh_idx_expected' : np.array(False)
+    },
+    { # Can handle all rank-0 array inputs
+     'Rsh' : np.array(20.),
+     'Rs' : np.array(0.1),
+     'nNsVth' : np.array(0.5),
+     'I' : np.array(3.),
+     'I0' : np.array(6.e-7),
+     'IL' : np.array(7.),
+     'V_expected' : np.array(7.5049875193450521),
+     'current_sum_at_diode_node_expected' : np.array(0.),
+     'inf_Rsh_idx_expected' : np.array(False)
+    },
+    { # Can handle all rank-1 singleton array inputs
+     'Rsh' : np.array([20.]),
+     'Rs' : np.array([0.1]),
+     'nNsVth' : np.array([0.5]),
+     'I' : np.array([3.]),
+     'I0' : np.array([6.e-7]),
+     'IL' : np.array([7.]),
+     'V_expected' : np.array([7.5049875193450521]),
+     'current_sum_at_diode_node_expected' : np.array([0.]),
+     'inf_Rsh_idx_expected' : np.array([False])
+    },
+    { # Can handle all rank-1 non-singleton array inputs with infinite shunt
+      #  resistance, Rsh=inf gives V=Voc=nNsVth*(np.log(IL + I0) - np.log(I0) 
+      #  at I=0
+     'Rsh' : np.array([np.inf, 20.]),
+     'Rs' : np.array([0.1, 0.1]),
+     'nNsVth' : np.array([0.5, 0.5]),
+     'I' : np.array([0., 3.]),
+     'I0' : np.array([6.e-7, 6.e-7]),
+     'IL' : np.array([7., 7.]),
+     'V_expected' : np.array([0.5*(np.log(7. + 6.e-7) - np.log(6.e-7)), 7.5049875193450521]),
+     'current_sum_at_diode_node_expected' : np.array([0., 0.]),
+     'inf_Rsh_idx_expected' : np.array([True, False])
+    },
+    { # Can handle mixed inputs with a rank-2 array with infinite shunt
+      #  resistance, Rsh=inf gives V=Voc=nNsVth*(np.log(IL + I0) - np.log(I0) 
+      #  at I=0
+     'Rsh' : np.array([[np.inf, np.inf], [np.inf, np.inf]]),
+     'Rs' : np.array([0.1]),
+     'nNsVth' : np.array(0.5),
+     'I' : 0.,
+     'I0' : np.array([6.e-7]),
+     'IL' : np.array([7.]),
+     'V_expected' : 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))*np.ones((2,2)),
+     'current_sum_at_diode_node_expected' : np.array([[0., 0.], [0., 0.]]),
+     'inf_Rsh_idx_expected' : np.array([[True, True], [True, True]])
+    },
+    { # Can handle ideal series and shunt, Rsh=inf and Rs=0 give 
+      #  V = nNsVth*(np.log(IL - I + I0) - np.log(I0))
+     'Rsh' : np.inf,
+     'Rs' : 0.,
+     'nNsVth' : 0.5,
+     'I' : np.array([7., 7./2., 0.]),
+     'I0' : 6.e-7,
+     'IL' : 7.,
+     'V_expected' : np.array([0., 0.5*(np.log(7. - 7./2. + 6.e-7) - np.log(6.e-7)), 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))]),
+     'current_sum_at_diode_node_expected' : np.array([0., 0., 0.]),
+     'inf_Rsh_idx_expected' : np.array([True, True, True])
+    },
+    { # Can handle only ideal series resistance, no closed form solution
+     'Rsh' : 20.,
+     'Rs' : 0.,
+     'nNsVth' : 0.5,
+     'I' : 3.,
+     'I0' : 6.e-7,
+     'IL' : 7.,
+     'V_expected' : pvsystem.v_from_i_alt(Rsh=20., Rs=0., nNsVth=0.5, I=3., I0=6.e-7, IL=7.),
+     'current_sum_at_diode_node_expected' : np.array(0.),
+     'inf_Rsh_idx_expected' : np.array(False)
+    },
+    { # Can handle all python scalar inputs with big LambertW arg
+     'Rsh' : 500.,
+     'Rs' : 10.,
+     'nNsVth' : 4.06,
+     'I' : 0.,
+     'I0' : 6.e-10,
+     'IL' : 1.2,
+     'V_expected' : np.array(86.320000493521079),
+     'current_sum_at_diode_node_expected' : np.array(0.),
+     'inf_Rsh_idx_expected' : np.array(False)
+    },
+    { # Can handle all python scalar inputs with bigger LambertW arg
+      #  1000 W/m^2 on a Canadian Solar 220M with 20 C ambient temp
+      #  github issue 225
+     'Rsh' : 190.,
+     'Rs' : 1.065,
+     'nNsVth' : 2.89,
+     'I' : 0.,
+     'I0' : 7.05196029e-08,
+     'IL' : 10.491262,
+     'V_expected' : np.array(54.303958833791455),
+     'current_sum_at_diode_node_expected' : np.array(0.),
+     'inf_Rsh_idx_expected' : np.array(False)
+    }])
+def fixture_v_from_i_alt(request):
+    return request.param
+
 @requires_scipy
-def test_i_from_v_alt():
-    # Solution set of Python scalars
-    Rsh = 20.
-    Rs = 0.1
-    nNsVth = 0.5
-    V = 40.
-    I0 = 6.e-7
-    IL = 7.
-    I = -299.746389916
+def test_v_from_i_alt(fixture_v_from_i_alt):
+    # Solution set from fixture
+    Rsh = fixture_v_from_i_alt['Rsh']
+    Rs = fixture_v_from_i_alt['Rs']
+    nNsVth = fixture_v_from_i_alt['nNsVth']
+    I = fixture_v_from_i_alt['I']
+    I0 = fixture_v_from_i_alt['I0']
+    IL = fixture_v_from_i_alt['IL']
+    V_expected = fixture_v_from_i_alt['V_expected']
+    current_sum_at_diode_node_expected = fixture_v_from_i_alt['current_sum_at_diode_node_expected']
+    inf_Rsh_idx_expected = fixture_v_from_i_alt['inf_Rsh_idx_expected']
     
     # Convergence criteria
     atol = 1.e-11
     
-    # Can handle all python scalar inputs
-    I_out = pvsystem.i_from_v_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL)
-    I_expected = np.array(I)
-    assert(isinstance(I_out, type(I_expected)))
-    assert_allclose(I_out, I_expected)
-    _, meta_dict = pvsystem.i_from_v_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL, return_meta_dict=True)
-    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
-
-    # Can handle all rank-0 array inputs
-    I_out = pvsystem.i_from_v_alt(Rsh=np.array(Rsh), Rs=np.array(Rs), nNsVth=np.array(nNsVth), V=np.array(V), I0=np.array(I0), IL=np.array(IL))
-    I_expected = np.array(I)
-    assert(isinstance(I_out, type(I_expected)))
-    assert_allclose(I_out, I_expected)
-    _, meta_dict = pvsystem.i_from_v_alt(Rsh=np.array(Rsh), Rs=np.array(Rs), nNsVth=np.array(nNsVth), V=np.array(V), I0=np.array(I0), IL=np.array(IL), return_meta_dict=True)
-    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
-
-    # Can handle all rank-1 singleton array inputs
-    I_out = pvsystem.i_from_v_alt(Rsh=np.array([Rsh]), Rs=np.array([Rs]), nNsVth=np.array([nNsVth]), V=np.array([V]), I0=np.array([I0]), IL=np.array([IL]))
-    I_expected = np.array([I])
-    assert(isinstance(I_out, type(I_expected)))
-    assert_allclose(I_out, I_expected)
-    _, meta_dict = pvsystem.i_from_v_alt(Rsh=np.array([Rsh]), Rs=np.array([Rs]), nNsVth=np.array([nNsVth]), V=np.array([V]), I0=np.array([I0]), IL=np.array([IL]), return_meta_dict=True)
-    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
-
-    # Can handle all rank-1 non-singleton array inputs with a zero series resistance (Rs=0 gives I=IL=Isc at V=0)
-    I_out = pvsystem.i_from_v_alt(Rsh=np.array([Rsh, Rsh]), Rs=np.array([0., Rs]), nNsVth=np.array([nNsVth, nNsVth]), V=np.array([0., V]), I0=np.array([I0, I0]), IL=np.array([IL, IL]))
-    I_expected = np.array([IL, I])
-    assert(isinstance(I_out, type(I_expected)))
-    assert_allclose(I_out, I_expected)
-    _, meta_dict = pvsystem.i_from_v_alt(Rsh=np.array([Rsh, Rsh]), Rs=np.array([0., Rs]), nNsVth=np.array([nNsVth, nNsVth]), V=np.array([0., V]), I0=np.array([I0, I0]), IL=np.array([IL, IL]), return_meta_dict=True)    
-    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
-    
-    # Can handle mixed inputs with a rank-2 array with zero series resistance (Rs=0 gives I=IL=Isc at V=0)
-    I_out = pvsystem.i_from_v_alt(Rsh=np.array([Rsh]), Rs=np.array([[0., 0.], [0., 0.]]), nNsVth=np.array(nNsVth), V=np.array(0.), I0=np.array([I0]), IL=np.array([IL]))
-    I_expected = np.array([[IL, IL], [IL, IL]])
-    assert(isinstance(I_out, type(I_expected)))
-    assert_allclose(I_out, I_expected)
-    _, meta_dict = pvsystem.i_from_v_alt(Rsh=np.array([Rsh]), Rs=np.array([[0., 0.], [0., 0.]]), nNsVth=np.array(nNsVth), V=np.array(0.), I0=np.array([I0]), IL=np.array([IL]), return_meta_dict=True)
-    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
-    
-    # Can handle ideal series and shunt
-    V_oc = nNsVth*(np.log(IL + I0) - np.log(I0))
-    I_out = pvsystem.i_from_v_alt(Rsh=np.inf, Rs=0., nNsVth=nNsVth, V=np.array([0., V_oc/2., V_oc]), I0=I0, IL=IL)
-    I_expected = np.array([IL, IL - I0*np.expm1(V_oc/2./nNsVth), IL - I0*np.expm1(V_oc/nNsVth)])
-    assert(isinstance(I_out, type(I_expected)))
-    assert_allclose(I_out, I_expected)
-    _, meta_dict = pvsystem.i_from_v_alt(Rsh=np.inf, Rs=0., nNsVth=nNsVth, V=np.array([0., V_oc/2., V_oc]), I0=I0, IL=IL, return_meta_dict=True)
-    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
-    
-    # Can handle only ideal shunt resistance
-    I_out = pvsystem.i_from_v_alt(Rsh=np.inf, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL)
-    assert(isinstance(I_out, np.ndarray))
-    assert(np.isfinite(I_out)) # No exact expression to evaluate
-    _, meta_dict = pvsystem.i_from_v_alt(Rsh=np.inf, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL, return_meta_dict=True)
-    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
+    V = pvsystem.v_from_i_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, I=I, I0=I0, IL=IL)
+    assert(isinstance(V, type(V_expected)))
+    assert(isinstance(V.dtype, type(V_expected.dtype)))
+    assert_allclose(V, V_expected, atol=atol)
+    _, meta_dict = pvsystem.v_from_i_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, I=I, I0=I0, IL=IL, return_meta_dict=True)
+    assert(isinstance(meta_dict['current_sum_at_diode_node'], type(V)))
+    assert(isinstance(meta_dict['current_sum_at_diode_node'].dtype, type(V.dtype)))
+    assert_allclose(meta_dict['current_sum_at_diode_node'], current_sum_at_diode_node_expected, atol=atol)
+    assert(isinstance(meta_dict['inf_Rsh_idx'], type(V)))
+    assert(isinstance(meta_dict['inf_Rsh_idx'].dtype, type(V.dtype)))
+    assert_array_equal(meta_dict['inf_Rsh_idx'], inf_Rsh_idx_expected)
     
     # TODO Stability as Rs->0^+ and/or Rsh->inf
 
 
+@pytest.fixture(params=[
+    { # Can handle all python scalar inputs
+     'Rsh' : 20.,
+     'Rs' : 0.1,
+     'nNsVth' : 0.5,
+     'V' : 40.,
+     'I0' : 6.e-7,
+     'IL' : 7.,
+     'I_expected' : np.array(-299.746389916),
+     'current_sum_at_diode_node_expected' : np.array(0.),
+     'zero_Rs_idx_expected' : np.array(False)
+    },
+    { # Can handle all rank-0 array inputs
+     'Rsh' : np.array(20.),
+     'Rs' : np.array(0.1),
+     'nNsVth' : np.array(0.5),
+     'V' : np.array(40.),
+     'I0' : np.array(6.e-7),
+     'IL' : np.array(7.),
+     'I_expected' : np.array(-299.746389916),
+     'current_sum_at_diode_node_expected' : np.array(0.),
+     'zero_Rs_idx_expected' : np.array(False)
+    },
+    { # Can handle all rank-1 singleton array inputs
+     'Rsh' : np.array([20.]),
+     'Rs' : np.array([0.1]),
+     'nNsVth' : np.array([0.5]),
+     'V' : np.array([40.]),
+     'I0' : np.array([6.e-7]),
+     'IL' : np.array([7.]),
+     'I_expected' : np.array([-299.746389916]),
+     'current_sum_at_diode_node_expected' : np.array([0.]),
+     'zero_Rs_idx_expected' : np.array([False])
+    },
+    { # Can handle all rank-1 non-singleton array inputs with a zero 
+      #  series resistance, Rs=0 gives I=IL=Isc at V=0
+     'Rsh' : np.array([20., 20.]),
+     'Rs' : np.array([0., 0.1]),
+     'nNsVth' : np.array([0.5, 0.5]),
+     'V' : np.array([0., 40.]),
+     'I0' : np.array([6.e-7, 6.e-7]),
+     'IL' : np.array([7., 7.]),
+     'I_expected' : np.array([7., -299.746389916]),
+     'current_sum_at_diode_node_expected' : np.array([0., 0.]),
+     'zero_Rs_idx_expected' : np.array([True, False])
+    },
+    { # Can handle mixed inputs with a rank-2 array with zero series 
+      #  resistance, Rs=0 gives I=IL=Isc at V=0
+     'Rsh' : np.array([20.]),
+     'Rs' : np.array([[0., 0.], [0., 0.]]),
+     'nNsVth' : np.array(0.5),
+     'V' : 0.,
+     'I0' : np.array([6.e-7]),
+     'IL' : np.array([7.]),
+     'I_expected' : np.array([[7., 7.], [7., 7.]]),
+     'current_sum_at_diode_node_expected' : np.array([[0., 0.], [0., 0.]]),
+     'zero_Rs_idx_expected' : np.array([[True, True], [True, True]])
+    },
+    { # Can handle ideal series and shunt, Rsh=inf and Rs=0 give 
+      #  V_oc = nNsVth*(np.log(IL + I0) - np.log(I0))
+     'Rsh' : np.inf,
+     'Rs' : 0.,
+     'nNsVth' : 0.5,
+     'V' : np.array([0., 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))/2., 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))]),
+     'I0' : 6.e-7,
+     'IL' : 7.,
+     'I_expected' : np.array([7., 7. - 6.e-7*np.expm1((np.log(7. + 6.e-7) - np.log(6.e-7))/2.), 0.]),
+     'current_sum_at_diode_node_expected' : np.array([0., 0., 0.]),
+     'zero_Rs_idx_expected' : np.array([True, True, True])
+    },
+    { # Can handle only ideal shunt resistance, no closed form solution
+     'Rsh' : np.inf,
+     'Rs' : 0.1,
+     'nNsVth' : 0.5,
+     'V' : 40.,
+     'I0' : 6.e-7,
+     'IL' : 7.,
+     'I_expected' : pvsystem.i_from_v_alt(Rsh=np.inf, Rs=0.1, nNsVth=0.5, V=40., I0=6.e-7, IL=7.),
+     'current_sum_at_diode_node_expected' : np.array(0.),
+     'zero_Rs_idx_expected' : np.array(False)
+    }])
+def fixture_i_from_v_alt(request):
+    return request.param
+
 @requires_scipy
-def test_v_from_i_alt():
-    # Solution set of Python scalars
-    Rsh = 20.
-    Rs = 0.1
-    nNsVth = 0.5
-    I = 3.
-    I0 = 6.e-7
-    IL = 7.
-    V = 7.5049875193450521
-    
+def test_i_from_v_alt(fixture_i_from_v_alt):
+    # Solution set from fixture
+    Rsh = fixture_i_from_v_alt['Rsh']
+    Rs = fixture_i_from_v_alt['Rs']
+    nNsVth = fixture_i_from_v_alt['nNsVth']
+    V = fixture_i_from_v_alt['V']
+    I0 = fixture_i_from_v_alt['I0']
+    IL = fixture_i_from_v_alt['IL']
+    I_expected = fixture_i_from_v_alt['I_expected']
+    current_sum_at_diode_node_expected = fixture_i_from_v_alt['current_sum_at_diode_node_expected']
+    zero_Rs_idx_expected = fixture_i_from_v_alt['zero_Rs_idx_expected']
+
     # Convergence criteria
     atol = 1.e-11
     
-    # Can handle all python scalar inputs
-    V_out = pvsystem.v_from_i_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, I=I, I0=I0, IL=IL)
-    V_expected = np.array(V)
-    assert(isinstance(V_out, type(V_expected)))
-    assert_allclose(V_out, V_expected)
-    _, meta_dict = pvsystem.v_from_i_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, I=I, I0=I0, IL=IL, return_meta_dict=True)
-    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
-
-    # Can handle all rank-0 array inputs
-    V_out = pvsystem.v_from_i_alt(Rsh=np.array(Rsh), Rs=np.array(Rs), nNsVth=np.array(nNsVth), I=np.array(I), I0=np.array(I0), IL=np.array(IL))
-    V_expected = np.array(V)
-    assert(isinstance(V_out, type(V_expected)))
-    assert_allclose(V_out, V_expected)
-    _, meta_dict = pvsystem.v_from_i_alt(Rsh=np.array(Rsh), Rs=np.array(Rs), nNsVth=np.array(nNsVth), I=np.array(I), I0=np.array(I0), IL=np.array(IL), return_meta_dict=True)
-    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
-
-    # Can handle all rank-1 singleton array inputs
-    V_out = pvsystem.v_from_i_alt(Rsh=np.array([Rsh]), Rs=np.array([Rs]), nNsVth=np.array([nNsVth]), I=np.array([I]), I0=np.array([I0]), IL=np.array([IL]))
-    V_expected = np.array([V])
-    assert(isinstance(V_out, type(V_expected)))
-    assert_allclose(V_out, V_expected)
-    _, meta_dict = pvsystem.v_from_i_alt(Rsh=np.array([Rsh]), Rs=np.array([Rs]), nNsVth=np.array([nNsVth]), I=np.array([I]), I0=np.array([I0]), IL=np.array([IL]), return_meta_dict=True)
-    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
-
-    # Can handle all rank-1 non-singleton array inputs with infinite shunt resistance (Rsh=inf gives V=Voc at I=0)
-    V_oc = nNsVth*(np.log(IL + I0) - np.log(I0))
-    V_out = pvsystem.v_from_i_alt(Rsh=np.array([np.inf, Rsh]), Rs=np.array([Rs, Rs]), nNsVth=np.array([nNsVth, nNsVth]), I=np.array([0., I]), I0=np.array([I0, I0]), IL=np.array([IL, IL]))
-    V_expected = np.array([V_oc, V])
-    assert(isinstance(V_out, type(V_expected)))
-    assert_allclose(V_out, V_expected)
-    _, meta_dict = pvsystem.v_from_i_alt(Rsh=np.array([np.inf, Rsh]), Rs=np.array([Rs, Rs]), nNsVth=np.array([nNsVth, nNsVth]), I=np.array([0., I]), I0=np.array([I0, I0]), IL=np.array([IL, IL]), return_meta_dict=True)    
-    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
-    
-    # Can handle mixed inputs with a rank-2 array with infinite shunt resistance (Rsh=inf gives V=Voc at I=0)
-    V_oc = nNsVth*(np.log(IL + I0) - np.log(I0))
-    V_out = pvsystem.v_from_i_alt(Rsh=np.array([[np.inf, np.inf], [np.inf, np.inf]]), Rs=np.array([Rs]), nNsVth=np.array(nNsVth), I=np.array(0.), I0=np.array([I0]), IL=np.array([IL]))
-    V_expected = np.array([[V_oc, V_oc], [V_oc, V_oc]])
-    assert(isinstance(V_out, type(V_expected)))
-    assert_allclose(V_out, V_expected)
-    _, meta_dict = pvsystem.v_from_i_alt(Rsh=np.array([[np.inf, np.inf], [np.inf, np.inf]]), Rs=np.array([Rs]), nNsVth=np.array(nNsVth), I=np.array(0.), I0=np.array([I0]), IL=np.array([IL]), return_meta_dict=True)
-    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
-    
-    # Can handle ideal series and shunt
-    V_oc = nNsVth*(np.log(IL + I0) - np.log(I0))
-    V_out = pvsystem.v_from_i_alt(Rsh=np.inf, Rs=0., nNsVth=nNsVth, I=np.array([IL, IL/2., 0.]), I0=I0, IL=IL)
-    V_expected = np.array([0., nNsVth*(np.log(IL - IL/2. + I0) - np.log(I0)), V_oc])
-    assert(isinstance(V_out, type(V_expected)))
-    assert_allclose(V_out, V_expected)
-    _, meta_dict = pvsystem.v_from_i_alt(Rsh=np.inf, Rs=0., nNsVth=nNsVth, I=np.array([IL, IL/2., 0.]), I0=I0, IL=IL, return_meta_dict=True)
-    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
-    
-    # Can handle only ideal series resistance
-    V_out = pvsystem.v_from_i_alt(Rsh=Rsh, Rs=0., nNsVth=nNsVth, I=I, I0=I0, IL=IL)
-    assert(isinstance(V_out, np.ndarray))
-    assert(np.isfinite(V_out)) # No exact expression to evaluate
-    _, meta_dict = pvsystem.v_from_i_alt(Rsh=Rsh, Rs=0., nNsVth=nNsVth, I=I, I0=I0, IL=IL, return_meta_dict=True)
-    assert_allclose(meta_dict['current_sum_at_diode_node'], 0., atol=atol)
+    I = pvsystem.i_from_v_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL)
+    assert(isinstance(I, type(I_expected)))
+    assert(isinstance(I.dtype, type(I_expected.dtype)))
+    assert_allclose(I, I_expected, atol=atol)
+    _, meta_dict = pvsystem.i_from_v_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL, return_meta_dict=True)
+    assert(isinstance(meta_dict['current_sum_at_diode_node'], type(I)))
+    assert(isinstance(meta_dict['current_sum_at_diode_node'].dtype, type(I.dtype)))
+    assert_allclose(meta_dict['current_sum_at_diode_node'], current_sum_at_diode_node_expected, atol=atol)
+    assert(isinstance(meta_dict['zero_Rs_idx'], type(I)))
+    assert(isinstance(meta_dict['zero_Rs_idx'].dtype, type(I.dtype)))
+    assert_array_equal(meta_dict['zero_Rs_idx'], zero_Rs_idx_expected)
     
     # TODO Stability as Rs->0^+ and/or Rsh->inf
 

From 135ce4c93bbba37ad1de5eb689664d5735933558 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Fri, 23 Jun 2017 12:03:34 -0600
Subject: [PATCH 06/22] Add @requires_scipy to test fixtures

---
 pvlib/test/test_pvsystem.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index 19c15884f3..6114f21c29 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -397,7 +397,7 @@ def test_current_sum_at_diode_node():
     results_inf_Rsh_2 = IL - I0*np.expm1((V + I*Rs)/nNsVth) - I
     assert_array_equal(results_inf_Rsh_1, results_inf_Rsh_2)
 
-
+@requires_scipy
 @pytest.fixture(params=[
     { # Can handle all python scalar inputs
      'Rsh' : 20.,
@@ -539,6 +539,7 @@ def test_v_from_i_alt(fixture_v_from_i_alt):
     # TODO Stability as Rs->0^+ and/or Rsh->inf
 
 
+@requires_scipy
 @pytest.fixture(params=[
     { # Can handle all python scalar inputs
      'Rsh' : 20.,

From f927aef87071cd0bb7bd556dbdce2b39fbf27a7b Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Fri, 23 Jun 2017 13:19:24 -0600
Subject: [PATCH 07/22] More current_sum_at_diode_node() tests and using
 fixture

---
 pvlib/test/test_pvsystem.py | 155 +++++++++++++++++++++++++-----------
 1 file changed, 107 insertions(+), 48 deletions(-)

diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index 6114f21c29..39798d0958 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -353,51 +353,111 @@ def test_PVSystem_calcparams_desoto(cec_module_params):
     assert_allclose(nNsVth, 0.473)
 
 
-def test_current_sum_at_diode_node():
-    V = np.array([40., 0. , 0])
-    I = np.array([0., 0., 3.])
-    IL = 7.
-    I0 = 6.e-7
-    nNsVth = 0.5 
-    Rs = 0.1
-    Rsh = 20.
-    
-    results_1 = np.full_like(V, np.nan)
-    results_2 = np.full_like(V, np.nan)
-    
-    results_1[0] = pvsystem.current_sum_at_diode_node(V=V[0], I=I[0], IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
-    results_2[0] = IL - I0*np.expm1(V[0]/nNsVth) - V[0]/Rsh
-    
-    results_1[1] = pvsystem.current_sum_at_diode_node(V=V[1], I=I[1], IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
-    results_2[1] = IL
-    
-    results_1[2] = pvsystem.current_sum_at_diode_node(V=V[2], I=I[2], IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
-    results_2[2] = IL - I0*np.expm1(I[2]*Rs/nNsVth) - I[2]*Rs/Rsh - I[2]
-    
-    assert_array_equal(results_1, results_2)
+@pytest.fixture(params=[# Not necessarily I-V curve solutions
+    { # Can handle all python scalar inputs
+     'V' : 40.,
+     'I' : 3.,
+     'IL' : 7.,
+     'I0' : 6.e-7,
+     'nNsVth' : 0.5,
+     'Rs' : 0.1,
+     'Rsh' : 20.,
+     'current_sum_at_diode_node_expected' : np.float64(7. - 6.e-7*np.expm1((40. + 3.*0.1)/0.5) - (40. + 3.*0.1)/20. - 3.)
+    },
+    { # Can handle all rank-0 array inputs
+     'V' : np.array(40.),
+     'I' : np.array(3.),
+     'IL' : np.array(7.),
+     'I0' : np.array(6.e-7),
+     'nNsVth' : np.array(0.5),
+     'Rs' : np.array(0.1),
+     'Rsh' : np.array(20.),
+     'current_sum_at_diode_node_expected' : np.float64(7. - 6.e-7*np.expm1((40. + 3.*0.1)/0.5) - (40. + 3.*0.1)/20. - 3.)
+    },
+    { # Can handle all rank-1 singleton array inputs
+     'V' : np.array([40.]),
+     'I' : np.array([3.]),
+     'IL' : np.array([7.]),
+     'I0' : np.array([6.e-7]),
+     'nNsVth' : np.array([0.5]),
+     'Rs' : np.array([0.1]),
+     'Rsh' : np.array([20.]),
+     'current_sum_at_diode_node_expected' : np.array([7. - 6.e-7*np.expm1((40. + 3.*0.1)/0.5) - (40. + 3.*0.1)/20. - 3.])
+    },
+    { # Can handle all rank-1 non-singleton array inputs
+     'V' : np.array([40., 0. , 0.]),
+     'I' : np.array([0., 0., 3.]),
+     'IL' : np.array([7., 7., 7.]),
+     'I0' : np.array([6.e-7, 6.e-7, 6.e-7]),
+     'nNsVth' : np.array([0.5, 0.5, 0.5]),
+     'Rs' : np.array([0.1, 0.1, 0.1]),
+     'Rsh' : np.array([20., 20., 20.]),
+     'current_sum_at_diode_node_expected' : np.array([7. - 6.e-7*np.expm1(40./0.5) - 40./0.1, 7., 7. - 6.e-7*np.expm1(3.*0.1/0.5) - 3.*0.1/20. - 3.])
+    },
+    { # Can handle mixed inputs with non-singleton Pandas Series
+     'V' : pd.Series([40., 0. , 0.]),
+     'I' : pd.Series([0., 0., 3.]),
+     'IL' : 7.,
+     'I0' : 6.e-7,
+     'nNsVth' : 0.5,
+     'Rs' : 0.1,
+     'Rsh' : 20.,
+     'current_sum_at_diode_node_expected' : pd.Series([7. - 6.e-7*np.expm1(40./0.5) - 40./0.1, 7., 7. - 6.e-7*np.expm1(3.*0.1/0.5) - 3.*0.1/20. - 3.])
+    },
+    { # Can handle mixed inputs with rank-2 arrays
+     'V' : np.array([[40., 0. , 0.], [0., 0. , 40.]]),
+     'I' : np.array([[0., 0., 3.], [3., 0., 0.]]),
+     'IL' : 7.,
+     'I0' : np.full((1,3), 6.e-7),
+     'nNsVth' : np.array(0.5),
+     'Rs' : np.array([0.1]),
+     'Rsh' : np.full((2,3), 20.),
+     'current_sum_at_diode_node_expected' : np.array([[7. - 6.e-7*np.expm1(40./0.5) - 40./0.1, 7., 7. - 6.e-7*np.expm1(3.*0.1/0.5) - 3.*0.1/20. - 3.], \
+                                                      [ 7. - 6.e-7*np.expm1(3.*0.1/0.5) - 3.*0.1/20. - 3., 7., 7. - 6.e-7*np.expm1(40./0.5) - 40./0.1]])
+    },
+    { # Can handle infinite shunt resistance with positive series resistance
+     'V' : 40.,
+     'I' : 3.,
+     'IL' : 7.,
+     'I0' : 6.e-7,
+     'nNsVth' : 0.5,
+     'Rs' : 0.1,
+     'Rsh' : np.inf,
+     'current_sum_at_diode_node_expected' : np.float64(7. - 6.e-7*np.expm1((40. + 3.*0.1)/0.5) - 3.)
+    },
+    { # Can handle infinite shunt resistance with zero series resistance
+     'V' : 40.,
+     'I' : 3.,
+     'IL' : 7.,
+     'I0' : 6.e-7,
+     'nNsVth' : 0.5,
+     'Rs' : 0.,
+     'Rsh' : np.inf,
+     'current_sum_at_diode_node_expected' : np.float64(7. - 6.e-7*np.expm1(40./0.5) - 3.)
+    }])
+def fixture_current_sum_at_diode_node(request):
+    return request.param
+
+def test_current_sum_at_diode_node(fixture_current_sum_at_diode_node):
+    # Note: The computation of this function is so straight forward that we do
+    #  NOT extensively verify ufunc behavior
     
-    results_vec = pvsystem.current_sum_at_diode_node(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
-    assert_array_equal(results_vec, results_2)
-
-    V = 0.
-    I = 0.
-    nNsVth = np.asarray([0.45, 0.5, 0.55])
-    results_vec = pvsystem.current_sum_at_diode_node(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
-    assert_array_equal(results_vec, np.asarray([IL, IL, IL]))
-
-    nNsVth = pd.Series([0.45, 0.5, 0.55])
-    results_series = pvsystem.current_sum_at_diode_node(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
-    assert_series_equal(results_series, pd.Series([IL, IL, IL]))    
-
-    V = 40.
-    I = 3.
-    nNsVth = 0.5
-    Rsh = np.inf
-    results_inf_Rsh_1 = pvsystem.current_sum_at_diode_node(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
-    results_inf_Rsh_2 = IL - I0*np.expm1((V + I*Rs)/nNsVth) - I
-    assert_array_equal(results_inf_Rsh_1, results_inf_Rsh_2)
+    # Solution set loaded from fixture
+    V = fixture_current_sum_at_diode_node['V']
+    I = fixture_current_sum_at_diode_node['I']
+    IL = fixture_current_sum_at_diode_node['IL']
+    I0 = fixture_current_sum_at_diode_node['I0']
+    nNsVth = fixture_current_sum_at_diode_node['nNsVth']
+    Rs = fixture_current_sum_at_diode_node['Rs']
+    Rsh = fixture_current_sum_at_diode_node['Rsh']
+    current_sum_at_diode_node_expected = fixture_current_sum_at_diode_node['current_sum_at_diode_node_expected']
+
+    current_sum_at_diode_node = pvsystem.current_sum_at_diode_node(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    assert(isinstance(current_sum_at_diode_node, type(current_sum_at_diode_node_expected)))
+    assert(isinstance(current_sum_at_diode_node.dtype, type(current_sum_at_diode_node_expected.dtype)))
+    assert_array_equal(current_sum_at_diode_node, current_sum_at_diode_node_expected)
+
 
-@requires_scipy
 @pytest.fixture(params=[
     { # Can handle all python scalar inputs
      'Rsh' : 20.,
@@ -477,7 +537,7 @@ def test_current_sum_at_diode_node():
      'I' : 3.,
      'I0' : 6.e-7,
      'IL' : 7.,
-     'V_expected' : pvsystem.v_from_i_alt(Rsh=20., Rs=0., nNsVth=0.5, I=3., I0=6.e-7, IL=7.),
+     'V_expected' : np.array(7.804987519345062),
      'current_sum_at_diode_node_expected' : np.array(0.),
      'inf_Rsh_idx_expected' : np.array(False)
     },
@@ -510,7 +570,7 @@ def fixture_v_from_i_alt(request):
 
 @requires_scipy
 def test_v_from_i_alt(fixture_v_from_i_alt):
-    # Solution set from fixture
+    # Solution set loaded from fixture
     Rsh = fixture_v_from_i_alt['Rsh']
     Rs = fixture_v_from_i_alt['Rs']
     nNsVth = fixture_v_from_i_alt['nNsVth']
@@ -539,7 +599,6 @@ def test_v_from_i_alt(fixture_v_from_i_alt):
     # TODO Stability as Rs->0^+ and/or Rsh->inf
 
 
-@requires_scipy
 @pytest.fixture(params=[
     { # Can handle all python scalar inputs
      'Rsh' : 20.,
@@ -617,7 +676,7 @@ def test_v_from_i_alt(fixture_v_from_i_alt):
      'V' : 40.,
      'I0' : 6.e-7,
      'IL' : 7.,
-     'I_expected' : pvsystem.i_from_v_alt(Rsh=np.inf, Rs=0.1, nNsVth=0.5, V=40., I0=6.e-7, IL=7.),
+     'I_expected' : np.array(-299.7383436645412),
      'current_sum_at_diode_node_expected' : np.array(0.),
      'zero_Rs_idx_expected' : np.array(False)
     }])
@@ -626,7 +685,7 @@ def fixture_i_from_v_alt(request):
 
 @requires_scipy
 def test_i_from_v_alt(fixture_i_from_v_alt):
-    # Solution set from fixture
+    # Solution set loaded from fixture
     Rsh = fixture_i_from_v_alt['Rsh']
     Rs = fixture_i_from_v_alt['Rs']
     nNsVth = fixture_i_from_v_alt['nNsVth']

From 741a7f47a7b20330cc33f865788de55ea32f491a Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Fri, 23 Jun 2017 14:27:27 -0600
Subject: [PATCH 08/22] Naming, documentation, and formatting

---
 pvlib/pvsystem.py           | 204 ++++++++++++++++++++++++------------
 pvlib/test/test_pvsystem.py | 146 +++++++++++++-------------
 2 files changed, 208 insertions(+), 142 deletions(-)

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 044fa80989..148e9c92a5 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1851,16 +1851,17 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
     I = np.asarray(current, np.float64)
     I0 = np.asarray(saturation_current, np.float64)
     IL = np.asarray(photocurrent, np.float64)
-    
+
     # This transforms any ideal Rsh=np.inf into Gsh=0., which is generally more numerically stable
     Gsh = 1./Rsh
-    
+
     # argW cannot be float128
     argW = I0/(Gsh*nNsVth)*np.exp((-I + IL + I0)/(Gsh*nNsVth))
     lambertwterm = lambertw(argW).real
 
     # Calculate using log(argW) in case argW is really big
-    logargW = (np.log(I0) - np.log(Gsh) - np.log(nNsVth) + (-I + IL + I0)/(Gsh*nNsVth))
+    logargW = \
+(np.log(I0) - np.log(Gsh) - np.log(nNsVth) + (-I + IL + I0)/(Gsh*nNsVth))
 
     # Three iterations of Newton-Raphson method to solve
     # w+log(w)=logargW. The initial guess is w=logargW. Where direct
@@ -1871,7 +1872,8 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
         w = w * (1 - np.log(w) + logargW) / (1 + w)
     lambertwterm_log = w
 
-    lambertwterm = np.where(np.isfinite(lambertwterm), lambertwterm, lambertwterm_log)
+    lambertwterm = \
+np.where(np.isfinite(lambertwterm), lambertwterm, lambertwterm_log)
 
     # Eqn. 3 in Jain and Kapoor, 2004
     V = (IL + I0 - I)/Gsh - I*Rs - nNsVth*lambertwterm
@@ -1937,88 +1939,122 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
     V = np.asarray(voltage, np.float64)
     I0 = np.asarray(saturation_current, np.float64)
     IL = np.asarray(photocurrent, np.float64)
-    
+
     # This transforms any ideal Rsh=np.inf into Gsh=0., which is generally more numerically stable
     Gsh = 1./Rsh
-    
+
     # argW cannot be float128
-    argW = Rs*I0/(nNsVth*(Rs*Gsh + 1.))*np.exp((Rs*(IL + I0) + V)/(nNsVth*(Rs*Gsh + 1.)))
+    argW = \
+Rs*I0/(nNsVth*(Rs*Gsh + 1.))*np.exp((Rs*(IL + I0) + V)/(nNsVth*(Rs*Gsh + 1.)))
     lambertwterm = lambertw(argW).real
 
     # Eqn. 4 in Jain and Kapoor, 2004
     I = (IL + I0 - V*Gsh)/(Rs*Gsh + 1.) - (nNsVth/Rs)*lambertwterm
-    
+
     return I
 
 
-def current_sum_at_diode_node(V, I, IL, I0, nNsVth, Rs, Rsh):
+def sdm_current_sum(V, I, IL, I0, nNsVth, Rs, Rsh):
     '''
-    TODO Description
+    Computes the sum of currents at the diode node at electrical steady state 
+    using the standard single diode model (SDM) as described in, e.g., Jain and
+    Kapoor 2004 [1]. An ideal device is specified by Rs=0 and Rsh=numpy.inf.
+    This function behaves as a numpy ufunc, including pandas.Series inputs.
 
     Parameters
     ----------
-    V : float64 numpy.ndarray, scalar, or pandas.series
+    V : numeric
         Device voltage [V]
-    I : float64 numpy.ndarray, scalar, or pandas.series
+
+    I : numeric
         Device current [A]
-    IL : float64 numpy.ndarray, scalar, or pandas.series
+
+    IL : numeric
         Device photocurrent [A]
-    I0 : float64 numpy.ndarray, scalar, or pandas.series
+
+    I0 : numeric
         Device saturation current [A]
-    nNsVth : float64 numpy.ndarray, scalar, or pandas.series
+
+    nNsVth : numeric
         Device thermal voltage [V]
-    Rs : float64 numpy.ndarray, scalar, or pandas.series
+
+    Rs : numeric
         Device series resistance [Ohm]
-    Rsh : float64 numpy.ndarray, scalar, or pandas.series
+
+    Rsh : numeric
         Device shunt resistance [Ohm]
-    
+
     Returns
     -------
-    current_sum_at_diode_node : float64 numpy.ndarray, scalar, or pandas.series
+    sdm_current_sum : numeric
         Sum of currents at the diode node in equivalent circuit model [A]
+
+    References
+    ----------
+    [1] A. Jain, A. Kapoor, "Exact analytical solutions of the parameters of 
+    real solar cells using Lambert W-function", Solar Energy Materials and 
+    Solar Cells, 81 (2004) 269-277.
     '''
-    
-    # TODO Ensure that this qualifies as a numpy ufunc
-    
-    # current_sum_at_diode_node
+
+    # sdm_current_sum
     return IL - I0*np.expm1((V + I*Rs)/nNsVth) - (V + I*Rs)/Rsh - I
 
 
-def v_from_i_alt(I, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
+def sdm_v_from_i(I, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
     '''
-    TODO Description
+    Computes the device voltage at the given device current using the standard 
+    single diode model (SDM) as described in, e.g., Jain and Kapoor 2004 [1].
+    An ideal device is specified by Rs=0 and Rsh=numpy.inf and the solution is 
+    per Eq 3 of [1] unless Rsh=numpy.inf gives a more accurate (closed form) 
+    solution as determined by the sum of currents at the diode node, which 
+    should be zero at electrical steady state.
+    Inputs to this function can include scalars and pandas.Series, but it
+    always outputs a float64 numpy.ndarray regardless of input type(s).
 
     Parameters
     ----------
-    I : float64 numpy.ndarray or scalar
+    I : numeric
         Device current [A]
-    IL : float64 numpy.ndarray or scalar
+
+    IL : numeric
         Device photocurrent [A]
-    I0 : float64 numpy.ndarray or scalar
+
+    I0 : numeric
         Device saturation current [A]
-    nNsVth : float64 numpy.ndarray or scalar
+
+    nNsVth : numeric
         Device thermal voltage [V]
-    Rs : float64 numpy.ndarray or scalar
+
+    Rs : numeric
         Device series resistance [Ohm]
-    Rsh : float64 numpy.ndarray or scalar
+
+    Rsh : numeric
         Device shunt resistance [Ohm]
+
     return_meta_dict : boolean scalar
         Return additional computation metadata dictionary
-    
+
     Returns
     -------
     V : float64 numpy.ndarray
         Device voltage [V]
-    meta_dict : python dictionary (optional, returned when return_meta_dict=True)
+
+    meta_dict : dictionary (optional, returned when return_meta_dict=True)
         Metadata for computation
-            meta_dict['current_sum_at_diode_node'] : float64 numpy.ndarray like V
-                Sum of currents at the diode node in equivalent circuit model [A]
+
+            meta_dict['sdm_current_sum'] : float64 numpy.ndarray like V
+                Sum of currents at diode node in equivalent circuit model [A]
+
             meta_dict['inf_Rsh_idx'] : boolean numpy.ndarray like V
                 Indices where infinite shunt resistance gives best solution
+
+    References
+    ----------
+    [1] A. Jain, A. Kapoor, "Exact analytical solutions of the parameters of 
+    real solar cells using Lambert W-function", Solar Energy Materials and 
+    Solar Cells, 81 (2004) 269-277.
     '''
 
-    # TODO Check if this qualifies as a numpy ufunc
-        
     # Ensure inputs are all np.ndarray with np.float64 type
     I = np.asarray(I, np.float64)
     IL = np.asarray(IL, np.float64)
@@ -2026,64 +2062,89 @@ def v_from_i_alt(I, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
     nNsVth = np.asarray(nNsVth, np.float64)
     Rs = np.asarray(Rs, np.float64)
     Rsh = np.asarray(Rsh, np.float64)
-  
+
     # Default computation of V using zero_term keeps Rsh shape info
     zero_term = np.zeros_like(I*Rs/Rsh)
     V = nNsVth*(np.log(IL - I - zero_term + I0) - np.log(I0)) - I*Rs
-    current_sum_at_diode_node_out = current_sum_at_diode_node(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
-    
+    sdm_current_sum_out = \
+sdm_current_sum(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+
     # Computation of V using LambertW for provided Rsh
     V_LambertW = v_from_i(Rsh, Rs, nNsVth, I, I0, IL)
-    current_sum_at_diode_node_LambertW = current_sum_at_diode_node(V=V_LambertW, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
-    
+    sdm_current_sum_LambertW = \
+sdm_current_sum(V=V_LambertW, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+
     # Compute selection indices (may be a scalar boolean)
-    finite_Rsh_idx = np.logical_and(np.isfinite(current_sum_at_diode_node_LambertW), np.absolute(current_sum_at_diode_node_LambertW) <= np.absolute(current_sum_at_diode_node_out))
-    
+    finite_Rsh_idx = np.logical_and(np.isfinite(sdm_current_sum_LambertW), \
+np.absolute(sdm_current_sum_LambertW) <= np.absolute(sdm_current_sum_out))
+
     # These are always np.ndarray
     V = np.where(finite_Rsh_idx, V_LambertW, V)
-    current_sum_at_diode_node_out = np.where(finite_Rsh_idx, current_sum_at_diode_node_LambertW, current_sum_at_diode_node_out)
-        
+    sdm_current_sum_out = \
+np.where(finite_Rsh_idx, sdm_current_sum_LambertW, sdm_current_sum_out)
+
     if return_meta_dict:
-        return V, {'current_sum_at_diode_node' : current_sum_at_diode_node_out, 'inf_Rsh_idx' : np.array(np.logical_not(finite_Rsh_idx))}
+        return V, {'sdm_current_sum' : sdm_current_sum_out, \
+'inf_Rsh_idx' : np.array(np.logical_not(finite_Rsh_idx))}
     else:
         return V
 
 
-def i_from_v_alt(V, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
+def sdm_i_from_v(V, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
     '''
-    TODO Description
-    
+    Computes the device current at the given device voltage using the standard 
+    single diode model (SDM) as described in, e.g., Jain and Kapoor 2004 [1].
+    An ideal device is specified by Rs=0 and Rsh=numpy.inf and the solution is 
+    per Eq 2 of [1] unless Rs=0 gives a more accurate (closed form) 
+    solution as determined by the sum of currents at the diode node, which 
+    should be zero at electrical steady state.
+    Inputs to this function can include scalars and pandas.Series, but it
+    always outputs a float64 numpy.ndarray regardless of input type(s).
+
     Parameters
     ----------
-    V : float64 numpy.ndarray or scalar
+    V : numeric
         Device voltage [V]
-    IL : float64 numpy.ndarray or scalar
+
+    IL : numeric
         Device photocurrent [A]
-    I0 : float64 numpy.ndarray or scalar
+
+    I0 : numeric
         Device saturation current [A]
-    nNsVth : float64 numpy.ndarray or scalar
+
+    nNsVth : numeric
         Device thermal voltage [V]
-    Rs : float64 numpy.ndarray or scalar
+
+    Rs : numeric
         Device series resistance [Ohm]
-    Rsh : float64 numpy.ndarray or scalar
+
+    Rsh : numeric
         Device shunt resistance [Ohm]
+
     return_meta_dict : boolean scalar
         Return additional computation metadata dictionary
-    
+
     Returns
     -------
     I : float64 numpy.ndarray
         Device current [A]
-    meta_dict : python dictionary (optional, returned when return_meta_dict=True)
+
+    meta_dict : dictionary (optional, returned when return_meta_dict=True)
         Metadata for computation
-            meta_dict['current_sum_at_diode_node'] : float64 numpy.ndarray like I
-                Sum of currents at the diode node in equivalent circuit model [A]
+
+            meta_dict['sdm_current_sum'] : float64 numpy.ndarray like I
+                Sum of currents at diode node in equivalent circuit model [A]
+
             meta_dict['zero_Rs_idx'] : boolean numpy.ndarray like I
                 Indices where zero series resistance gives best solution
+
+    References
+    ----------
+    [1] A. Jain, A. Kapoor, "Exact analytical solutions of the parameters of 
+    real solar cells using Lambert W-function", Solar Energy Materials and 
+    Solar Cells, 81 (2004) 269-277.
     '''
 
-    # TODO Check if this qualifies as a numpy ufunc
-        
     # Ensure inputs are all np.ndarray with np.float64 type
     V = np.asarray(V, np.float64)
     IL = np.asarray(IL, np.float64)
@@ -2095,21 +2156,26 @@ def i_from_v_alt(V, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
     # Default computation of I using zero_term keeps Rs shape info
     zero_term = np.zeros_like(Rs)
     I = IL - I0*np.expm1((V + zero_term)/nNsVth) - (V + zero_term)/Rsh
-    current_sum_at_diode_node_out = current_sum_at_diode_node(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    sdm_current_sum_out = \
+sdm_current_sum(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
 
     # Computation of I using LambertW for provided Rs
     I_LambertW = i_from_v(Rsh, Rs, nNsVth, V, I0, IL)
-    current_sum_at_diode_node_LambertW = current_sum_at_diode_node(V=V, I=I_LambertW, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    sdm_current_sum_LambertW = \
+sdm_current_sum(V=V, I=I_LambertW, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
 
     # Compute selection indices (may be a scalar boolean)
-    nonzero_Rs_idx = np.logical_and(np.isfinite(current_sum_at_diode_node_LambertW), np.absolute(current_sum_at_diode_node_LambertW) <= np.absolute(current_sum_at_diode_node_out))
-    
+    nonzero_Rs_idx = np.logical_and(np.isfinite(sdm_current_sum_LambertW), \
+np.absolute(sdm_current_sum_LambertW) <= np.absolute(sdm_current_sum_out))
+
     # These are always np.ndarray
     I = np.where(nonzero_Rs_idx, I_LambertW, I)
-    current_sum_at_diode_node_out = np.where(nonzero_Rs_idx, current_sum_at_diode_node_LambertW, current_sum_at_diode_node_out)
-    
+    sdm_current_sum_out = \
+np.where(nonzero_Rs_idx, sdm_current_sum_LambertW, sdm_current_sum_out)
+
     if return_meta_dict:
-        return I, {'current_sum_at_diode_node' : current_sum_at_diode_node_out, 'zero_Rs_idx' : np.array(np.logical_not(nonzero_Rs_idx))}
+        return I, {'sdm_current_sum' : sdm_current_sum_out, \
+'zero_Rs_idx' : np.array(np.logical_not(nonzero_Rs_idx))}
     else:
         return I
 
diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index 39798d0958..9650d98cfb 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -362,7 +362,7 @@ def test_PVSystem_calcparams_desoto(cec_module_params):
      'nNsVth' : 0.5,
      'Rs' : 0.1,
      'Rsh' : 20.,
-     'current_sum_at_diode_node_expected' : np.float64(7. - 6.e-7*np.expm1((40. + 3.*0.1)/0.5) - (40. + 3.*0.1)/20. - 3.)
+     'sdm_current_sum_expected' : np.float64(7. - 6.e-7*np.expm1((40. + 3.*0.1)/0.5) - (40. + 3.*0.1)/20. - 3.)
     },
     { # Can handle all rank-0 array inputs
      'V' : np.array(40.),
@@ -372,7 +372,7 @@ def test_PVSystem_calcparams_desoto(cec_module_params):
      'nNsVth' : np.array(0.5),
      'Rs' : np.array(0.1),
      'Rsh' : np.array(20.),
-     'current_sum_at_diode_node_expected' : np.float64(7. - 6.e-7*np.expm1((40. + 3.*0.1)/0.5) - (40. + 3.*0.1)/20. - 3.)
+     'sdm_current_sum_expected' : np.float64(7. - 6.e-7*np.expm1((40. + 3.*0.1)/0.5) - (40. + 3.*0.1)/20. - 3.)
     },
     { # Can handle all rank-1 singleton array inputs
      'V' : np.array([40.]),
@@ -382,7 +382,7 @@ def test_PVSystem_calcparams_desoto(cec_module_params):
      'nNsVth' : np.array([0.5]),
      'Rs' : np.array([0.1]),
      'Rsh' : np.array([20.]),
-     'current_sum_at_diode_node_expected' : np.array([7. - 6.e-7*np.expm1((40. + 3.*0.1)/0.5) - (40. + 3.*0.1)/20. - 3.])
+     'sdm_current_sum_expected' : np.array([7. - 6.e-7*np.expm1((40. + 3.*0.1)/0.5) - (40. + 3.*0.1)/20. - 3.])
     },
     { # Can handle all rank-1 non-singleton array inputs
      'V' : np.array([40., 0. , 0.]),
@@ -392,7 +392,7 @@ def test_PVSystem_calcparams_desoto(cec_module_params):
      'nNsVth' : np.array([0.5, 0.5, 0.5]),
      'Rs' : np.array([0.1, 0.1, 0.1]),
      'Rsh' : np.array([20., 20., 20.]),
-     'current_sum_at_diode_node_expected' : np.array([7. - 6.e-7*np.expm1(40./0.5) - 40./0.1, 7., 7. - 6.e-7*np.expm1(3.*0.1/0.5) - 3.*0.1/20. - 3.])
+     'sdm_current_sum_expected' : np.array([7. - 6.e-7*np.expm1(40./0.5) - 40./0.1, 7., 7. - 6.e-7*np.expm1(3.*0.1/0.5) - 3.*0.1/20. - 3.])
     },
     { # Can handle mixed inputs with non-singleton Pandas Series
      'V' : pd.Series([40., 0. , 0.]),
@@ -402,7 +402,7 @@ def test_PVSystem_calcparams_desoto(cec_module_params):
      'nNsVth' : 0.5,
      'Rs' : 0.1,
      'Rsh' : 20.,
-     'current_sum_at_diode_node_expected' : pd.Series([7. - 6.e-7*np.expm1(40./0.5) - 40./0.1, 7., 7. - 6.e-7*np.expm1(3.*0.1/0.5) - 3.*0.1/20. - 3.])
+     'sdm_current_sum_expected' : pd.Series([7. - 6.e-7*np.expm1(40./0.5) - 40./0.1, 7., 7. - 6.e-7*np.expm1(3.*0.1/0.5) - 3.*0.1/20. - 3.])
     },
     { # Can handle mixed inputs with rank-2 arrays
      'V' : np.array([[40., 0. , 0.], [0., 0. , 40.]]),
@@ -412,8 +412,8 @@ def test_PVSystem_calcparams_desoto(cec_module_params):
      'nNsVth' : np.array(0.5),
      'Rs' : np.array([0.1]),
      'Rsh' : np.full((2,3), 20.),
-     'current_sum_at_diode_node_expected' : np.array([[7. - 6.e-7*np.expm1(40./0.5) - 40./0.1, 7., 7. - 6.e-7*np.expm1(3.*0.1/0.5) - 3.*0.1/20. - 3.], \
-                                                      [ 7. - 6.e-7*np.expm1(3.*0.1/0.5) - 3.*0.1/20. - 3., 7., 7. - 6.e-7*np.expm1(40./0.5) - 40./0.1]])
+     'sdm_current_sum_expected' : np.array([[7. - 6.e-7*np.expm1(40./0.5) - 40./0.1, 7., 7. - 6.e-7*np.expm1(3.*0.1/0.5) - 3.*0.1/20. - 3.], \
+                                            [ 7. - 6.e-7*np.expm1(3.*0.1/0.5) - 3.*0.1/20. - 3., 7., 7. - 6.e-7*np.expm1(40./0.5) - 40./0.1]])
     },
     { # Can handle infinite shunt resistance with positive series resistance
      'V' : 40.,
@@ -423,7 +423,7 @@ def test_PVSystem_calcparams_desoto(cec_module_params):
      'nNsVth' : 0.5,
      'Rs' : 0.1,
      'Rsh' : np.inf,
-     'current_sum_at_diode_node_expected' : np.float64(7. - 6.e-7*np.expm1((40. + 3.*0.1)/0.5) - 3.)
+     'sdm_current_sum_expected' : np.float64(7. - 6.e-7*np.expm1((40. + 3.*0.1)/0.5) - 3.)
     },
     { # Can handle infinite shunt resistance with zero series resistance
      'V' : 40.,
@@ -433,29 +433,29 @@ def test_PVSystem_calcparams_desoto(cec_module_params):
      'nNsVth' : 0.5,
      'Rs' : 0.,
      'Rsh' : np.inf,
-     'current_sum_at_diode_node_expected' : np.float64(7. - 6.e-7*np.expm1(40./0.5) - 3.)
+     'sdm_current_sum_expected' : np.float64(7. - 6.e-7*np.expm1(40./0.5) - 3.)
     }])
-def fixture_current_sum_at_diode_node(request):
+def fixture_sdm_current_sum(request):
     return request.param
 
-def test_current_sum_at_diode_node(fixture_current_sum_at_diode_node):
+def test_sdm_current_sum(fixture_sdm_current_sum):
     # Note: The computation of this function is so straight forward that we do
     #  NOT extensively verify ufunc behavior
     
     # Solution set loaded from fixture
-    V = fixture_current_sum_at_diode_node['V']
-    I = fixture_current_sum_at_diode_node['I']
-    IL = fixture_current_sum_at_diode_node['IL']
-    I0 = fixture_current_sum_at_diode_node['I0']
-    nNsVth = fixture_current_sum_at_diode_node['nNsVth']
-    Rs = fixture_current_sum_at_diode_node['Rs']
-    Rsh = fixture_current_sum_at_diode_node['Rsh']
-    current_sum_at_diode_node_expected = fixture_current_sum_at_diode_node['current_sum_at_diode_node_expected']
+    V = fixture_sdm_current_sum['V']
+    I = fixture_sdm_current_sum['I']
+    IL = fixture_sdm_current_sum['IL']
+    I0 = fixture_sdm_current_sum['I0']
+    nNsVth = fixture_sdm_current_sum['nNsVth']
+    Rs = fixture_sdm_current_sum['Rs']
+    Rsh = fixture_sdm_current_sum['Rsh']
+    sdm_current_sum_expected = fixture_sdm_current_sum['sdm_current_sum_expected']
 
-    current_sum_at_diode_node = pvsystem.current_sum_at_diode_node(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
-    assert(isinstance(current_sum_at_diode_node, type(current_sum_at_diode_node_expected)))
-    assert(isinstance(current_sum_at_diode_node.dtype, type(current_sum_at_diode_node_expected.dtype)))
-    assert_array_equal(current_sum_at_diode_node, current_sum_at_diode_node_expected)
+    sdm_current_sum = pvsystem.sdm_current_sum(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    assert(isinstance(sdm_current_sum, type(sdm_current_sum_expected)))
+    assert(isinstance(sdm_current_sum.dtype, type(sdm_current_sum_expected.dtype)))
+    assert_array_equal(sdm_current_sum, sdm_current_sum_expected)
 
 
 @pytest.fixture(params=[
@@ -467,7 +467,7 @@ def test_current_sum_at_diode_node(fixture_current_sum_at_diode_node):
      'I0' : 6.e-7,
      'IL' : 7.,
      'V_expected' : np.array(7.5049875193450521),
-     'current_sum_at_diode_node_expected' : np.array(0.),
+     'sdm_current_sum_expected' : np.array(0.),
      'inf_Rsh_idx_expected' : np.array(False)
     },
     { # Can handle all rank-0 array inputs
@@ -478,7 +478,7 @@ def test_current_sum_at_diode_node(fixture_current_sum_at_diode_node):
      'I0' : np.array(6.e-7),
      'IL' : np.array(7.),
      'V_expected' : np.array(7.5049875193450521),
-     'current_sum_at_diode_node_expected' : np.array(0.),
+     'sdm_current_sum_expected' : np.array(0.),
      'inf_Rsh_idx_expected' : np.array(False)
     },
     { # Can handle all rank-1 singleton array inputs
@@ -489,7 +489,7 @@ def test_current_sum_at_diode_node(fixture_current_sum_at_diode_node):
      'I0' : np.array([6.e-7]),
      'IL' : np.array([7.]),
      'V_expected' : np.array([7.5049875193450521]),
-     'current_sum_at_diode_node_expected' : np.array([0.]),
+     'sdm_current_sum_expected' : np.array([0.]),
      'inf_Rsh_idx_expected' : np.array([False])
     },
     { # Can handle all rank-1 non-singleton array inputs with infinite shunt
@@ -502,7 +502,7 @@ def test_current_sum_at_diode_node(fixture_current_sum_at_diode_node):
      'I0' : np.array([6.e-7, 6.e-7]),
      'IL' : np.array([7., 7.]),
      'V_expected' : np.array([0.5*(np.log(7. + 6.e-7) - np.log(6.e-7)), 7.5049875193450521]),
-     'current_sum_at_diode_node_expected' : np.array([0., 0.]),
+     'sdm_current_sum_expected' : np.array([0., 0.]),
      'inf_Rsh_idx_expected' : np.array([True, False])
     },
     { # Can handle mixed inputs with a rank-2 array with infinite shunt
@@ -515,7 +515,7 @@ def test_current_sum_at_diode_node(fixture_current_sum_at_diode_node):
      'I0' : np.array([6.e-7]),
      'IL' : np.array([7.]),
      'V_expected' : 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))*np.ones((2,2)),
-     'current_sum_at_diode_node_expected' : np.array([[0., 0.], [0., 0.]]),
+     'sdm_current_sum_expected' : np.array([[0., 0.], [0., 0.]]),
      'inf_Rsh_idx_expected' : np.array([[True, True], [True, True]])
     },
     { # Can handle ideal series and shunt, Rsh=inf and Rs=0 give 
@@ -527,7 +527,7 @@ def test_current_sum_at_diode_node(fixture_current_sum_at_diode_node):
      'I0' : 6.e-7,
      'IL' : 7.,
      'V_expected' : np.array([0., 0.5*(np.log(7. - 7./2. + 6.e-7) - np.log(6.e-7)), 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))]),
-     'current_sum_at_diode_node_expected' : np.array([0., 0., 0.]),
+     'sdm_current_sum_expected' : np.array([0., 0., 0.]),
      'inf_Rsh_idx_expected' : np.array([True, True, True])
     },
     { # Can handle only ideal series resistance, no closed form solution
@@ -538,7 +538,7 @@ def test_current_sum_at_diode_node(fixture_current_sum_at_diode_node):
      'I0' : 6.e-7,
      'IL' : 7.,
      'V_expected' : np.array(7.804987519345062),
-     'current_sum_at_diode_node_expected' : np.array(0.),
+     'sdm_current_sum_expected' : np.array(0.),
      'inf_Rsh_idx_expected' : np.array(False)
     },
     { # Can handle all python scalar inputs with big LambertW arg
@@ -549,7 +549,7 @@ def test_current_sum_at_diode_node(fixture_current_sum_at_diode_node):
      'I0' : 6.e-10,
      'IL' : 1.2,
      'V_expected' : np.array(86.320000493521079),
-     'current_sum_at_diode_node_expected' : np.array(0.),
+     'sdm_current_sum_expected' : np.array(0.),
      'inf_Rsh_idx_expected' : np.array(False)
     },
     { # Can handle all python scalar inputs with bigger LambertW arg
@@ -562,41 +562,41 @@ def test_current_sum_at_diode_node(fixture_current_sum_at_diode_node):
      'I0' : 7.05196029e-08,
      'IL' : 10.491262,
      'V_expected' : np.array(54.303958833791455),
-     'current_sum_at_diode_node_expected' : np.array(0.),
+     'sdm_current_sum_expected' : np.array(0.),
      'inf_Rsh_idx_expected' : np.array(False)
     }])
-def fixture_v_from_i_alt(request):
+def fixture_sdm_v_from_i(request):
     return request.param
 
 @requires_scipy
-def test_v_from_i_alt(fixture_v_from_i_alt):
+def test_sdm_v_from_i(fixture_sdm_v_from_i):
     # Solution set loaded from fixture
-    Rsh = fixture_v_from_i_alt['Rsh']
-    Rs = fixture_v_from_i_alt['Rs']
-    nNsVth = fixture_v_from_i_alt['nNsVth']
-    I = fixture_v_from_i_alt['I']
-    I0 = fixture_v_from_i_alt['I0']
-    IL = fixture_v_from_i_alt['IL']
-    V_expected = fixture_v_from_i_alt['V_expected']
-    current_sum_at_diode_node_expected = fixture_v_from_i_alt['current_sum_at_diode_node_expected']
-    inf_Rsh_idx_expected = fixture_v_from_i_alt['inf_Rsh_idx_expected']
+    Rsh = fixture_sdm_v_from_i['Rsh']
+    Rs = fixture_sdm_v_from_i['Rs']
+    nNsVth = fixture_sdm_v_from_i['nNsVth']
+    I = fixture_sdm_v_from_i['I']
+    I0 = fixture_sdm_v_from_i['I0']
+    IL = fixture_sdm_v_from_i['IL']
+    V_expected = fixture_sdm_v_from_i['V_expected']
+    sdm_current_sum_expected = fixture_sdm_v_from_i['sdm_current_sum_expected']
+    inf_Rsh_idx_expected = fixture_sdm_v_from_i['inf_Rsh_idx_expected']
     
     # Convergence criteria
     atol = 1.e-11
     
-    V = pvsystem.v_from_i_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, I=I, I0=I0, IL=IL)
+    V = pvsystem.sdm_v_from_i(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, I=I, I0=I0, IL=IL)
     assert(isinstance(V, type(V_expected)))
     assert(isinstance(V.dtype, type(V_expected.dtype)))
     assert_allclose(V, V_expected, atol=atol)
-    _, meta_dict = pvsystem.v_from_i_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, I=I, I0=I0, IL=IL, return_meta_dict=True)
-    assert(isinstance(meta_dict['current_sum_at_diode_node'], type(V)))
-    assert(isinstance(meta_dict['current_sum_at_diode_node'].dtype, type(V.dtype)))
-    assert_allclose(meta_dict['current_sum_at_diode_node'], current_sum_at_diode_node_expected, atol=atol)
+    _, meta_dict = pvsystem.sdm_v_from_i(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, I=I, I0=I0, IL=IL, return_meta_dict=True)
+    assert(isinstance(meta_dict['sdm_current_sum'], type(V)))
+    assert(isinstance(meta_dict['sdm_current_sum'].dtype, type(V.dtype)))
+    assert_allclose(meta_dict['sdm_current_sum'], sdm_current_sum_expected, atol=atol)
     assert(isinstance(meta_dict['inf_Rsh_idx'], type(V)))
     assert(isinstance(meta_dict['inf_Rsh_idx'].dtype, type(V.dtype)))
     assert_array_equal(meta_dict['inf_Rsh_idx'], inf_Rsh_idx_expected)
     
-    # TODO Stability as Rs->0^+ and/or Rsh->inf
+    # TODO Stability as Rs->0^+ and/or Rsh->inf and benchmarks
 
 
 @pytest.fixture(params=[
@@ -608,7 +608,7 @@ def test_v_from_i_alt(fixture_v_from_i_alt):
      'I0' : 6.e-7,
      'IL' : 7.,
      'I_expected' : np.array(-299.746389916),
-     'current_sum_at_diode_node_expected' : np.array(0.),
+     'sdm_current_sum_expected' : np.array(0.),
      'zero_Rs_idx_expected' : np.array(False)
     },
     { # Can handle all rank-0 array inputs
@@ -619,7 +619,7 @@ def test_v_from_i_alt(fixture_v_from_i_alt):
      'I0' : np.array(6.e-7),
      'IL' : np.array(7.),
      'I_expected' : np.array(-299.746389916),
-     'current_sum_at_diode_node_expected' : np.array(0.),
+     'sdm_current_sum_expected' : np.array(0.),
      'zero_Rs_idx_expected' : np.array(False)
     },
     { # Can handle all rank-1 singleton array inputs
@@ -630,7 +630,7 @@ def test_v_from_i_alt(fixture_v_from_i_alt):
      'I0' : np.array([6.e-7]),
      'IL' : np.array([7.]),
      'I_expected' : np.array([-299.746389916]),
-     'current_sum_at_diode_node_expected' : np.array([0.]),
+     'sdm_current_sum_expected' : np.array([0.]),
      'zero_Rs_idx_expected' : np.array([False])
     },
     { # Can handle all rank-1 non-singleton array inputs with a zero 
@@ -642,7 +642,7 @@ def test_v_from_i_alt(fixture_v_from_i_alt):
      'I0' : np.array([6.e-7, 6.e-7]),
      'IL' : np.array([7., 7.]),
      'I_expected' : np.array([7., -299.746389916]),
-     'current_sum_at_diode_node_expected' : np.array([0., 0.]),
+     'sdm_current_sum_expected' : np.array([0., 0.]),
      'zero_Rs_idx_expected' : np.array([True, False])
     },
     { # Can handle mixed inputs with a rank-2 array with zero series 
@@ -654,7 +654,7 @@ def test_v_from_i_alt(fixture_v_from_i_alt):
      'I0' : np.array([6.e-7]),
      'IL' : np.array([7.]),
      'I_expected' : np.array([[7., 7.], [7., 7.]]),
-     'current_sum_at_diode_node_expected' : np.array([[0., 0.], [0., 0.]]),
+     'sdm_current_sum_expected' : np.array([[0., 0.], [0., 0.]]),
      'zero_Rs_idx_expected' : np.array([[True, True], [True, True]])
     },
     { # Can handle ideal series and shunt, Rsh=inf and Rs=0 give 
@@ -666,7 +666,7 @@ def test_v_from_i_alt(fixture_v_from_i_alt):
      'I0' : 6.e-7,
      'IL' : 7.,
      'I_expected' : np.array([7., 7. - 6.e-7*np.expm1((np.log(7. + 6.e-7) - np.log(6.e-7))/2.), 0.]),
-     'current_sum_at_diode_node_expected' : np.array([0., 0., 0.]),
+     'sdm_current_sum_expected' : np.array([0., 0., 0.]),
      'zero_Rs_idx_expected' : np.array([True, True, True])
     },
     { # Can handle only ideal shunt resistance, no closed form solution
@@ -677,41 +677,41 @@ def test_v_from_i_alt(fixture_v_from_i_alt):
      'I0' : 6.e-7,
      'IL' : 7.,
      'I_expected' : np.array(-299.7383436645412),
-     'current_sum_at_diode_node_expected' : np.array(0.),
+     'sdm_current_sum_expected' : np.array(0.),
      'zero_Rs_idx_expected' : np.array(False)
     }])
-def fixture_i_from_v_alt(request):
+def fixture_sdm_i_from_v(request):
     return request.param
 
 @requires_scipy
-def test_i_from_v_alt(fixture_i_from_v_alt):
+def test_sdm_i_from_v(fixture_sdm_i_from_v):
     # Solution set loaded from fixture
-    Rsh = fixture_i_from_v_alt['Rsh']
-    Rs = fixture_i_from_v_alt['Rs']
-    nNsVth = fixture_i_from_v_alt['nNsVth']
-    V = fixture_i_from_v_alt['V']
-    I0 = fixture_i_from_v_alt['I0']
-    IL = fixture_i_from_v_alt['IL']
-    I_expected = fixture_i_from_v_alt['I_expected']
-    current_sum_at_diode_node_expected = fixture_i_from_v_alt['current_sum_at_diode_node_expected']
-    zero_Rs_idx_expected = fixture_i_from_v_alt['zero_Rs_idx_expected']
+    Rsh = fixture_sdm_i_from_v['Rsh']
+    Rs = fixture_sdm_i_from_v['Rs']
+    nNsVth = fixture_sdm_i_from_v['nNsVth']
+    V = fixture_sdm_i_from_v['V']
+    I0 = fixture_sdm_i_from_v['I0']
+    IL = fixture_sdm_i_from_v['IL']
+    I_expected = fixture_sdm_i_from_v['I_expected']
+    sdm_current_sum_expected = fixture_sdm_i_from_v['sdm_current_sum_expected']
+    zero_Rs_idx_expected = fixture_sdm_i_from_v['zero_Rs_idx_expected']
 
     # Convergence criteria
     atol = 1.e-11
     
-    I = pvsystem.i_from_v_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL)
+    I = pvsystem.sdm_i_from_v(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL)
     assert(isinstance(I, type(I_expected)))
     assert(isinstance(I.dtype, type(I_expected.dtype)))
     assert_allclose(I, I_expected, atol=atol)
-    _, meta_dict = pvsystem.i_from_v_alt(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL, return_meta_dict=True)
-    assert(isinstance(meta_dict['current_sum_at_diode_node'], type(I)))
-    assert(isinstance(meta_dict['current_sum_at_diode_node'].dtype, type(I.dtype)))
-    assert_allclose(meta_dict['current_sum_at_diode_node'], current_sum_at_diode_node_expected, atol=atol)
+    _, meta_dict = pvsystem.sdm_i_from_v(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL, return_meta_dict=True)
+    assert(isinstance(meta_dict['sdm_current_sum'], type(I)))
+    assert(isinstance(meta_dict['sdm_current_sum'].dtype, type(I.dtype)))
+    assert_allclose(meta_dict['sdm_current_sum'], sdm_current_sum_expected, atol=atol)
     assert(isinstance(meta_dict['zero_Rs_idx'], type(I)))
     assert(isinstance(meta_dict['zero_Rs_idx'].dtype, type(I.dtype)))
     assert_array_equal(meta_dict['zero_Rs_idx'], zero_Rs_idx_expected)
     
-    # TODO Stability as Rs->0^+ and/or Rsh->inf
+    # TODO Stability as Rs->0^+ and/or Rsh->inf and benchmarks
 
 
 @requires_scipy

From c9129dc56588d9aacea8212c22d8b6c8f1554407 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Wed, 26 Jul 2017 20:45:23 -0600
Subject: [PATCH 09/22] Deprecate replaced functions and flake8

---
 pvlib/pvsystem.py           | 78 +++++++++++++++++++------------------
 pvlib/test/test_pvsystem.py |  2 +-
 2 files changed, 41 insertions(+), 39 deletions(-)

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 148e9c92a5..065696e6a4 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -437,18 +437,18 @@ def singlediode(self, photocurrent, saturation_current,
 
     def i_from_v(self, resistance_shunt, resistance_series, nNsVth, voltage,
                  saturation_current, photocurrent):
-        """Wrapper around the :py:func:`i_from_v` function.
+        """Wrapper around the :py:func:`sdm_i_from_v` function.
 
         Parameters
         ----------
-        See pvsystem.i_from_v for details
+        See pvsystem.sdm_i_from_v for details
 
         Returns
         -------
-        See pvsystem.i_from_v for details
+        See pvsystem.sdm_i_from_v for details
         """
-        return i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
-                        saturation_current, photocurrent)
+        return sdm_i_from_v(voltage, photocurrent, saturation_current,
+                            nNsVth, resistance_series, resistance_shunt)
 
     # inverter now specified by self.inverter_parameters
     def snlinverter(self, v_dc, p_dc):
@@ -1790,12 +1790,14 @@ def _pwr_optfcn(df, loc):
 
     I = i_from_v(df['r_sh'], df['r_s'], df['nNsVth'],
                  df[loc], df['i_0'], df['i_l'])
-    return I*df[loc]
+    return I * df[loc]
 
 
 def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
              saturation_current, photocurrent):
     '''
+    DEPRECATED: Use sdm_v_from_i() instead.
+
     Calculates voltage from current per Eq 3 Jain and Kapoor 2004 [1].
 
     Parameters
@@ -1884,6 +1886,8 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
 def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
              saturation_current, photocurrent):
     '''
+    DEPRECATED: Use sdm_i_from_v() instead.
+
     Calculates current from voltage per Eq 2 Jain and Kapoor 2004 [1].
 
     Parameters
@@ -1996,8 +2000,8 @@ def sdm_current_sum(V, I, IL, I0, nNsVth, Rs, Rsh):
     Solar Cells, 81 (2004) 269-277.
     '''
 
-    # sdm_current_sum
-    return IL - I0*np.expm1((V + I*Rs)/nNsVth) - (V + I*Rs)/Rsh - I
+    VplusItimesRs = V + I * Rs
+    return IL - I0 * np.expm1(VplusItimesRs / nNsVth) - VplusItimesRs / Rsh - I
 
 
 def sdm_v_from_i(I, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
@@ -2031,7 +2035,7 @@ def sdm_v_from_i(I, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
     Rsh : numeric
         Device shunt resistance [Ohm]
 
-    return_meta_dict : boolean scalar
+    return_meta_dict : boolean scalar (default is False)
         Return additional computation metadata dictionary
 
     Returns
@@ -2063,29 +2067,28 @@ def sdm_v_from_i(I, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
     Rs = np.asarray(Rs, np.float64)
     Rsh = np.asarray(Rsh, np.float64)
 
-    # Default computation of V using zero_term keeps Rsh shape info
-    zero_term = np.zeros_like(I*Rs/Rsh)
-    V = nNsVth*(np.log(IL - I - zero_term + I0) - np.log(I0)) - I*Rs
-    sdm_current_sum_out = \
-sdm_current_sum(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    # Default computation of V using np.zeros_like keeps Rsh shape info
+    V = nNsVth * (np.log(IL - I - np.zeros_like(I * Rs / Rsh) + I0) -
+                  np.log(I0)) - I * Rs
+    sdm_current_sum_out = sdm_current_sum(V, I, IL, I0, nNsVth, Rs, Rsh)
 
     # Computation of V using LambertW for provided Rsh
-    V_LambertW = v_from_i(Rsh, Rs, nNsVth, I, I0, IL)
-    sdm_current_sum_LambertW = \
-sdm_current_sum(V=V_LambertW, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    V_lw = v_from_i(Rsh, Rs, nNsVth, I, I0, IL)
+    sdm_current_sum_lw = sdm_current_sum(V_lw, I, IL, I0, nNsVth, Rs, Rsh)
 
     # Compute selection indices (may be a scalar boolean)
-    finite_Rsh_idx = np.logical_and(np.isfinite(sdm_current_sum_LambertW), \
-np.absolute(sdm_current_sum_LambertW) <= np.absolute(sdm_current_sum_out))
+    finite_Rsh_idx = np.logical_and(np.isfinite(sdm_current_sum_lw),
+                                    np.absolute(sdm_current_sum_lw) <=
+                                    np.absolute(sdm_current_sum_out))
 
     # These are always np.ndarray
-    V = np.where(finite_Rsh_idx, V_LambertW, V)
-    sdm_current_sum_out = \
-np.where(finite_Rsh_idx, sdm_current_sum_LambertW, sdm_current_sum_out)
+    V = np.where(finite_Rsh_idx, V_lw, V)
+    sdm_current_sum_out = np.where(finite_Rsh_idx, sdm_current_sum_lw,
+                                   sdm_current_sum_out)
 
     if return_meta_dict:
-        return V, {'sdm_current_sum' : sdm_current_sum_out, \
-'inf_Rsh_idx' : np.array(np.logical_not(finite_Rsh_idx))}
+        return V, {'sdm_current_sum': sdm_current_sum_out,
+                   'inf_Rsh_idx': np.array(np.logical_not(finite_Rsh_idx))}
     else:
         return V
 
@@ -2121,7 +2124,7 @@ def sdm_i_from_v(V, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
     Rsh : numeric
         Device shunt resistance [Ohm]
 
-    return_meta_dict : boolean scalar
+    return_meta_dict : boolean scalar (default is False)
         Return additional computation metadata dictionary
 
     Returns
@@ -2155,27 +2158,26 @@ def sdm_i_from_v(V, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
 
     # Default computation of I using zero_term keeps Rs shape info
     zero_term = np.zeros_like(Rs)
-    I = IL - I0*np.expm1((V + zero_term)/nNsVth) - (V + zero_term)/Rsh
-    sdm_current_sum_out = \
-sdm_current_sum(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    I = IL - I0 * np.expm1((V + zero_term) / nNsVth) - (V + zero_term) / Rsh
+    sdm_current_sum_out = sdm_current_sum(V, I, IL, I0, nNsVth, Rs, Rsh)
 
     # Computation of I using LambertW for provided Rs
-    I_LambertW = i_from_v(Rsh, Rs, nNsVth, V, I0, IL)
-    sdm_current_sum_LambertW = \
-sdm_current_sum(V=V, I=I_LambertW, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
+    I_lw = i_from_v(Rsh, Rs, nNsVth, V, I0, IL)
+    sdm_current_sum_lw = sdm_current_sum(V, I_lw, IL, I0, nNsVth, Rs, Rsh)
 
     # Compute selection indices (may be a scalar boolean)
-    nonzero_Rs_idx = np.logical_and(np.isfinite(sdm_current_sum_LambertW), \
-np.absolute(sdm_current_sum_LambertW) <= np.absolute(sdm_current_sum_out))
+    nonzero_Rs_idx = np.logical_and(np.isfinite(sdm_current_sum_lw),
+                                    np.absolute(sdm_current_sum_lw) <=
+                                    np.absolute(sdm_current_sum_out))
 
     # These are always np.ndarray
-    I = np.where(nonzero_Rs_idx, I_LambertW, I)
-    sdm_current_sum_out = \
-np.where(nonzero_Rs_idx, sdm_current_sum_LambertW, sdm_current_sum_out)
+    I = np.where(nonzero_Rs_idx, I_lw, I)
+    sdm_current_sum_out = np.where(nonzero_Rs_idx, sdm_current_sum_lw,
+                                   sdm_current_sum_out)
 
     if return_meta_dict:
-        return I, {'sdm_current_sum' : sdm_current_sum_out, \
-'zero_Rs_idx' : np.array(np.logical_not(nonzero_Rs_idx))}
+        return I, {'sdm_current_sum': sdm_current_sum_out,
+                   'zero_Rs_idx': np.array(np.logical_not(nonzero_Rs_idx))}
     else:
         return I
 
diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index 9650d98cfb..b1f2c3f74e 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -439,7 +439,7 @@ def fixture_sdm_current_sum(request):
     return request.param
 
 def test_sdm_current_sum(fixture_sdm_current_sum):
-    # Note: The computation of this function is so straight forward that we do
+    # Note: The computation of this function is so straightforward that we do
     #  NOT extensively verify ufunc behavior
     
     # Solution set loaded from fixture

From 85e957d6bb3d34846ad4001f6ed1f16d63eec613 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Wed, 26 Jul 2017 22:29:58 -0600
Subject: [PATCH 10/22] Add release documentation and flake8 again

---
 docs/sphinx/source/api.rst             |  9 ++++++---
 docs/sphinx/source/whatsnew/v0.4.6.rst | 16 ++++++++++++++++
 pvlib/atmosphere.py                    |  2 +-
 pvlib/pvsystem.py                      |  9 +++++----
 pvlib/test/test_irradiance.py          |  2 +-
 pvlib/test/test_modelchain.py          |  1 +
 pvlib/test/test_pvsystem.py            |  4 ----
 7 files changed, 30 insertions(+), 13 deletions(-)

diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst
index 5fd65b8733..d1c54fd97e 100644
--- a/docs/sphinx/source/api.rst
+++ b/docs/sphinx/source/api.rst
@@ -200,10 +200,13 @@ Functions relevant for the single diode model.
 .. autosummary::
    :toctree: generated/
 
-   pvsystem.singlediode
    pvsystem.calcparams_desoto
-   pvsystem.v_from_i
-   pvsystem.i_from_v
+   pvsystem.i_from_v (DEPRECATED: Use pvsystem.sdm_i_from_v instead)
+   pvsystem.sdm_current_sum
+   pvsystem.sdm_i_from_v
+   pvsystem.sdm_v_from_i
+   pvsystem.singlediode
+   pvsystem.v_from_i (DEPRECATED: Use pvsystem.sdm_v_from_i instead)
 
 SAPM model
 ----------
diff --git a/docs/sphinx/source/whatsnew/v0.4.6.rst b/docs/sphinx/source/whatsnew/v0.4.6.rst
index 37f968d851..b1b3bbdaa1 100644
--- a/docs/sphinx/source/whatsnew/v0.4.6.rst
+++ b/docs/sphinx/source/whatsnew/v0.4.6.rst
@@ -17,10 +17,22 @@ Bug fixes
 Enhancements
 ~~~~~~~~~~~~
 * Added default values to docstrings of all functions (:issue:`336`)
+* Ideal devices supported in single diode model, e.g., 
+  resistance_series = 0 and/or resistance_shunt = numpy.inf (:issue:`340`)
+* Computations for very near ideal devices are more numerically stable in 
+  single diode model (:issue:`340`)
 
 API Changes
 ~~~~~~~~~~~
 * Removed parameter w from _calc_d (:issue:`344`)
+* Deprecated `pvsystem.v_from_i` and `pvsystem.i_from_v` functions, which will 
+  be removed in a future release. Instead use `pvsystem.sdm_v_from_i` and
+  `pvsystem.sdm_i_from_v`, which are ~2x slower but can handle ideal devices. 
+  Note that these functions always return a numpy.ndarray (:issue:`340`)
+* Added `pvsystem.sdm_current_sum` function, which computes the sum of currents 
+  at the diode node in the single diode model and is used to improve numerical 
+  stability in `pvsystem.sdm_v_from_i` and `pvsystem.sdm_i_from_v` for very 
+  near ideal devices (:issue:`340`)
 
 Documentation
 ~~~~~~~~~~~~~
@@ -31,6 +43,9 @@ Testing
 
 * Added explicit tests for aoi and aoi_projection functions.
 * Update test of `ModelChain.__repr__` to take in account :issue:`352`
+* Significant new test cases added for `pvsystem.sdm_v_from_i` and
+  `pvsystem.sdm_i_from_v`, beyond but including cases for deprecated 
+  `pvsystem.v_from_i` and `pvsystem.i_from_v` (:issue:`340`)
 
 
 Contributors
@@ -41,3 +56,4 @@ Contributors
 * Alaina Kafkes
 * Birgit Schachler
 * Jonathan Gaffiot
+* Mark Campanelli
diff --git a/pvlib/atmosphere.py b/pvlib/atmosphere.py
index a69541c93a..4878cc796d 100644
--- a/pvlib/atmosphere.py
+++ b/pvlib/atmosphere.py
@@ -380,7 +380,7 @@ def first_solar_spectral_correction(pw, airmass_absolute, module_type=None,
         The module used to calculate the spectral correction
         coefficients corresponds to the Mult-crystalline silicon
         Manufacturer 2 Model C from [3]_. Spectral Response (SR) of CIGS
-        and a-Si modules used to derive coefficients can be found in [4]_ 
+        and a-Si modules used to derive coefficients can be found in [4]_
 
     coefficients : None or array-like, default None
         allows for entry of user defined spectral correction
diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index b82638f1a3..22477cde50 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -588,10 +588,9 @@ def __init__(self, pvsystem=None, location=None, **kwargs):
         Location.__init__(self, **new_kwargs)
 
     def __repr__(self):
-        attrs = [
-            'name', 'latitude', 'longitude', 'altitude', 'tz', 'surface_tilt',
-            'surface_azimuth', 'module', 'inverter', 'albedo', 'racking_model'
-                 ]
+        attrs = ['name', 'latitude', 'longitude', 'altitude', 'tz',
+                 'surface_tilt', 'surface_azimuth', 'module', 'inverter',
+                 'albedo', 'racking_model']
         return ('LocalizedPVSystem: \n  ' + '\n  '.join(
             ('{}: {}'.format(attr, getattr(self, attr)) for attr in attrs)))
 
@@ -2076,6 +2075,7 @@ def sdm_v_from_i(I, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
     sdm_current_sum_out = sdm_current_sum(V, I, IL, I0, nNsVth, Rs, Rsh)
 
     # Computation of V using LambertW for provided Rsh
+    # v_from_i is deprecated: move v_from_i code here when/if it is removed.
     V_lw = v_from_i(Rsh, Rs, nNsVth, I, I0, IL)
     sdm_current_sum_lw = sdm_current_sum(V_lw, I, IL, I0, nNsVth, Rs, Rsh)
 
@@ -2165,6 +2165,7 @@ def sdm_i_from_v(V, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
     sdm_current_sum_out = sdm_current_sum(V, I, IL, I0, nNsVth, Rs, Rsh)
 
     # Computation of I using LambertW for provided Rs
+    # i_from_v is deprecated: move i_from_v code here when/if it is removed.
     I_lw = i_from_v(Rsh, Rs, nNsVth, V, I0, IL)
     sdm_current_sum_lw = sdm_current_sum(V, I_lw, IL, I0, nNsVth, Rs, Rsh)
 
diff --git a/pvlib/test/test_irradiance.py b/pvlib/test/test_irradiance.py
index ba4a1ba74c..8c87f29eea 100755
--- a/pvlib/test/test_irradiance.py
+++ b/pvlib/test/test_irradiance.py
@@ -432,7 +432,7 @@ def test_dni():
     (90, 0, 30, 60, 75.5224878, 0.25),
     (90, 0, 30, 170, 119.4987042, -0.4924038)])
 def test_aoi_and_aoi_projection(surface_tilt, surface_azimuth, solar_zenith,
-                                solar_azimuth,  aoi_expected,
+                                solar_azimuth, aoi_expected,
                                 aoi_proj_expected):
     aoi = irradiance.aoi(surface_tilt, surface_azimuth, solar_zenith,
                          solar_azimuth)
diff --git a/pvlib/test/test_modelchain.py b/pvlib/test/test_modelchain.py
index e2199a283f..e3ab0ef61d 100644
--- a/pvlib/test/test_modelchain.py
+++ b/pvlib/test/test_modelchain.py
@@ -443,6 +443,7 @@ def test_ModelChain___repr__(system, location, strategy, strategy_str):
 
     assert mc.__repr__() == expected
 
+
 @requires_scipy
 def test_weather_irradiance_input(system, location):
     """Test will raise a warning and should be removed in future versions."""
diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index 5d446ec3a7..31bfb0122b 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -595,8 +595,6 @@ def test_sdm_v_from_i(fixture_sdm_v_from_i):
     assert(isinstance(meta_dict['inf_Rsh_idx'], type(V)))
     assert(isinstance(meta_dict['inf_Rsh_idx'].dtype, type(V.dtype)))
     assert_array_equal(meta_dict['inf_Rsh_idx'], inf_Rsh_idx_expected)
-    
-    # TODO Stability as Rs->0^+ and/or Rsh->inf and benchmarks
 
 
 @pytest.fixture(params=[
@@ -710,8 +708,6 @@ def test_sdm_i_from_v(fixture_sdm_i_from_v):
     assert(isinstance(meta_dict['zero_Rs_idx'], type(I)))
     assert(isinstance(meta_dict['zero_Rs_idx'].dtype, type(I.dtype)))
     assert_array_equal(meta_dict['zero_Rs_idx'], zero_Rs_idx_expected)
-    
-    # TODO Stability as Rs->0^+ and/or Rsh->inf and benchmarks
 
 
 @requires_scipy

From fb22171eb96df492db48070a1d9c43da7383f198 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Wed, 26 Jul 2017 23:07:33 -0600
Subject: [PATCH 11/22] Replace deprecated function usages and update
 test_singlediode_series_ivcurve

---
 pvlib/pvsystem.py           | 31 +++++++++++++-------------
 pvlib/test/test_pvsystem.py | 44 ++++++++++++++++++-------------------
 2 files changed, 36 insertions(+), 39 deletions(-)

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 22477cde50..5e3ec3194c 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1661,12 +1661,12 @@ def singlediode(photocurrent, saturation_current, resistance_series,
     '''
 
     # Find short circuit current using Lambert W
-    i_sc = i_from_v(resistance_shunt, resistance_series, nNsVth, 0.01,
-                    saturation_current, photocurrent)
+    i_sc = sdm_i_from_v(0., photocurrent, saturation_current, nNsVth,
+                        resistance_series, resistance_shunt)
 
     # Find open circuit voltage using Lambert W
-    v_oc = v_from_i(resistance_shunt, resistance_series, nNsVth, 0.0,
-                    saturation_current, photocurrent)
+    v_oc = sdm_v_from_i(0., photocurrent, saturation_current, nNsVth,
+                        resistance_series, resistance_shunt)
 
     params = {'r_sh': resistance_shunt,
               'r_s': resistance_series,
@@ -1674,19 +1674,19 @@ def singlediode(photocurrent, saturation_current, resistance_series,
               'i_0': saturation_current,
               'i_l': photocurrent}
 
-    p_mp, v_mp = _golden_sect_DataFrame(params, 0, v_oc*1.14, _pwr_optfcn)
+    p_mp, v_mp = _golden_sect_DataFrame(params, 0., v_oc * 1.14, _pwr_optfcn)
 
     # Invert the Power-Current curve. Find the current where the inverted power
     # is minimized. This is i_mp. Start the optimization at v_oc/2
-    i_mp = i_from_v(resistance_shunt, resistance_series, nNsVth, v_mp,
-                    saturation_current, photocurrent)
+    i_mp = sdm_i_from_v(v_mp, photocurrent, saturation_current, nNsVth,
+                        resistance_series, resistance_shunt)
 
     # Find Ix and Ixx using Lambert W
-    i_x = i_from_v(resistance_shunt, resistance_series, nNsVth,
-                   0.5*v_oc, saturation_current, photocurrent)
+    i_x = sdm_i_from_v(0.5 * v_oc, photocurrent, saturation_current, nNsVth,
+                       resistance_series, resistance_shunt)
 
-    i_xx = i_from_v(resistance_shunt, resistance_series, nNsVth,
-                    0.5*(v_oc+v_mp), saturation_current, photocurrent)
+    i_xx = sdm_i_from_v(0.5 * (v_oc + v_mp), photocurrent, saturation_current,
+                        nNsVth, resistance_series, resistance_shunt)
 
     out = OrderedDict()
     out['i_sc'] = i_sc
@@ -1701,9 +1701,8 @@ def singlediode(photocurrent, saturation_current, resistance_series,
     if ivcurve_pnts:
         ivcurve_v = (np.asarray(v_oc)[..., np.newaxis] *
                      np.linspace(0, 1, ivcurve_pnts))
-        ivcurve_i = i_from_v(
-            resistance_shunt, resistance_series, nNsVth, ivcurve_v.T,
-            saturation_current, photocurrent).T
+        ivcurve_i = sdm_i_from_v(ivcurve_v.T, photocurrent, saturation_current,
+                                 nNsVth, resistance_series, resistance_shunt).T
         out['v'] = ivcurve_v
         out['i'] = ivcurve_i
 
@@ -1790,8 +1789,8 @@ def _pwr_optfcn(df, loc):
     Function to find power from ``i_from_v``.
     '''
 
-    I = i_from_v(df['r_sh'], df['r_s'], df['nNsVth'],
-                 df[loc], df['i_0'], df['i_l'])
+    I = sdm_i_from_v(df[loc], df['i_l'], df['i_0'], df['nNsVth'], df['r_s'],
+                     df['r_sh'])
     return I * df[loc]
 
 
diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index 31bfb0122b..18f226d9f5 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -822,30 +822,28 @@ def test_singlediode_series_ivcurve(cec_module_params):
     times = pd.DatetimeIndex(start='2015-06-01', periods=3, freq='6H')
     poa_data = pd.Series([0, 400, 800], index=times)
     IL, I0, Rs, Rsh, nNsVth = pvsystem.calcparams_desoto(
-                                         poa_data,
-                                         temp_cell=25,
-                                         alpha_isc=cec_module_params['alpha_sc'],
-                                         module_parameters=cec_module_params,
-                                         EgRef=1.121,
-                                         dEgdT=-0.0002677)
-    
+                                  poa_data, temp_cell=25,
+                                  alpha_isc=cec_module_params['alpha_sc'],
+                                  module_parameters=cec_module_params,
+                                  EgRef=1.121, dEgdT=-0.0002677)
+
     out = pvsystem.singlediode(IL, I0, Rs, Rsh, nNsVth, ivcurve_pnts=3)
-    
-    expected = OrderedDict([('i_sc', array([        0.,  3.01054475,  6.00675648])),
-             ('v_oc', array([         nan,   9.96886962,  10.29530483])),
-             ('i_mp', array([        nan,  2.65191983,  5.28594672])),
-             ('v_mp', array([        nan,  8.33392491,  8.4159707 ])),
-             ('p_mp', array([         nan,  22.10090078,  44.48637274])),
-             ('i_x', array([        nan,  2.88414114,  5.74622046])),
-             ('i_xx', array([        nan,  2.04340914,  3.90007956])),
-             ('v',
-              array([[         nan,          nan,          nan],
-       [  0.        ,   4.98443481,   9.96886962],
-       [  0.        ,   5.14765242,  10.29530483]])),
-             ('i',
-              array([[             nan,              nan,              nan],
-       [  3.01079860e+00,   2.88414114e+00,   3.10862447e-14],
-       [  6.00726296e+00,   5.74622046e+00,   0.00000000e+00]]))])
+
+    expected = OrderedDict([('i_sc', array([0., 3.01054475, 6.00675648])),
+                            ('v_oc', array([0., 9.96886962, 10.29530483])),
+                            ('i_mp', array([0., 2.65191983, 5.28594672])),
+                            ('v_mp', array([0., 8.33392491, 8.4159707])),
+                            ('p_mp', array([0., 22.10090078, 44.48637274])),
+                            ('i_x', array([0., 2.88414114, 5.74622046])),
+                            ('i_xx', array([0., 2.04340914, 3.90007956])),
+                            ('v', array([[0., 0., 0.],
+                                         [0., 4.98443481, 9.96886962],
+                                         [0., 5.14765242, 10.29530483]])),
+                            ('i', array([[0., 0., 0.],
+                                         [3.01079860e+00, 2.88414114e+00,
+                                          3.10862447e-14],
+                                         [6.00726296e+00, 5.74622046e+00,
+                                          0.00000000e+00]]))])
 
     for k, v in out.items():
         assert_allclose(expected[k], v, atol=1e-2)

From cdd2bd980decec2febb819908c7f34b3b0f26427 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Sun, 10 Sep 2017 18:48:41 -0600
Subject: [PATCH 12/22] Conform to existing API

---
 docs/sphinx/source/api.rst             |   7 +-
 docs/sphinx/source/whatsnew/v0.5.1.rst |  26 +-
 pvlib/pvsystem.py                      | 459 +++++++++----------------
 pvlib/test/test_pvsystem.py            | 307 ++++-------------
 4 files changed, 250 insertions(+), 549 deletions(-)

diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst
index 237944458c..4070ed56c3 100644
--- a/docs/sphinx/source/api.rst
+++ b/docs/sphinx/source/api.rst
@@ -202,12 +202,9 @@ Functions relevant for the single diode model.
    :toctree: generated/
 
    pvsystem.calcparams_desoto
-   pvsystem.i_from_v (DEPRECATED: Use pvsystem.sdm_i_from_v instead)
-   pvsystem.sdm_current_sum
-   pvsystem.sdm_i_from_v
-   pvsystem.sdm_v_from_i
+   pvsystem.i_from_v
    pvsystem.singlediode
-   pvsystem.v_from_i (DEPRECATED: Use pvsystem.sdm_v_from_i instead)
+   pvsystem.v_from_i
 
 SAPM model
 ----------
diff --git a/docs/sphinx/source/whatsnew/v0.5.1.rst b/docs/sphinx/source/whatsnew/v0.5.1.rst
index 51e13faf1f..bf928337dc 100644
--- a/docs/sphinx/source/whatsnew/v0.5.1.rst
+++ b/docs/sphinx/source/whatsnew/v0.5.1.rst
@@ -16,33 +16,31 @@ Bug fixes
 Enhancements
 ~~~~~~~~~~~~
 * Improve clearsky.lookup_linke_turbidity speed. (:issue:`368`)
-* Ideal devices supported in single diode model, e.g., 
+* Ideal devices supported in single diode model, e.g.,
   resistance_series = 0 and/or resistance_shunt = numpy.inf (:issue:`340`)
-* Computations for very near ideal devices are more numerically stable in 
-  single diode model (:issue:`340`)
+* `pvsystem.v_from_i` and `pvsystem.i_from_v` computations for near ideal
+  devices are more numerically stable. However, very, very near ideal
+  resistance_series and/or resistance_shunt may still cause issues with the
+  implicit solver (:issue:`340`)
 
 API Changes
 ~~~~~~~~~~~
-* Deprecated `pvsystem.v_from_i` and `pvsystem.i_from_v` functions, which will 
-  be removed in a future release. Instead use `pvsystem.sdm_v_from_i` and
-  `pvsystem.sdm_i_from_v`, which are ~2x slower but can handle ideal devices. 
-  Note that these functions always return a numpy.ndarray (:issue:`340`)
-* Added `pvsystem.sdm_current_sum` function, which computes the sum of currents 
-  at the diode node in the single diode model and is used to improve numerical 
-  stability in `pvsystem.sdm_v_from_i` and `pvsystem.sdm_i_from_v` for very 
-  near ideal devices (:issue:`340`)
+* `pvsystem.v_from_i` and `pvsystem.i_from_v` functions now accept
+  resistance_series = 0 and/or resistance_shunt = numpy.inf as inputs
+  (:issue:`340`)
 
 Documentation
 ~~~~~~~~~~~~~
 * Doc string of modelchain.basic_chain was updated to describe args
   more accurately
+* Doc strings of `singlediode`, `pvsystem.v_from_i`, and `pvsystem.i_from_v`
+  were updated to describe acceptable input arg ranges
 
 Testing
 ~~~~~~~
 * Changed test for clearsky.haurwitz to operate on zenith angles
-* Significant new test cases added for `pvsystem.sdm_v_from_i` and
-  `pvsystem.sdm_i_from_v`, beyond but including cases for deprecated 
-  `pvsystem.v_from_i` and `pvsystem.i_from_v` (:issue:`340`)
+* Significant new test cases added for `pvsystem.v_from_i` and
+  `pvsystem.i_from_v` (:issue:`340`)
 
 Contributors
 ~~~~~~~~~~~~
diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index a534d2a632..0ff6a253df 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -440,18 +440,18 @@ def singlediode(self, photocurrent, saturation_current,
 
     def i_from_v(self, resistance_shunt, resistance_series, nNsVth, voltage,
                  saturation_current, photocurrent):
-        """Wrapper around the :py:func:`sdm_i_from_v` function.
+        """Wrapper around the :py:func:`i_from_v` function.
 
         Parameters
         ----------
-        See pvsystem.sdm_i_from_v for details
+        See pvsystem.i_from_v for details
 
         Returns
         -------
-        See pvsystem.sdm_i_from_v for details
+        See pvsystem.i_from_v for details
         """
-        return sdm_i_from_v(voltage, photocurrent, saturation_current,
-                            nNsVth, resistance_series, resistance_shunt)
+        return i_from_v(resistance_shunt, resistance_series, nNsVth, voltage, \
+                        saturation_current, photocurrent)
 
     # inverter now specified by self.inverter_parameters
     def snlinverter(self, v_dc, p_dc):
@@ -1598,18 +1598,22 @@ def singlediode(photocurrent, saturation_current, resistance_series,
     photocurrent : numeric
         Light-generated current (photocurrent) in amperes under desired
         IV curve conditions. Often abbreviated ``I_L``.
+        0 <= photocurrent
 
     saturation_current : numeric
         Diode saturation current in amperes under desired IV curve
         conditions. Often abbreviated ``I_0``.
+        0 < saturation_current
 
     resistance_series : numeric
         Series resistance in ohms under desired IV curve conditions.
         Often abbreviated ``Rs``.
+        0 <= resistance_series < numpy.inf
 
     resistance_shunt : numeric
         Shunt resistance in ohms under desired IV curve conditions.
         Often abbreviated ``Rsh``.
+        0 <= resistance_shunt <= numpy.inf
 
     nNsVth : numeric
         The product of three components. 1) The usual diode ideal factor
@@ -1619,6 +1623,7 @@ def singlediode(photocurrent, saturation_current, resistance_series,
         ``k*temp_cell/q``, where k is Boltzmann's constant (J/K),
         temp_cell is the temperature of the p-n junction in Kelvin, and
         q is the charge of an electron (coulombs).
+        0 < nNsVth
 
     ivcurve_pnts : None or int, default None
         Number of points in the desired IV curve. If None or 0, no
@@ -1674,13 +1679,13 @@ def singlediode(photocurrent, saturation_current, resistance_series,
     calcparams_desoto
     '''
 
-    # Find short circuit current using Lambert W
-    i_sc = sdm_i_from_v(0., photocurrent, saturation_current, nNsVth,
-                        resistance_series, resistance_shunt)
+    # Compute short circuit current
+    i_sc = i_from_v(resistance_shunt, resistance_series, nNsVth, 0., \
+                    saturation_current, photocurrent)
 
-    # Find open circuit voltage using Lambert W
-    v_oc = sdm_v_from_i(0., photocurrent, saturation_current, nNsVth,
-                        resistance_series, resistance_shunt)
+    # Compute open circuit voltage
+    v_oc = v_from_i(resistance_shunt, resistance_series, nNsVth, 0., \
+                    saturation_current, photocurrent)
 
     params = {'r_sh': resistance_shunt,
               'r_s': resistance_series,
@@ -1692,15 +1697,15 @@ def singlediode(photocurrent, saturation_current, resistance_series,
 
     # Invert the Power-Current curve. Find the current where the inverted power
     # is minimized. This is i_mp. Start the optimization at v_oc/2
-    i_mp = sdm_i_from_v(v_mp, photocurrent, saturation_current, nNsVth,
-                        resistance_series, resistance_shunt)
+    i_mp = i_from_v(resistance_shunt, resistance_series, nNsVth, v_mp, \
+                    saturation_current, photocurrent)
 
     # Find Ix and Ixx using Lambert W
-    i_x = sdm_i_from_v(0.5 * v_oc, photocurrent, saturation_current, nNsVth,
-                       resistance_series, resistance_shunt)
+    i_x = i_from_v(resistance_shunt, resistance_series, nNsVth, 0.5 * v_oc, \
+                   saturation_current, photocurrent)
 
-    i_xx = sdm_i_from_v(0.5 * (v_oc + v_mp), photocurrent, saturation_current,
-                        nNsVth, resistance_series, resistance_shunt)
+    i_xx = i_from_v(resistance_shunt, resistance_series, nNsVth, \
+                    0.5 * (v_oc + v_mp), saturation_current, photocurrent)
 
     out = OrderedDict()
     out['i_sc'] = i_sc
@@ -1713,10 +1718,14 @@ def singlediode(photocurrent, saturation_current, resistance_series,
 
     # create ivcurve
     if ivcurve_pnts:
+        print(v_oc)
         ivcurve_v = (np.asarray(v_oc)[..., np.newaxis] *
                      np.linspace(0, 1, ivcurve_pnts))
-        ivcurve_i = sdm_i_from_v(ivcurve_v.T, photocurrent, saturation_current,
-                                 nNsVth, resistance_series, resistance_shunt).T
+        print(ivcurve_v)
+        ivcurve_i = i_from_v(resistance_shunt, resistance_series, nNsVth, \
+                             ivcurve_v.T, saturation_current, photocurrent).T
+        print(ivcurve_i)
+
         out['v'] = ivcurve_v
         out['i'] = ivcurve_i
 
@@ -1803,27 +1812,35 @@ def _pwr_optfcn(df, loc):
     Function to find power from ``i_from_v``.
     '''
 
-    I = sdm_i_from_v(df[loc], df['i_l'], df['i_0'], df['nNsVth'], df['r_s'],
-                     df['r_sh'])
+    I = i_from_v(df['r_sh'], df['r_s'], df['nNsVth'], df[loc], df['i_0'], \
+                 df['i_l'])
+
     return I * df[loc]
 
 
 def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
              saturation_current, photocurrent):
     '''
-    DEPRECATED: Use sdm_v_from_i() instead.
-
-    Calculates voltage from current per Eq 3 Jain and Kapoor 2004 [1].
+    Computes the device voltage at the given device current using the standard
+     single diode model (SDM) as described in, e.g., Jain and Kapoor 2004 [1].
+    The solution is per Eq 3 of [1] except when resistance_shunt=numpy.inf, in
+     which case the explict solution for voltage is used.
+    Inputs to this function can include scalars and pandas.Series, but it
+     always outputs a float64 numpy.ndarray regardless of input type(s).
+    Ideal device parameters are specified by resistance_shunt=np.inf and
+     resistance_series=0.
 
     Parameters
     ----------
     resistance_shunt : numeric
         Shunt resistance in ohms under desired IV curve conditions.
         Often abbreviated ``Rsh``.
+        0 <= resistance_shunt <= numpy.inf
 
     resistance_series : numeric
         Series resistance in ohms under desired IV curve conditions.
         Often abbreviated ``Rs``.
+        0 <= resistance_series < numpy.inf
 
     nNsVth : numeric
         The product of three components. 1) The usual diode ideal factor
@@ -1833,6 +1850,7 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
         ``k*temp_cell/q``, where k is Boltzmann's constant (J/K),
         temp_cell is the temperature of the p-n junction in Kelvin, and
         q is the charge of an electron (coulombs).
+        0 < nNsVth
 
     current : numeric
         The current in amperes under desired IV curve conditions.
@@ -1840,14 +1858,16 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
     saturation_current : numeric
         Diode saturation current in amperes under desired IV curve
         conditions. Often abbreviated ``I_0``.
+        0 < saturation_current
 
     photocurrent : numeric
         Light-generated current (photocurrent) in amperes under desired
         IV curve conditions. Often abbreviated ``I_L``.
+        0 < photocurrent
 
     Returns
     -------
-    current : np.ndarray or np.float64
+    current : np.ndarray or scalar
 
     References
     ----------
@@ -1861,7 +1881,8 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
         raise ImportError('This function requires scipy')
 
     # asarray turns Series into arrays so that we don't have to worry
-    # about multidimensional broadcasting failing
+    #  about multidimensional broadcasting failing
+    # Note that lambertw function doesn't support float128
     Rsh = np.asarray(resistance_shunt, np.float64)
     Rs = np.asarray(resistance_series, np.float64)
     nNsVth = np.asarray(nNsVth, np.float64)
@@ -1872,48 +1893,95 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
     # This transforms any ideal Rsh=np.inf into Gsh=0., which is generally more numerically stable
     Gsh = 1./Rsh
 
-    # argW cannot be float128
-    argW = I0/(Gsh*nNsVth)*np.exp((-I + IL + I0)/(Gsh*nNsVth))
-    lambertwterm = lambertw(argW).real
-
-    # Calculate using log(argW) in case argW is really big
-    logargW = \
-(np.log(I0) - np.log(Gsh) - np.log(nNsVth) + (-I + IL + I0)/(Gsh*nNsVth))
-
-    # 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 i in range(0, 3):
-        w = w * (1 - np.log(w) + logargW) / (1 + w)
-    lambertwterm_log = w
-
-    lambertwterm = \
-np.where(np.isfinite(lambertwterm), lambertwterm, lambertwterm_log)
-
-    # Eqn. 3 in Jain and Kapoor, 2004
-    V = (IL + I0 - I)/Gsh - I*Rs - nNsVth*lambertwterm
-
-    return V
+    # Intitalize output V (including shape) by solving explicit model with
+    #   Gsh=0, multiplying by np.ones_like(Gsh) identity in order to also
+    #   capture shape of Gsh.
+    V = (nNsVth*np.log1p((IL - I)/I0) - I*Rs)*np.ones_like(Gsh)
+
+    # Record if inputs were all scalar
+    output_is_scalar = np.isscalar(V)
+
+    # Multiply by np.atleast_1d in order to convert scalars to arrays, because
+    #  we need to guarantee the ability to index arrays
+    V = np.atleast_1d(V)
+
+    # Expand Gsh input shape to match output V (if needed)
+    Gsh = Gsh*np.ones_like(V)
+
+    # Determine indices where Gsh>0 requires implicit model solution
+    idx = 0. < Gsh
+
+    # Only compute using LambertW if there are cases with Gsh>0
+    if np.any(idx):
+        # Expand remaining inputs to accomodate common indexing
+        Rs = Rs*np.ones_like(V)
+        nNsVth = nNsVth*np.ones_like(V)
+        I = I*np.ones_like(V)
+        I0 = I0*np.ones_like(V)
+        IL = IL*np.ones_like(V)
+
+        # LambertW argument, argW cannot be float128
+        argW = I0[idx]/(Gsh[idx]*nNsVth[idx])* \
+               np.exp((-I[idx] + IL[idx] + I0[idx])/(Gsh[idx]*nNsVth[idx]))
+
+        # lambertw typically returns complex value with zero imaginary part
+        lambertwterm = lambertw(argW).real
+
+        # Record indices where LambertW input overflowed output
+        idx_w = np.logical_not(np.isfinite(lambertwterm))
+
+        # Only re-compute LambertW if it overflowed
+        if np.any(idx_w):
+            # Calculate using log(argW) in case argW is really big
+            logargW = (np.log(I0[idx]) - np.log(Gsh[idx]) - \
+                      np.log(nNsVth[idx]) + \
+                      (-I[idx] + IL[idx] + I0[idx])/ \
+                      (Gsh[idx]*nNsVth[idx]))[idx_w]
+
+            # 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 i in range(0, 3):
+                w = w * (1 - np.log(w) + logargW) / (1 + w)
+            lambertwterm[idx_w] = w
+
+        # Eqn. 3 in Jain and Kapoor, 2004
+        #  V = -I*(Rs + Rsh) + IL*Rsh - nNsVth*lambertwterm + I0*Rsh
+        # Recasted in terms of Gsh=1/Rsh for better numerical stability.
+        V[idx] = (IL[idx] + I0[idx] - I[idx])/Gsh[idx] - I[idx]*Rs[idx] - \
+                 nNsVth[idx]*lambertwterm
+
+    if output_is_scalar:
+        return np.asscalar(V)
+    else:
+        return V
 
 
 def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
              saturation_current, photocurrent):
     '''
-    DEPRECATED: Use sdm_i_from_v() instead.
-
-    Calculates current from voltage per Eq 2 Jain and Kapoor 2004 [1].
+    Computes the device current at the given device voltage using the standard
+     single diode model (SDM) as described in, e.g., Jain and Kapoor 2004 [1].
+    The solution is per Eq 2 of [1] except when resistance_series=0, in
+     which case the explict solution for current is used.
+    Inputs to this function can include scalars and pandas.Series, but it
+     always outputs a float64 numpy.ndarray regardless of input type(s).
+    Ideal device parameters are specified by resistance_shunt=np.inf and
+     resistance_series=0.
 
     Parameters
     ----------
     resistance_shunt : numeric
         Shunt resistance in ohms under desired IV curve conditions.
         Often abbreviated ``Rsh``.
+        0 <= resistance_shunt <= numpy.inf
 
     resistance_series : numeric
         Series resistance in ohms under desired IV curve conditions.
         Often abbreviated ``Rs``.
+        0 <= resistance_series < numpy.inf
 
     nNsVth : numeric
         The product of three components. 1) The usual diode ideal factor
@@ -1923,6 +1991,7 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
         ``k*temp_cell/q``, where k is Boltzmann's constant (J/K),
         temp_cell is the temperature of the p-n junction in Kelvin, and
         q is the charge of an electron (coulombs).
+        0 < nNsVth
 
     voltage : numeric
         The voltage in Volts under desired IV curve conditions.
@@ -1930,14 +1999,16 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
     saturation_current : numeric
         Diode saturation current in amperes under desired IV curve
         conditions. Often abbreviated ``I_0``.
+        0 < saturation_current
 
     photocurrent : numeric
         Light-generated current (photocurrent) in amperes under desired
         IV curve conditions. Often abbreviated ``I_L``.
+        0 <= photocurrent
 
     Returns
     -------
-    current : np.ndarray or np.float64
+    current : np.ndarray or scalar
 
     References
     ----------
@@ -1951,7 +2022,8 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
         raise ImportError('This function requires scipy')
 
     # asarray turns Series into arrays so that we don't have to worry
-    # about multidimensional broadcasting failing
+    #  about multidimensional broadcasting failing
+    # Note that lambertw function doesn't support float128
     Rsh = np.asarray(resistance_shunt, np.float64)
     Rs = np.asarray(resistance_series, np.float64)
     nNsVth = np.asarray(nNsVth, np.float64)
@@ -1962,239 +2034,50 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
     # This transforms any ideal Rsh=np.inf into Gsh=0., which is generally more numerically stable
     Gsh = 1./Rsh
 
-    # argW cannot be float128
-    argW = \
-Rs*I0/(nNsVth*(Rs*Gsh + 1.))*np.exp((Rs*(IL + I0) + V)/(nNsVth*(Rs*Gsh + 1.)))
-    lambertwterm = lambertw(argW).real
-
-    # Eqn. 4 in Jain and Kapoor, 2004
-    I = (IL + I0 - V*Gsh)/(Rs*Gsh + 1.) - (nNsVth/Rs)*lambertwterm
-
-    return I
-
-
-def sdm_current_sum(V, I, IL, I0, nNsVth, Rs, Rsh):
-    '''
-    Computes the sum of currents at the diode node at electrical steady state 
-    using the standard single diode model (SDM) as described in, e.g., Jain and
-    Kapoor 2004 [1]. An ideal device is specified by Rs=0 and Rsh=numpy.inf.
-    This function behaves as a numpy ufunc, including pandas.Series inputs.
-
-    Parameters
-    ----------
-    V : numeric
-        Device voltage [V]
-
-    I : numeric
-        Device current [A]
-
-    IL : numeric
-        Device photocurrent [A]
-
-    I0 : numeric
-        Device saturation current [A]
-
-    nNsVth : numeric
-        Device thermal voltage [V]
-
-    Rs : numeric
-        Device series resistance [Ohm]
-
-    Rsh : numeric
-        Device shunt resistance [Ohm]
-
-    Returns
-    -------
-    sdm_current_sum : numeric
-        Sum of currents at the diode node in equivalent circuit model [A]
-
-    References
-    ----------
-    [1] A. Jain, A. Kapoor, "Exact analytical solutions of the parameters of 
-    real solar cells using Lambert W-function", Solar Energy Materials and 
-    Solar Cells, 81 (2004) 269-277.
-    '''
-
-    VplusItimesRs = V + I * Rs
-    return IL - I0 * np.expm1(VplusItimesRs / nNsVth) - VplusItimesRs / Rsh - I
-
-
-def sdm_v_from_i(I, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
-    '''
-    Computes the device voltage at the given device current using the standard 
-    single diode model (SDM) as described in, e.g., Jain and Kapoor 2004 [1].
-    An ideal device is specified by Rs=0 and Rsh=numpy.inf and the solution is 
-    per Eq 3 of [1] unless Rsh=numpy.inf gives a more accurate (closed form) 
-    solution as determined by the sum of currents at the diode node, which 
-    should be zero at electrical steady state.
-    Inputs to this function can include scalars and pandas.Series, but it
-    always outputs a float64 numpy.ndarray regardless of input type(s).
-
-    Parameters
-    ----------
-    I : numeric
-        Device current [A]
-
-    IL : numeric
-        Device photocurrent [A]
-
-    I0 : numeric
-        Device saturation current [A]
-
-    nNsVth : numeric
-        Device thermal voltage [V]
-
-    Rs : numeric
-        Device series resistance [Ohm]
-
-    Rsh : numeric
-        Device shunt resistance [Ohm]
-
-    return_meta_dict : boolean scalar (default is False)
-        Return additional computation metadata dictionary
-
-    Returns
-    -------
-    V : float64 numpy.ndarray
-        Device voltage [V]
-
-    meta_dict : dictionary (optional, returned when return_meta_dict=True)
-        Metadata for computation
-
-            meta_dict['sdm_current_sum'] : float64 numpy.ndarray like V
-                Sum of currents at diode node in equivalent circuit model [A]
-
-            meta_dict['inf_Rsh_idx'] : boolean numpy.ndarray like V
-                Indices where infinite shunt resistance gives best solution
-
-    References
-    ----------
-    [1] A. Jain, A. Kapoor, "Exact analytical solutions of the parameters of 
-    real solar cells using Lambert W-function", Solar Energy Materials and 
-    Solar Cells, 81 (2004) 269-277.
-    '''
-
-    # Ensure inputs are all np.ndarray with np.float64 type
-    I = np.asarray(I, np.float64)
-    IL = np.asarray(IL, np.float64)
-    I0 = np.asarray(I0, np.float64)
-    nNsVth = np.asarray(nNsVth, np.float64)
-    Rs = np.asarray(Rs, np.float64)
-    Rsh = np.asarray(Rsh, np.float64)
-
-    # Default computation of V using np.zeros_like keeps Rsh shape info
-    V = nNsVth * (np.log(IL - I - np.zeros_like(I * Rs / Rsh) + I0) -
-                  np.log(I0)) - I * Rs
-    sdm_current_sum_out = sdm_current_sum(V, I, IL, I0, nNsVth, Rs, Rsh)
-
-    # Computation of V using LambertW for provided Rsh
-    # v_from_i is deprecated: move v_from_i code here when/if it is removed.
-    V_lw = v_from_i(Rsh, Rs, nNsVth, I, I0, IL)
-    sdm_current_sum_lw = sdm_current_sum(V_lw, I, IL, I0, nNsVth, Rs, Rsh)
-
-    # Compute selection indices (may be a scalar boolean)
-    finite_Rsh_idx = np.logical_and(np.isfinite(sdm_current_sum_lw),
-                                    np.absolute(sdm_current_sum_lw) <=
-                                    np.absolute(sdm_current_sum_out))
-
-    # These are always np.ndarray
-    V = np.where(finite_Rsh_idx, V_lw, V)
-    sdm_current_sum_out = np.where(finite_Rsh_idx, sdm_current_sum_lw,
-                                   sdm_current_sum_out)
-
-    if return_meta_dict:
-        return V, {'sdm_current_sum': sdm_current_sum_out,
-                   'inf_Rsh_idx': np.array(np.logical_not(finite_Rsh_idx))}
-    else:
-        return V
-
-
-def sdm_i_from_v(V, IL, I0, nNsVth, Rs, Rsh, return_meta_dict=False):
-    '''
-    Computes the device current at the given device voltage using the standard 
-    single diode model (SDM) as described in, e.g., Jain and Kapoor 2004 [1].
-    An ideal device is specified by Rs=0 and Rsh=numpy.inf and the solution is 
-    per Eq 2 of [1] unless Rs=0 gives a more accurate (closed form) 
-    solution as determined by the sum of currents at the diode node, which 
-    should be zero at electrical steady state.
-    Inputs to this function can include scalars and pandas.Series, but it
-    always outputs a float64 numpy.ndarray regardless of input type(s).
-
-    Parameters
-    ----------
-    V : numeric
-        Device voltage [V]
-
-    IL : numeric
-        Device photocurrent [A]
-
-    I0 : numeric
-        Device saturation current [A]
-
-    nNsVth : numeric
-        Device thermal voltage [V]
-
-    Rs : numeric
-        Device series resistance [Ohm]
-
-    Rsh : numeric
-        Device shunt resistance [Ohm]
-
-    return_meta_dict : boolean scalar (default is False)
-        Return additional computation metadata dictionary
-
-    Returns
-    -------
-    I : float64 numpy.ndarray
-        Device current [A]
-
-    meta_dict : dictionary (optional, returned when return_meta_dict=True)
-        Metadata for computation
-
-            meta_dict['sdm_current_sum'] : float64 numpy.ndarray like I
-                Sum of currents at diode node in equivalent circuit model [A]
-
-            meta_dict['zero_Rs_idx'] : boolean numpy.ndarray like I
-                Indices where zero series resistance gives best solution
-
-    References
-    ----------
-    [1] A. Jain, A. Kapoor, "Exact analytical solutions of the parameters of 
-    real solar cells using Lambert W-function", Solar Energy Materials and 
-    Solar Cells, 81 (2004) 269-277.
-    '''
-
-    # Ensure inputs are all np.ndarray with np.float64 type
-    V = np.asarray(V, np.float64)
-    IL = np.asarray(IL, np.float64)
-    I0 = np.asarray(I0, np.float64)
-    nNsVth = np.asarray(nNsVth, np.float64)
-    Rs = np.asarray(Rs, np.float64)
-    Rsh = np.asarray(Rsh, np.float64)
-
-    # Default computation of I using zero_term keeps Rs shape info
-    zero_term = np.zeros_like(Rs)
-    I = IL - I0 * np.expm1((V + zero_term) / nNsVth) - (V + zero_term) / Rsh
-    sdm_current_sum_out = sdm_current_sum(V, I, IL, I0, nNsVth, Rs, Rsh)
-
-    # Computation of I using LambertW for provided Rs
-    # i_from_v is deprecated: move i_from_v code here when/if it is removed.
-    I_lw = i_from_v(Rsh, Rs, nNsVth, V, I0, IL)
-    sdm_current_sum_lw = sdm_current_sum(V, I_lw, IL, I0, nNsVth, Rs, Rsh)
-
-    # Compute selection indices (may be a scalar boolean)
-    nonzero_Rs_idx = np.logical_and(np.isfinite(sdm_current_sum_lw),
-                                    np.absolute(sdm_current_sum_lw) <=
-                                    np.absolute(sdm_current_sum_out))
-
-    # These are always np.ndarray
-    I = np.where(nonzero_Rs_idx, I_lw, I)
-    sdm_current_sum_out = np.where(nonzero_Rs_idx, sdm_current_sum_lw,
-                                   sdm_current_sum_out)
-
-    if return_meta_dict:
-        return I, {'sdm_current_sum': sdm_current_sum_out,
-                   'zero_Rs_idx': np.array(np.logical_not(nonzero_Rs_idx))}
+    # Intitalize output I (including shape) by solving explicit model with
+    #   Rs=0, multiplying by np.ones_like(Rs) identity in order to also
+    #   capture shape of Rs.
+    I = (IL - I0*np.expm1(V/nNsVth) - Gsh*V)*np.ones_like(Rs)
+
+    # Record if inputs were all scalar
+    output_is_scalar = np.isscalar(I)
+
+    # Multiply by np.atleast_1d in order to convert scalars to arrays, because
+    #  we need to guarantee the ability to index arrays
+    I = np.atleast_1d(I)
+
+    # Expand Rs input shape to match output I (if needed)
+    Rs = Rs*np.ones_like(I)
+
+    # Determine indices where Rs>0 requires implicit model solution
+    idx = 0. < Rs
+
+    # Only compute using LambertW if there are cases with Rs>0
+    if np.any(idx):
+        # Expand remaining inputs to accomodate common indexing
+        Gsh = Gsh*np.ones_like(I)
+        nNsVth = nNsVth*np.ones_like(I)
+        V = V*np.ones_like(I)
+        I0 = I0*np.ones_like(I)
+        IL = IL*np.ones_like(I)
+
+        # LambertW argument, argW cannot be float128
+        argW = Rs[idx]*I0[idx]/(nNsVth[idx]*(Rs[idx]*Gsh[idx] + 1.))* \
+               np.exp((Rs[idx]*(IL[idx] + I0[idx]) + V[idx])/ \
+               (nNsVth[idx]*(Rs[idx]*Gsh[idx] + 1.)))
+
+        # lambertw typically returns complex value with zero imaginary part
+        lambertwterm = lambertw(argW).real
+
+        # Eqn. 2 in Jain and Kapoor, 2004
+        #  I = -V/(Rs + Rsh) - (nNsVth/Rs)*lambertwterm + \
+        #      Rsh*(IL + I0)/(Rs + Rsh)
+        # Recasted in terms of Gsh=1/Rsh for better numerical stability.
+        I[idx] = (IL[idx] + I0[idx] - V[idx]*Gsh[idx])/(Rs[idx]*Gsh[idx] + \
+                 1.) - (nNsVth[idx]/Rs[idx])*lambertwterm
+
+    if output_is_scalar:
+        return np.asscalar(I)
     else:
         return I
 
diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index c654411a31..115e058574 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -353,111 +353,6 @@ def test_PVSystem_calcparams_desoto(cec_module_params):
     assert_allclose(nNsVth, 0.473)
 
 
-@pytest.fixture(params=[# Not necessarily I-V curve solutions
-    { # Can handle all python scalar inputs
-     'V' : 40.,
-     'I' : 3.,
-     'IL' : 7.,
-     'I0' : 6.e-7,
-     'nNsVth' : 0.5,
-     'Rs' : 0.1,
-     'Rsh' : 20.,
-     'sdm_current_sum_expected' : np.float64(7. - 6.e-7*np.expm1((40. + 3.*0.1)/0.5) - (40. + 3.*0.1)/20. - 3.)
-    },
-    { # Can handle all rank-0 array inputs
-     'V' : np.array(40.),
-     'I' : np.array(3.),
-     'IL' : np.array(7.),
-     'I0' : np.array(6.e-7),
-     'nNsVth' : np.array(0.5),
-     'Rs' : np.array(0.1),
-     'Rsh' : np.array(20.),
-     'sdm_current_sum_expected' : np.float64(7. - 6.e-7*np.expm1((40. + 3.*0.1)/0.5) - (40. + 3.*0.1)/20. - 3.)
-    },
-    { # Can handle all rank-1 singleton array inputs
-     'V' : np.array([40.]),
-     'I' : np.array([3.]),
-     'IL' : np.array([7.]),
-     'I0' : np.array([6.e-7]),
-     'nNsVth' : np.array([0.5]),
-     'Rs' : np.array([0.1]),
-     'Rsh' : np.array([20.]),
-     'sdm_current_sum_expected' : np.array([7. - 6.e-7*np.expm1((40. + 3.*0.1)/0.5) - (40. + 3.*0.1)/20. - 3.])
-    },
-    { # Can handle all rank-1 non-singleton array inputs
-     'V' : np.array([40., 0. , 0.]),
-     'I' : np.array([0., 0., 3.]),
-     'IL' : np.array([7., 7., 7.]),
-     'I0' : np.array([6.e-7, 6.e-7, 6.e-7]),
-     'nNsVth' : np.array([0.5, 0.5, 0.5]),
-     'Rs' : np.array([0.1, 0.1, 0.1]),
-     'Rsh' : np.array([20., 20., 20.]),
-     'sdm_current_sum_expected' : np.array([7. - 6.e-7*np.expm1(40./0.5) - 40./0.1, 7., 7. - 6.e-7*np.expm1(3.*0.1/0.5) - 3.*0.1/20. - 3.])
-    },
-    { # Can handle mixed inputs with non-singleton Pandas Series
-     'V' : pd.Series([40., 0. , 0.]),
-     'I' : pd.Series([0., 0., 3.]),
-     'IL' : 7.,
-     'I0' : 6.e-7,
-     'nNsVth' : 0.5,
-     'Rs' : 0.1,
-     'Rsh' : 20.,
-     'sdm_current_sum_expected' : pd.Series([7. - 6.e-7*np.expm1(40./0.5) - 40./0.1, 7., 7. - 6.e-7*np.expm1(3.*0.1/0.5) - 3.*0.1/20. - 3.])
-    },
-    { # Can handle mixed inputs with rank-2 arrays
-     'V' : np.array([[40., 0. , 0.], [0., 0. , 40.]]),
-     'I' : np.array([[0., 0., 3.], [3., 0., 0.]]),
-     'IL' : 7.,
-     'I0' : np.full((1,3), 6.e-7),
-     'nNsVth' : np.array(0.5),
-     'Rs' : np.array([0.1]),
-     'Rsh' : np.full((2,3), 20.),
-     'sdm_current_sum_expected' : np.array([[7. - 6.e-7*np.expm1(40./0.5) - 40./0.1, 7., 7. - 6.e-7*np.expm1(3.*0.1/0.5) - 3.*0.1/20. - 3.], \
-                                            [ 7. - 6.e-7*np.expm1(3.*0.1/0.5) - 3.*0.1/20. - 3., 7., 7. - 6.e-7*np.expm1(40./0.5) - 40./0.1]])
-    },
-    { # Can handle infinite shunt resistance with positive series resistance
-     'V' : 40.,
-     'I' : 3.,
-     'IL' : 7.,
-     'I0' : 6.e-7,
-     'nNsVth' : 0.5,
-     'Rs' : 0.1,
-     'Rsh' : np.inf,
-     'sdm_current_sum_expected' : np.float64(7. - 6.e-7*np.expm1((40. + 3.*0.1)/0.5) - 3.)
-    },
-    { # Can handle infinite shunt resistance with zero series resistance
-     'V' : 40.,
-     'I' : 3.,
-     'IL' : 7.,
-     'I0' : 6.e-7,
-     'nNsVth' : 0.5,
-     'Rs' : 0.,
-     'Rsh' : np.inf,
-     'sdm_current_sum_expected' : np.float64(7. - 6.e-7*np.expm1(40./0.5) - 3.)
-    }])
-def fixture_sdm_current_sum(request):
-    return request.param
-
-def test_sdm_current_sum(fixture_sdm_current_sum):
-    # Note: The computation of this function is so straightforward that we do
-    #  NOT extensively verify ufunc behavior
-    
-    # Solution set loaded from fixture
-    V = fixture_sdm_current_sum['V']
-    I = fixture_sdm_current_sum['I']
-    IL = fixture_sdm_current_sum['IL']
-    I0 = fixture_sdm_current_sum['I0']
-    nNsVth = fixture_sdm_current_sum['nNsVth']
-    Rs = fixture_sdm_current_sum['Rs']
-    Rsh = fixture_sdm_current_sum['Rsh']
-    sdm_current_sum_expected = fixture_sdm_current_sum['sdm_current_sum_expected']
-
-    sdm_current_sum = pvsystem.sdm_current_sum(V=V, I=I, IL=IL, I0=I0, nNsVth=nNsVth, Rs=Rs, Rsh=Rsh)
-    assert(isinstance(sdm_current_sum, type(sdm_current_sum_expected)))
-    assert(isinstance(sdm_current_sum.dtype, type(sdm_current_sum_expected.dtype)))
-    assert_array_equal(sdm_current_sum, sdm_current_sum_expected)
-
-
 @pytest.fixture(params=[
     { # Can handle all python scalar inputs
      'Rsh' : 20.,
@@ -466,9 +361,7 @@ def test_sdm_current_sum(fixture_sdm_current_sum):
      'I' : 3.,
      'I0' : 6.e-7,
      'IL' : 7.,
-     'V_expected' : np.array(7.5049875193450521),
-     'sdm_current_sum_expected' : np.array(0.),
-     'inf_Rsh_idx_expected' : np.array(False)
+     'V_expected' : 7.5049875193450521
     },
     { # Can handle all rank-0 array inputs
      'Rsh' : np.array(20.),
@@ -477,9 +370,7 @@ def test_sdm_current_sum(fixture_sdm_current_sum):
      'I' : np.array(3.),
      'I0' : np.array(6.e-7),
      'IL' : np.array(7.),
-     'V_expected' : np.array(7.5049875193450521),
-     'sdm_current_sum_expected' : np.array(0.),
-     'inf_Rsh_idx_expected' : np.array(False)
+     'V_expected' : 7.5049875193450521
     },
     { # Can handle all rank-1 singleton array inputs
      'Rsh' : np.array([20.]),
@@ -488,12 +379,10 @@ def test_sdm_current_sum(fixture_sdm_current_sum):
      'I' : np.array([3.]),
      'I0' : np.array([6.e-7]),
      'IL' : np.array([7.]),
-     'V_expected' : np.array([7.5049875193450521]),
-     'sdm_current_sum_expected' : np.array([0.]),
-     'inf_Rsh_idx_expected' : np.array([False])
+     'V_expected' : np.array([7.5049875193450521])
     },
     { # Can handle all rank-1 non-singleton array inputs with infinite shunt
-      #  resistance, Rsh=inf gives V=Voc=nNsVth*(np.log(IL + I0) - np.log(I0) 
+      #  resistance, Rsh=inf gives V=Voc=nNsVth*(np.log(IL + I0) - np.log(I0)
       #  at I=0
      'Rsh' : np.array([np.inf, 20.]),
      'Rs' : np.array([0.1, 0.1]),
@@ -501,12 +390,10 @@ def test_sdm_current_sum(fixture_sdm_current_sum):
      'I' : np.array([0., 3.]),
      'I0' : np.array([6.e-7, 6.e-7]),
      'IL' : np.array([7., 7.]),
-     'V_expected' : np.array([0.5*(np.log(7. + 6.e-7) - np.log(6.e-7)), 7.5049875193450521]),
-     'sdm_current_sum_expected' : np.array([0., 0.]),
-     'inf_Rsh_idx_expected' : np.array([True, False])
+     'V_expected' : np.array([0.5*(np.log(7. + 6.e-7) - np.log(6.e-7)), 7.5049875193450521])
     },
     { # Can handle mixed inputs with a rank-2 array with infinite shunt
-      #  resistance, Rsh=inf gives V=Voc=nNsVth*(np.log(IL + I0) - np.log(I0) 
+      #  resistance, Rsh=inf gives V=Voc=nNsVth*(np.log(IL + I0) - np.log(I0)
       #  at I=0
      'Rsh' : np.array([[np.inf, np.inf], [np.inf, np.inf]]),
      'Rs' : np.array([0.1]),
@@ -514,11 +401,9 @@ def test_sdm_current_sum(fixture_sdm_current_sum):
      'I' : 0.,
      'I0' : np.array([6.e-7]),
      'IL' : np.array([7.]),
-     'V_expected' : 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))*np.ones((2,2)),
-     'sdm_current_sum_expected' : np.array([[0., 0.], [0., 0.]]),
-     'inf_Rsh_idx_expected' : np.array([[True, True], [True, True]])
+     'V_expected' : 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))*np.ones((2,2))
     },
-    { # Can handle ideal series and shunt, Rsh=inf and Rs=0 give 
+    { # Can handle ideal series and shunt, Rsh=inf and Rs=0 give
       #  V = nNsVth*(np.log(IL - I + I0) - np.log(I0))
      'Rsh' : np.inf,
      'Rs' : 0.,
@@ -526,9 +411,7 @@ def test_sdm_current_sum(fixture_sdm_current_sum):
      'I' : np.array([7., 7./2., 0.]),
      'I0' : 6.e-7,
      'IL' : 7.,
-     'V_expected' : np.array([0., 0.5*(np.log(7. - 7./2. + 6.e-7) - np.log(6.e-7)), 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))]),
-     'sdm_current_sum_expected' : np.array([0., 0., 0.]),
-     'inf_Rsh_idx_expected' : np.array([True, True, True])
+     'V_expected' : np.array([0., 0.5*(np.log(7. - 7./2. + 6.e-7) - np.log(6.e-7)), 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))])
     },
     { # Can handle only ideal series resistance, no closed form solution
      'Rsh' : 20.,
@@ -537,9 +420,7 @@ def test_sdm_current_sum(fixture_sdm_current_sum):
      'I' : 3.,
      'I0' : 6.e-7,
      'IL' : 7.,
-     'V_expected' : np.array(7.804987519345062),
-     'sdm_current_sum_expected' : np.array(0.),
-     'inf_Rsh_idx_expected' : np.array(False)
+     'V_expected' : 7.804987519345062
     },
     { # Can handle all python scalar inputs with big LambertW arg
      'Rsh' : 500.,
@@ -548,9 +429,7 @@ def test_sdm_current_sum(fixture_sdm_current_sum):
      'I' : 0.,
      'I0' : 6.e-10,
      'IL' : 1.2,
-     'V_expected' : np.array(86.320000493521079),
-     'sdm_current_sum_expected' : np.array(0.),
-     'inf_Rsh_idx_expected' : np.array(False)
+     'V_expected' : 86.320000493521079
     },
     { # Can handle all python scalar inputs with bigger LambertW arg
       #  1000 W/m^2 on a Canadian Solar 220M with 20 C ambient temp
@@ -561,40 +440,31 @@ def test_sdm_current_sum(fixture_sdm_current_sum):
      'I' : 0.,
      'I0' : 7.05196029e-08,
      'IL' : 10.491262,
-     'V_expected' : np.array(54.303958833791455),
-     'sdm_current_sum_expected' : np.array(0.),
-     'inf_Rsh_idx_expected' : np.array(False)
+     'V_expected' : 54.303958833791455
     }])
-def fixture_sdm_v_from_i(request):
+def fixture_v_from_i(request):
     return request.param
 
 @requires_scipy
-def test_sdm_v_from_i(fixture_sdm_v_from_i):
+def test_v_from_i(fixture_v_from_i):
     # Solution set loaded from fixture
-    Rsh = fixture_sdm_v_from_i['Rsh']
-    Rs = fixture_sdm_v_from_i['Rs']
-    nNsVth = fixture_sdm_v_from_i['nNsVth']
-    I = fixture_sdm_v_from_i['I']
-    I0 = fixture_sdm_v_from_i['I0']
-    IL = fixture_sdm_v_from_i['IL']
-    V_expected = fixture_sdm_v_from_i['V_expected']
-    sdm_current_sum_expected = fixture_sdm_v_from_i['sdm_current_sum_expected']
-    inf_Rsh_idx_expected = fixture_sdm_v_from_i['inf_Rsh_idx_expected']
-    
+    Rsh = fixture_v_from_i['Rsh']
+    Rs = fixture_v_from_i['Rs']
+    nNsVth = fixture_v_from_i['nNsVth']
+    I = fixture_v_from_i['I']
+    I0 = fixture_v_from_i['I0']
+    IL = fixture_v_from_i['IL']
+    V_expected = fixture_v_from_i['V_expected']
+
     # Convergence criteria
     atol = 1.e-11
-    
-    V = pvsystem.sdm_v_from_i(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, I=I, I0=I0, IL=IL)
+
+    V = pvsystem.v_from_i(Rsh, Rs, nNsVth, I, I0, IL)
     assert(isinstance(V, type(V_expected)))
-    assert(isinstance(V.dtype, type(V_expected.dtype)))
+    if isinstance(V, type(np.ndarray)):
+        assert(isinstance(V.dtype, type(V_expected.dtype)))
+        assert(V.shape == V_expected.shape)
     assert_allclose(V, V_expected, atol=atol)
-    _, meta_dict = pvsystem.sdm_v_from_i(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, I=I, I0=I0, IL=IL, return_meta_dict=True)
-    assert(isinstance(meta_dict['sdm_current_sum'], type(V)))
-    assert(isinstance(meta_dict['sdm_current_sum'].dtype, type(V.dtype)))
-    assert_allclose(meta_dict['sdm_current_sum'], sdm_current_sum_expected, atol=atol)
-    assert(isinstance(meta_dict['inf_Rsh_idx'], type(V)))
-    assert(isinstance(meta_dict['inf_Rsh_idx'].dtype, type(V.dtype)))
-    assert_array_equal(meta_dict['inf_Rsh_idx'], inf_Rsh_idx_expected)
 
 
 @pytest.fixture(params=[
@@ -605,9 +475,7 @@ def test_sdm_v_from_i(fixture_sdm_v_from_i):
      'V' : 40.,
      'I0' : 6.e-7,
      'IL' : 7.,
-     'I_expected' : np.array(-299.746389916),
-     'sdm_current_sum_expected' : np.array(0.),
-     'zero_Rs_idx_expected' : np.array(False)
+     'I_expected' : -299.746389916
     },
     { # Can handle all rank-0 array inputs
      'Rsh' : np.array(20.),
@@ -616,9 +484,7 @@ def test_sdm_v_from_i(fixture_sdm_v_from_i):
      'V' : np.array(40.),
      'I0' : np.array(6.e-7),
      'IL' : np.array(7.),
-     'I_expected' : np.array(-299.746389916),
-     'sdm_current_sum_expected' : np.array(0.),
-     'zero_Rs_idx_expected' : np.array(False)
+     'I_expected' : -299.746389916
     },
     { # Can handle all rank-1 singleton array inputs
      'Rsh' : np.array([20.]),
@@ -627,11 +493,9 @@ def test_sdm_v_from_i(fixture_sdm_v_from_i):
      'V' : np.array([40.]),
      'I0' : np.array([6.e-7]),
      'IL' : np.array([7.]),
-     'I_expected' : np.array([-299.746389916]),
-     'sdm_current_sum_expected' : np.array([0.]),
-     'zero_Rs_idx_expected' : np.array([False])
+     'I_expected' : np.array([-299.746389916])
     },
-    { # Can handle all rank-1 non-singleton array inputs with a zero 
+    { # Can handle all rank-1 non-singleton array inputs with a zero
       #  series resistance, Rs=0 gives I=IL=Isc at V=0
      'Rsh' : np.array([20., 20.]),
      'Rs' : np.array([0., 0.1]),
@@ -639,11 +503,9 @@ def test_sdm_v_from_i(fixture_sdm_v_from_i):
      'V' : np.array([0., 40.]),
      'I0' : np.array([6.e-7, 6.e-7]),
      'IL' : np.array([7., 7.]),
-     'I_expected' : np.array([7., -299.746389916]),
-     'sdm_current_sum_expected' : np.array([0., 0.]),
-     'zero_Rs_idx_expected' : np.array([True, False])
+     'I_expected' : np.array([7., -299.746389916])
     },
-    { # Can handle mixed inputs with a rank-2 array with zero series 
+    { # Can handle mixed inputs with a rank-2 array with zero series
       #  resistance, Rs=0 gives I=IL=Isc at V=0
      'Rsh' : np.array([20.]),
      'Rs' : np.array([[0., 0.], [0., 0.]]),
@@ -651,11 +513,9 @@ def test_sdm_v_from_i(fixture_sdm_v_from_i):
      'V' : 0.,
      'I0' : np.array([6.e-7]),
      'IL' : np.array([7.]),
-     'I_expected' : np.array([[7., 7.], [7., 7.]]),
-     'sdm_current_sum_expected' : np.array([[0., 0.], [0., 0.]]),
-     'zero_Rs_idx_expected' : np.array([[True, True], [True, True]])
+     'I_expected' : np.array([[7., 7.], [7., 7.]])
     },
-    { # Can handle ideal series and shunt, Rsh=inf and Rs=0 give 
+    { # Can handle ideal series and shunt, Rsh=inf and Rs=0 give
       #  V_oc = nNsVth*(np.log(IL + I0) - np.log(I0))
      'Rsh' : np.inf,
      'Rs' : 0.,
@@ -663,9 +523,7 @@ def test_sdm_v_from_i(fixture_sdm_v_from_i):
      'V' : np.array([0., 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))/2., 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))]),
      'I0' : 6.e-7,
      'IL' : 7.,
-     'I_expected' : np.array([7., 7. - 6.e-7*np.expm1((np.log(7. + 6.e-7) - np.log(6.e-7))/2.), 0.]),
-     'sdm_current_sum_expected' : np.array([0., 0., 0.]),
-     'zero_Rs_idx_expected' : np.array([True, True, True])
+     'I_expected' : np.array([7., 7. - 6.e-7*np.expm1((np.log(7. + 6.e-7) - np.log(6.e-7))/2.), 0.])
     },
     { # Can handle only ideal shunt resistance, no closed form solution
      'Rsh' : np.inf,
@@ -674,73 +532,38 @@ def test_sdm_v_from_i(fixture_sdm_v_from_i):
      'V' : 40.,
      'I0' : 6.e-7,
      'IL' : 7.,
-     'I_expected' : np.array(-299.7383436645412),
-     'sdm_current_sum_expected' : np.array(0.),
-     'zero_Rs_idx_expected' : np.array(False)
+     'I_expected' : -299.7383436645412
     }])
-def fixture_sdm_i_from_v(request):
+def fixture_i_from_v(request):
     return request.param
 
 @requires_scipy
-def test_sdm_i_from_v(fixture_sdm_i_from_v):
+def test_i_from_v(fixture_i_from_v):
     # Solution set loaded from fixture
-    Rsh = fixture_sdm_i_from_v['Rsh']
-    Rs = fixture_sdm_i_from_v['Rs']
-    nNsVth = fixture_sdm_i_from_v['nNsVth']
-    V = fixture_sdm_i_from_v['V']
-    I0 = fixture_sdm_i_from_v['I0']
-    IL = fixture_sdm_i_from_v['IL']
-    I_expected = fixture_sdm_i_from_v['I_expected']
-    sdm_current_sum_expected = fixture_sdm_i_from_v['sdm_current_sum_expected']
-    zero_Rs_idx_expected = fixture_sdm_i_from_v['zero_Rs_idx_expected']
+    Rsh = fixture_i_from_v['Rsh']
+    Rs = fixture_i_from_v['Rs']
+    nNsVth = fixture_i_from_v['nNsVth']
+    V = fixture_i_from_v['V']
+    I0 = fixture_i_from_v['I0']
+    IL = fixture_i_from_v['IL']
+    I_expected = fixture_i_from_v['I_expected']
 
     # Convergence criteria
     atol = 1.e-11
-    
-    I = pvsystem.sdm_i_from_v(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL)
+
+    I = pvsystem.i_from_v(Rsh, Rs, nNsVth, V, I0, IL)
     assert(isinstance(I, type(I_expected)))
-    assert(isinstance(I.dtype, type(I_expected.dtype)))
+    if isinstance(I, type(np.ndarray)):
+        assert(isinstance(I.dtype, type(I_expected.dtype)))
+        assert(I.shape == I_expected.shape)
     assert_allclose(I, I_expected, atol=atol)
-    _, meta_dict = pvsystem.sdm_i_from_v(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL, return_meta_dict=True)
-    assert(isinstance(meta_dict['sdm_current_sum'], type(I)))
-    assert(isinstance(meta_dict['sdm_current_sum'].dtype, type(I.dtype)))
-    assert_allclose(meta_dict['sdm_current_sum'], sdm_current_sum_expected, atol=atol)
-    assert(isinstance(meta_dict['zero_Rs_idx'], type(I)))
-    assert(isinstance(meta_dict['zero_Rs_idx'].dtype, type(I.dtype)))
-    assert_array_equal(meta_dict['zero_Rs_idx'], zero_Rs_idx_expected)
-
-
-@requires_scipy
-def test_v_from_i():
-    output = pvsystem.v_from_i(20, .1, .5, 3, 6e-7, 7)
-    assert_allclose(7.5049875193450521, output, atol=1e-5)
-
-
-@requires_scipy
-def test_v_from_i_big():
-    output = pvsystem.v_from_i(500, 10, 4.06, 0, 6e-10, 1.2)
-    assert_allclose(86.320000493521079, output, atol=1e-5)
-
-
-@requires_scipy
-def test_v_from_i_bigger():
-    # 1000 W/m^2 on a Canadian Solar 220M with 20 C ambient temp
-    # github issue 225
-    output = pvsystem.v_from_i(190, 1.065, 2.89, 0, 7.05196029e-08, 10.491262)
-    assert_allclose(54.303958833791455, output, atol=1e-5)
-
-
-@requires_scipy
-def test_i_from_v():
-    output = pvsystem.i_from_v(20, .1, .5, 40, 6e-7, 7)
-    assert_allclose(-299.746389916, output, atol=1e-5)
 
 
 @requires_scipy
 def test_PVSystem_i_from_v():
     system = pvsystem.PVSystem()
     output = system.i_from_v(20, .1, .5, 40, 6e-7, 7)
-    assert_allclose(-299.746389916, output, atol=1e-5)
+    assert_allclose(output, -299.746389916, atol=1e-5)
 
 
 @requires_scipy
@@ -797,7 +620,7 @@ def test_singlediode_floats(sam_data):
         if k in ['i', 'v']:
             assert v is None
         else:
-            assert_allclose(expected[k], v, atol=3)
+            assert_allclose(v, expected[k], atol=3)
 
 
 @requires_scipy
@@ -814,7 +637,7 @@ def test_singlediode_floats_ivcurve():
                 'v': np.array([0.     ,  4.05315,  8.1063])}
     assert isinstance(out, dict)
     for k, v in out.items():
-        assert_allclose(expected[k], v, atol=3)
+        assert_allclose(v, expected[k], atol=3)
 
 
 @requires_scipy
@@ -846,7 +669,7 @@ def test_singlediode_series_ivcurve(cec_module_params):
                                           0.00000000e+00]]))])
 
     for k, v in out.items():
-        assert_allclose(expected[k], v, atol=1e-2)
+        assert_allclose(v, expected[k], atol=1e-2)
 
 
 def test_scale_voltage_current_power(sam_data):
@@ -876,16 +699,16 @@ def test_PVSystem_scale_voltage_current_power():
 
 def test_sapm_celltemp():
     default = pvsystem.sapm_celltemp(900, 5, 20)
-    assert_allclose(43.509, default['temp_cell'], 3)
-    assert_allclose(40.809, default['temp_module'], 3)
+    assert_allclose(default['temp_cell'], 43.509, 3)
+    assert_allclose(default['temp_module'], 40.809, 3)
     assert_frame_equal(default, pvsystem.sapm_celltemp(900, 5, 20,
                                                        [-3.47, -.0594, 3]))
 
 
 def test_sapm_celltemp_dict_like():
     default = pvsystem.sapm_celltemp(900, 5, 20)
-    assert_allclose(43.509, default['temp_cell'], 3)
-    assert_allclose(40.809, default['temp_module'], 3)
+    assert_allclose(default['temp_cell'], 43.509, 3)
+    assert_allclose(default['temp_module'], 40.809, 3)
     model = {'a':-3.47, 'b':-.0594, 'deltaT':3}
     assert_frame_equal(default, pvsystem.sapm_celltemp(900, 5, 20, model))
     model = pd.Series(model)
@@ -1114,7 +937,7 @@ def test_LocalizedPVSystem___repr__():
 def test_pvwatts_dc_scalars():
     expected = 88.65
     out = pvsystem.pvwatts_dc(900, 30, 100, -0.003)
-    assert_allclose(expected, out)
+    assert_allclose(out, expected)
 
 
 @needs_numpy_1_10
@@ -1126,7 +949,7 @@ def test_pvwatts_dc_arrays():
                          [   nan,    nan,    nan],
                          [   nan,  88.65,  88.65]])
     out = pvsystem.pvwatts_dc(irrad_trans, temp_cell, 100, -0.003)
-    assert_allclose(expected, out, equal_nan=True)
+    assert_allclose(out, expected, equal_nan=True)
 
 
 def test_pvwatts_dc_series():
@@ -1140,7 +963,7 @@ def test_pvwatts_dc_series():
 def test_pvwatts_ac_scalars():
     expected = 85.58556604752516
     out = pvsystem.pvwatts_ac(90, 100, 0.95)
-    assert_allclose(expected, out)
+    assert_allclose(out, expected)
 
 
 @needs_numpy_1_10
@@ -1151,7 +974,7 @@ def test_pvwatts_ac_arrays():
                          [ 47.60843624],
                          [ 95.        ]])
     out = pvsystem.pvwatts_ac(pdc, pdc0, 0.95)
-    assert_allclose(expected, out, equal_nan=True)
+    assert_allclose(out, expected, equal_nan=True)
 
 
 def test_pvwatts_ac_series():
@@ -1165,7 +988,7 @@ def test_pvwatts_ac_series():
 def test_pvwatts_losses_default():
     expected = 14.075660688264469
     out = pvsystem.pvwatts_losses()
-    assert_allclose(expected, out)
+    assert_allclose(out, expected)
 
 
 @needs_numpy_1_10
@@ -1173,7 +996,7 @@ def test_pvwatts_losses_arrays():
     expected = np.array([nan, 14.934904])
     age = np.array([nan, 1])
     out = pvsystem.pvwatts_losses(age=age)
-    assert_allclose(expected, out)
+    assert_allclose(out, expected)
 
 
 def test_pvwatts_losses_series():

From 223275a3fa5c1e5d39526c0acaa18542dda3ef49 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Sun, 10 Sep 2017 19:19:31 -0600
Subject: [PATCH 13/22] Run flake8

---
 pvlib/pvsystem.py           |  46 +++---
 pvlib/test/test_pvsystem.py | 283 ++++++++++++++++++------------------
 2 files changed, 169 insertions(+), 160 deletions(-)

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 0ff6a253df..1d5b698f2c 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -450,7 +450,7 @@ def i_from_v(self, resistance_shunt, resistance_series, nNsVth, voltage,
         -------
         See pvsystem.i_from_v for details
         """
-        return i_from_v(resistance_shunt, resistance_series, nNsVth, voltage, \
+        return i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
                         saturation_current, photocurrent)
 
     # inverter now specified by self.inverter_parameters
@@ -1680,11 +1680,11 @@ def singlediode(photocurrent, saturation_current, resistance_series,
     '''
 
     # Compute short circuit current
-    i_sc = i_from_v(resistance_shunt, resistance_series, nNsVth, 0., \
+    i_sc = i_from_v(resistance_shunt, resistance_series, nNsVth, 0.,
                     saturation_current, photocurrent)
 
     # Compute open circuit voltage
-    v_oc = v_from_i(resistance_shunt, resistance_series, nNsVth, 0., \
+    v_oc = v_from_i(resistance_shunt, resistance_series, nNsVth, 0.,
                     saturation_current, photocurrent)
 
     params = {'r_sh': resistance_shunt,
@@ -1697,14 +1697,14 @@ def singlediode(photocurrent, saturation_current, resistance_series,
 
     # Invert the Power-Current curve. Find the current where the inverted power
     # is minimized. This is i_mp. Start the optimization at v_oc/2
-    i_mp = i_from_v(resistance_shunt, resistance_series, nNsVth, v_mp, \
+    i_mp = i_from_v(resistance_shunt, resistance_series, nNsVth, v_mp,
                     saturation_current, photocurrent)
 
     # Find Ix and Ixx using Lambert W
-    i_x = i_from_v(resistance_shunt, resistance_series, nNsVth, 0.5 * v_oc, \
+    i_x = i_from_v(resistance_shunt, resistance_series, nNsVth, 0.5 * v_oc,
                    saturation_current, photocurrent)
 
-    i_xx = i_from_v(resistance_shunt, resistance_series, nNsVth, \
+    i_xx = i_from_v(resistance_shunt, resistance_series, nNsVth,
                     0.5 * (v_oc + v_mp), saturation_current, photocurrent)
 
     out = OrderedDict()
@@ -1722,7 +1722,7 @@ def singlediode(photocurrent, saturation_current, resistance_series,
         ivcurve_v = (np.asarray(v_oc)[..., np.newaxis] *
                      np.linspace(0, 1, ivcurve_pnts))
         print(ivcurve_v)
-        ivcurve_i = i_from_v(resistance_shunt, resistance_series, nNsVth, \
+        ivcurve_i = i_from_v(resistance_shunt, resistance_series, nNsVth,
                              ivcurve_v.T, saturation_current, photocurrent).T
         print(ivcurve_i)
 
@@ -1812,7 +1812,7 @@ def _pwr_optfcn(df, loc):
     Function to find power from ``i_from_v``.
     '''
 
-    I = i_from_v(df['r_sh'], df['r_s'], df['nNsVth'], df[loc], df['i_0'], \
+    I = i_from_v(df['r_sh'], df['r_s'], df['nNsVth'], df[loc], df['i_0'],
                  df['i_l'])
 
     return I * df[loc]
@@ -1890,7 +1890,8 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
     I0 = np.asarray(saturation_current, np.float64)
     IL = np.asarray(photocurrent, np.float64)
 
-    # This transforms any ideal Rsh=np.inf into Gsh=0., which is generally more numerically stable
+    # This transforms any ideal Rsh=np.inf into Gsh=0., which is generally
+    #  more numerically stable
     Gsh = 1./Rsh
 
     # Intitalize output V (including shape) by solving explicit model with
@@ -1921,8 +1922,8 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
         IL = IL*np.ones_like(V)
 
         # LambertW argument, argW cannot be float128
-        argW = I0[idx]/(Gsh[idx]*nNsVth[idx])* \
-               np.exp((-I[idx] + IL[idx] + I0[idx])/(Gsh[idx]*nNsVth[idx]))
+        argW = I0[idx] / (Gsh[idx]*nNsVth[idx]) * \
+            np.exp((-I[idx] + IL[idx] + I0[idx]) / (Gsh[idx]*nNsVth[idx]))
 
         # lambertw typically returns complex value with zero imaginary part
         lambertwterm = lambertw(argW).real
@@ -1933,10 +1934,10 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
         # Only re-compute LambertW if it overflowed
         if np.any(idx_w):
             # Calculate using log(argW) in case argW is really big
-            logargW = (np.log(I0[idx]) - np.log(Gsh[idx]) - \
-                      np.log(nNsVth[idx]) + \
-                      (-I[idx] + IL[idx] + I0[idx])/ \
-                      (Gsh[idx]*nNsVth[idx]))[idx_w]
+            logargW = (np.log(I0[idx]) - np.log(Gsh[idx]) -
+                       np.log(nNsVth[idx]) +
+                       (-I[idx] + IL[idx] + I0[idx]) /
+                       (Gsh[idx] * nNsVth[idx]))[idx_w]
 
             # Three iterations of Newton-Raphson method to solve
             #  w+log(w)=logargW. The initial guess is w=logargW. Where direct
@@ -1951,7 +1952,7 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
         #  V = -I*(Rs + Rsh) + IL*Rsh - nNsVth*lambertwterm + I0*Rsh
         # Recasted in terms of Gsh=1/Rsh for better numerical stability.
         V[idx] = (IL[idx] + I0[idx] - I[idx])/Gsh[idx] - I[idx]*Rs[idx] - \
-                 nNsVth[idx]*lambertwterm
+            nNsVth[idx]*lambertwterm
 
     if output_is_scalar:
         return np.asscalar(V)
@@ -2031,7 +2032,8 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
     I0 = np.asarray(saturation_current, np.float64)
     IL = np.asarray(photocurrent, np.float64)
 
-    # This transforms any ideal Rsh=np.inf into Gsh=0., which is generally more numerically stable
+    # This transforms any ideal Rsh=np.inf into Gsh=0., which is generally
+    #  more numerically stable
     Gsh = 1./Rsh
 
     # Intitalize output I (including shape) by solving explicit model with
@@ -2062,9 +2064,9 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
         IL = IL*np.ones_like(I)
 
         # LambertW argument, argW cannot be float128
-        argW = Rs[idx]*I0[idx]/(nNsVth[idx]*(Rs[idx]*Gsh[idx] + 1.))* \
-               np.exp((Rs[idx]*(IL[idx] + I0[idx]) + V[idx])/ \
-               (nNsVth[idx]*(Rs[idx]*Gsh[idx] + 1.)))
+        argW = Rs[idx]*I0[idx]/(nNsVth[idx]*(Rs[idx]*Gsh[idx] + 1.)) * \
+            np.exp((Rs[idx]*(IL[idx] + I0[idx]) + V[idx]) /
+                   (nNsVth[idx]*(Rs[idx]*Gsh[idx] + 1.)))
 
         # lambertw typically returns complex value with zero imaginary part
         lambertwterm = lambertw(argW).real
@@ -2073,8 +2075,8 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
         #  I = -V/(Rs + Rsh) - (nNsVth/Rs)*lambertwterm + \
         #      Rsh*(IL + I0)/(Rs + Rsh)
         # Recasted in terms of Gsh=1/Rsh for better numerical stability.
-        I[idx] = (IL[idx] + I0[idx] - V[idx]*Gsh[idx])/(Rs[idx]*Gsh[idx] + \
-                 1.) - (nNsVth[idx]/Rs[idx])*lambertwterm
+        I[idx] = (IL[idx] + I0[idx] - V[idx]*Gsh[idx]) / \
+            (Rs[idx]*Gsh[idx] + 1.) - (nNsVth[idx]/Rs[idx])*lambertwterm
 
     if output_is_scalar:
         return np.asscalar(I)
diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index 115e058574..70d6a4d25d 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -9,7 +9,7 @@
 
 import pytest
 from pandas.util.testing import assert_series_equal, assert_frame_equal
-from numpy.testing import assert_allclose, assert_array_equal
+from numpy.testing import assert_allclose
 
 from pvlib import tmy
 from pvlib import pvsystem
@@ -354,97 +354,101 @@ def test_PVSystem_calcparams_desoto(cec_module_params):
 
 
 @pytest.fixture(params=[
-    { # Can handle all python scalar inputs
-     'Rsh' : 20.,
-     'Rs' : 0.1,
-     'nNsVth' : 0.5,
-     'I' : 3.,
-     'I0' : 6.e-7,
-     'IL' : 7.,
-     'V_expected' : 7.5049875193450521
+    {  # Can handle all python scalar inputs
+     'Rsh': 20.,
+     'Rs': 0.1,
+     'nNsVth': 0.5,
+     'I': 3.,
+     'I0': 6.e-7,
+     'IL': 7.,
+     'V_expected': 7.5049875193450521
     },
-    { # Can handle all rank-0 array inputs
-     'Rsh' : np.array(20.),
-     'Rs' : np.array(0.1),
-     'nNsVth' : np.array(0.5),
-     'I' : np.array(3.),
-     'I0' : np.array(6.e-7),
-     'IL' : np.array(7.),
-     'V_expected' : 7.5049875193450521
+    {  # Can handle all rank-0 array inputs
+     'Rsh': np.array(20.),
+     'Rs': np.array(0.1),
+     'nNsVth': np.array(0.5),
+     'I': np.array(3.),
+     'I0': np.array(6.e-7),
+     'IL': np.array(7.),
+     'V_expected': 7.5049875193450521
     },
-    { # Can handle all rank-1 singleton array inputs
-     'Rsh' : np.array([20.]),
-     'Rs' : np.array([0.1]),
-     'nNsVth' : np.array([0.5]),
-     'I' : np.array([3.]),
-     'I0' : np.array([6.e-7]),
-     'IL' : np.array([7.]),
-     'V_expected' : np.array([7.5049875193450521])
+    {  # Can handle all rank-1 singleton array inputs
+     'Rsh': np.array([20.]),
+     'Rs': np.array([0.1]),
+     'nNsVth': np.array([0.5]),
+     'I': np.array([3.]),
+     'I0': np.array([6.e-7]),
+     'IL': np.array([7.]),
+     'V_expected': np.array([7.5049875193450521])
     },
-    { # Can handle all rank-1 non-singleton array inputs with infinite shunt
+    {  # Can handle all rank-1 non-singleton array inputs with infinite shunt
       #  resistance, Rsh=inf gives V=Voc=nNsVth*(np.log(IL + I0) - np.log(I0)
       #  at I=0
-     'Rsh' : np.array([np.inf, 20.]),
-     'Rs' : np.array([0.1, 0.1]),
-     'nNsVth' : np.array([0.5, 0.5]),
-     'I' : np.array([0., 3.]),
-     'I0' : np.array([6.e-7, 6.e-7]),
-     'IL' : np.array([7., 7.]),
-     'V_expected' : np.array([0.5*(np.log(7. + 6.e-7) - np.log(6.e-7)), 7.5049875193450521])
+      'Rsh': np.array([np.inf, 20.]),
+      'Rs': np.array([0.1, 0.1]),
+      'nNsVth': np.array([0.5, 0.5]),
+      'I': np.array([0., 3.]),
+      'I0': np.array([6.e-7, 6.e-7]),
+      'IL': np.array([7., 7.]),
+      'V_expected': np.array([0.5*(np.log(7. + 6.e-7) - np.log(6.e-7)),
+                             7.5049875193450521])
     },
-    { # Can handle mixed inputs with a rank-2 array with infinite shunt
+    {  # Can handle mixed inputs with a rank-2 array with infinite shunt
       #  resistance, Rsh=inf gives V=Voc=nNsVth*(np.log(IL + I0) - np.log(I0)
       #  at I=0
-     'Rsh' : np.array([[np.inf, np.inf], [np.inf, np.inf]]),
-     'Rs' : np.array([0.1]),
-     'nNsVth' : np.array(0.5),
-     'I' : 0.,
-     'I0' : np.array([6.e-7]),
-     'IL' : np.array([7.]),
-     'V_expected' : 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))*np.ones((2,2))
+      'Rsh': np.array([[np.inf, np.inf], [np.inf, np.inf]]),
+      'Rs': np.array([0.1]),
+      'nNsVth': np.array(0.5),
+      'I': 0.,
+      'I0': np.array([6.e-7]),
+      'IL': np.array([7.]),
+      'V_expected': 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))*np.ones((2, 2))
     },
-    { # Can handle ideal series and shunt, Rsh=inf and Rs=0 give
+    {  # Can handle ideal series and shunt, Rsh=inf and Rs=0 give
       #  V = nNsVth*(np.log(IL - I + I0) - np.log(I0))
-     'Rsh' : np.inf,
-     'Rs' : 0.,
-     'nNsVth' : 0.5,
-     'I' : np.array([7., 7./2., 0.]),
-     'I0' : 6.e-7,
-     'IL' : 7.,
-     'V_expected' : np.array([0., 0.5*(np.log(7. - 7./2. + 6.e-7) - np.log(6.e-7)), 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))])
+      'Rsh': np.inf,
+      'Rs': 0.,
+      'nNsVth': 0.5,
+      'I': np.array([7., 7./2., 0.]),
+      'I0': 6.e-7,
+      'IL': 7.,
+      'V_expected': np.array([0., 0.5*(np.log(7. - 7./2. + 6.e-7) -
+                              np.log(6.e-7)), 0.5*(np.log(7. + 6.e-7) -
+                              np.log(6.e-7))])
     },
-    { # Can handle only ideal series resistance, no closed form solution
-     'Rsh' : 20.,
-     'Rs' : 0.,
-     'nNsVth' : 0.5,
-     'I' : 3.,
-     'I0' : 6.e-7,
-     'IL' : 7.,
-     'V_expected' : 7.804987519345062
+    {  # Can handle only ideal series resistance, no closed form solution
+      'Rsh': 20.,
+      'Rs': 0.,
+      'nNsVth': 0.5,
+      'I': 3.,
+      'I0': 6.e-7,
+      'IL': 7.,
+      'V_expected': 7.804987519345062
     },
-    { # Can handle all python scalar inputs with big LambertW arg
-     'Rsh' : 500.,
-     'Rs' : 10.,
-     'nNsVth' : 4.06,
-     'I' : 0.,
-     'I0' : 6.e-10,
-     'IL' : 1.2,
-     'V_expected' : 86.320000493521079
+    {  # Can handle all python scalar inputs with big LambertW arg
+      'Rsh': 500.,
+      'Rs': 10.,
+      'nNsVth': 4.06,
+      'I': 0.,
+      'I0': 6.e-10,
+      'IL': 1.2,
+      'V_expected': 86.320000493521079
     },
-    { # Can handle all python scalar inputs with bigger LambertW arg
+    {  # Can handle all python scalar inputs with bigger LambertW arg
       #  1000 W/m^2 on a Canadian Solar 220M with 20 C ambient temp
       #  github issue 225
-     'Rsh' : 190.,
-     'Rs' : 1.065,
-     'nNsVth' : 2.89,
-     'I' : 0.,
-     'I0' : 7.05196029e-08,
-     'IL' : 10.491262,
-     'V_expected' : 54.303958833791455
+      'Rsh': 190.,
+      'Rs': 1.065,
+      'nNsVth': 2.89,
+      'I': 0.,
+      'I0': 7.05196029e-08,
+      'IL': 10.491262,
+      'V_expected': 54.303958833791455
     }])
 def fixture_v_from_i(request):
     return request.param
 
+
 @requires_scipy
 def test_v_from_i(fixture_v_from_i):
     # Solution set loaded from fixture
@@ -468,75 +472,78 @@ def test_v_from_i(fixture_v_from_i):
 
 
 @pytest.fixture(params=[
-    { # Can handle all python scalar inputs
-     'Rsh' : 20.,
-     'Rs' : 0.1,
-     'nNsVth' : 0.5,
-     'V' : 40.,
-     'I0' : 6.e-7,
-     'IL' : 7.,
-     'I_expected' : -299.746389916
+    {  # Can handle all python scalar inputs
+      'Rsh': 20.,
+      'Rs': 0.1,
+      'nNsVth': 0.5,
+      'V': 40.,
+      'I0': 6.e-7,
+      'IL': 7.,
+      'I_expected': -299.746389916
     },
-    { # Can handle all rank-0 array inputs
-     'Rsh' : np.array(20.),
-     'Rs' : np.array(0.1),
-     'nNsVth' : np.array(0.5),
-     'V' : np.array(40.),
-     'I0' : np.array(6.e-7),
-     'IL' : np.array(7.),
-     'I_expected' : -299.746389916
+    {  # Can handle all rank-0 array inputs
+      'Rsh': np.array(20.),
+      'Rs': np.array(0.1),
+      'nNsVth': np.array(0.5),
+      'V': np.array(40.),
+      'I0': np.array(6.e-7),
+      'IL': np.array(7.),
+      'I_expected': -299.746389916
     },
-    { # Can handle all rank-1 singleton array inputs
-     'Rsh' : np.array([20.]),
-     'Rs' : np.array([0.1]),
-     'nNsVth' : np.array([0.5]),
-     'V' : np.array([40.]),
-     'I0' : np.array([6.e-7]),
-     'IL' : np.array([7.]),
-     'I_expected' : np.array([-299.746389916])
+    {  # Can handle all rank-1 singleton array inputs
+      'Rsh': np.array([20.]),
+      'Rs': np.array([0.1]),
+      'nNsVth': np.array([0.5]),
+      'V': np.array([40.]),
+      'I0': np.array([6.e-7]),
+      'IL': np.array([7.]),
+      'I_expected': np.array([-299.746389916])
     },
-    { # Can handle all rank-1 non-singleton array inputs with a zero
+    {  # Can handle all rank-1 non-singleton array inputs with a zero
       #  series resistance, Rs=0 gives I=IL=Isc at V=0
-     'Rsh' : np.array([20., 20.]),
-     'Rs' : np.array([0., 0.1]),
-     'nNsVth' : np.array([0.5, 0.5]),
-     'V' : np.array([0., 40.]),
-     'I0' : np.array([6.e-7, 6.e-7]),
-     'IL' : np.array([7., 7.]),
-     'I_expected' : np.array([7., -299.746389916])
+      'Rsh': np.array([20., 20.]),
+      'Rs': np.array([0., 0.1]),
+      'nNsVth': np.array([0.5, 0.5]),
+      'V': np.array([0., 40.]),
+      'I0': np.array([6.e-7, 6.e-7]),
+      'IL': np.array([7., 7.]),
+      'I_expected': np.array([7., -299.746389916])
     },
-    { # Can handle mixed inputs with a rank-2 array with zero series
+    {  # Can handle mixed inputs with a rank-2 array with zero series
       #  resistance, Rs=0 gives I=IL=Isc at V=0
-     'Rsh' : np.array([20.]),
-     'Rs' : np.array([[0., 0.], [0., 0.]]),
-     'nNsVth' : np.array(0.5),
-     'V' : 0.,
-     'I0' : np.array([6.e-7]),
-     'IL' : np.array([7.]),
-     'I_expected' : np.array([[7., 7.], [7., 7.]])
+      'Rsh': np.array([20.]),
+      'Rs': np.array([[0., 0.], [0., 0.]]),
+      'nNsVth': np.array(0.5),
+      'V': 0.,
+      'I0': np.array([6.e-7]),
+      'IL': np.array([7.]),
+      'I_expected': np.array([[7., 7.], [7., 7.]])
     },
-    { # Can handle ideal series and shunt, Rsh=inf and Rs=0 give
+    {  # Can handle ideal series and shunt, Rsh=inf and Rs=0 give
       #  V_oc = nNsVth*(np.log(IL + I0) - np.log(I0))
-     'Rsh' : np.inf,
-     'Rs' : 0.,
-     'nNsVth' : 0.5,
-     'V' : np.array([0., 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))/2., 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))]),
-     'I0' : 6.e-7,
-     'IL' : 7.,
-     'I_expected' : np.array([7., 7. - 6.e-7*np.expm1((np.log(7. + 6.e-7) - np.log(6.e-7))/2.), 0.])
+      'Rsh': np.inf,
+      'Rs': 0.,
+      'nNsVth': 0.5,
+      'V': np.array([0., 0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))/2.,
+                     0.5*(np.log(7. + 6.e-7) - np.log(6.e-7))]),
+      'I0': 6.e-7,
+      'IL': 7.,
+      'I_expected': np.array([7., 7. - 6.e-7*np.expm1((np.log(7. + 6.e-7) -
+                              np.log(6.e-7))/2.), 0.])
     },
-    { # Can handle only ideal shunt resistance, no closed form solution
-     'Rsh' : np.inf,
-     'Rs' : 0.1,
-     'nNsVth' : 0.5,
-     'V' : 40.,
-     'I0' : 6.e-7,
-     'IL' : 7.,
-     'I_expected' : -299.7383436645412
+    {  # Can handle only ideal shunt resistance, no closed form solution
+      'Rsh': np.inf,
+      'Rs': 0.1,
+      'nNsVth': 0.5,
+      'V': 40.,
+      'I0': 6.e-7,
+      'IL': 7.,
+      'I_expected': -299.7383436645412
     }])
 def fixture_i_from_v(request):
     return request.param
 
+
 @requires_scipy
 def test_i_from_v(fixture_i_from_v):
     # Solution set loaded from fixture
@@ -633,8 +640,8 @@ def test_singlediode_floats_ivcurve():
                 'i_x': 6.7556075876880621,
                 'i_sc': 6.9646747613963198,
                 'v_mp': 6.221535886625464,
-                'i': np.array([6.965172e+00,   6.755882e+00,   2.575717e-14]),
-                'v': np.array([0.     ,  4.05315,  8.1063])}
+                'i': np.array([6.965172e+00, 6.755882e+00, 2.575717e-14]),
+                'v': np.array([0., 4.05315, 8.1063])}
     assert isinstance(out, dict)
     for k, v in out.items():
         assert_allclose(v, expected[k], atol=3)
@@ -709,7 +716,7 @@ def test_sapm_celltemp_dict_like():
     default = pvsystem.sapm_celltemp(900, 5, 20)
     assert_allclose(default['temp_cell'], 43.509, 3)
     assert_allclose(default['temp_module'], 40.809, 3)
-    model = {'a':-3.47, 'b':-.0594, 'deltaT':3}
+    model = {'a': -3.47, 'b': -.0594, 'deltaT': 3}
     assert_frame_equal(default, pvsystem.sapm_celltemp(900, 5, 20, model))
     model = pd.Series(model)
     assert_frame_equal(default, pvsystem.sapm_celltemp(900, 5, 20, model))
@@ -945,9 +952,9 @@ def test_pvwatts_dc_arrays():
     irrad_trans = np.array([np.nan, 900, 900])
     temp_cell = np.array([30, np.nan, 30])
     irrad_trans, temp_cell = np.meshgrid(irrad_trans, temp_cell)
-    expected = np.array([[   nan,  88.65,  88.65],
-                         [   nan,    nan,    nan],
-                         [   nan,  88.65,  88.65]])
+    expected = np.array([[nan,  88.65,  88.65],
+                         [nan,    nan,    nan],
+                         [nan,  88.65,  88.65]])
     out = pvsystem.pvwatts_dc(irrad_trans, temp_cell, 100, -0.003)
     assert_allclose(out, expected, equal_nan=True)
 
@@ -970,9 +977,9 @@ def test_pvwatts_ac_scalars():
 def test_pvwatts_ac_arrays():
     pdc = np.array([[np.nan], [50], [100]])
     pdc0 = 100
-    expected = np.array([[         nan],
-                         [ 47.60843624],
-                         [ 95.        ]])
+    expected = np.array([[nan],
+                         [47.60843624],
+                         [95.]])
     out = pvsystem.pvwatts_ac(pdc, pdc0, 0.95)
     assert_allclose(out, expected, equal_nan=True)
 

From 6b50ad9e3b31f40f08a4be39ca21417c933725fb Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Sun, 10 Sep 2017 19:38:17 -0600
Subject: [PATCH 14/22] Implement some code quality suggestions

---
 pvlib/pvsystem.py             |  4 +++-
 pvlib/test/test_irradiance.py | 13 +++++++------
 2 files changed, 10 insertions(+), 7 deletions(-)

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 1d5b698f2c..51af03fc08 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1823,6 +1823,7 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
     '''
     Computes the device voltage at the given device current using the standard
      single diode model (SDM) as described in, e.g., Jain and Kapoor 2004 [1].
+
     The solution is per Eq 3 of [1] except when resistance_shunt=numpy.inf, in
      which case the explict solution for voltage is used.
     Inputs to this function can include scalars and pandas.Series, but it
@@ -1944,7 +1945,7 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
             #  evaluation (above) results in NaN from overflow, 3 iterations
             #  of Newton's method gives approximately 8 digits of precision.
             w = logargW
-            for i in range(0, 3):
+            for _ in range(0, 3):
                 w = w * (1 - np.log(w) + logargW) / (1 + w)
             lambertwterm[idx_w] = w
 
@@ -1965,6 +1966,7 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
     '''
     Computes the device current at the given device voltage using the standard
      single diode model (SDM) as described in, e.g., Jain and Kapoor 2004 [1].
+
     The solution is per Eq 2 of [1] except when resistance_series=0, in
      which case the explict solution for current is used.
     Inputs to this function can include scalars and pandas.Series, but it
diff --git a/pvlib/test/test_irradiance.py b/pvlib/test/test_irradiance.py
index 8c87f29eea..4a6fe224db 100755
--- a/pvlib/test/test_irradiance.py
+++ b/pvlib/test/test_irradiance.py
@@ -425,12 +425,13 @@ def test_dni():
 
 
 @pytest.mark.parametrize(
-    'surface_tilt,surface_azimuth,solar_zenith,solar_azimuth,aoi_expected,aoi_proj_expected', [
-    (0, 0, 0, 0, 0, 1),
-    (30, 180, 30, 180, 0, 1),
-    (30, 180, 150, 0, 180, -1),
-    (90, 0, 30, 60, 75.5224878, 0.25),
-    (90, 0, 30, 170, 119.4987042, -0.4924038)])
+    'surface_tilt,surface_azimuth,solar_zenith,' +
+    'solar_azimuth,aoi_expected,aoi_proj_expected',
+    [(0, 0, 0, 0, 0, 1),
+     (30, 180, 30, 180, 0, 1),
+     (30, 180, 150, 0, 180, -1),
+     (90, 0, 30, 60, 75.5224878, 0.25),
+     (90, 0, 30, 170, 119.4987042, -0.4924038)])
 def test_aoi_and_aoi_projection(surface_tilt, surface_azimuth, solar_zenith,
                                 solar_azimuth, aoi_expected,
                                 aoi_proj_expected):

From e4e09be86b706533abfcde4f25778947c6417d5e Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Sun, 10 Sep 2017 19:44:57 -0600
Subject: [PATCH 15/22] Remove extraneous print statements

---
 pvlib/pvsystem.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 51af03fc08..831dd58628 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1718,13 +1718,11 @@ def singlediode(photocurrent, saturation_current, resistance_series,
 
     # create ivcurve
     if ivcurve_pnts:
-        print(v_oc)
         ivcurve_v = (np.asarray(v_oc)[..., np.newaxis] *
                      np.linspace(0, 1, ivcurve_pnts))
-        print(ivcurve_v)
+
         ivcurve_i = i_from_v(resistance_shunt, resistance_series, nNsVth,
                              ivcurve_v.T, saturation_current, photocurrent).T
-        print(ivcurve_i)
 
         out['v'] = ivcurve_v
         out['i'] = ivcurve_i

From 37484f128888ea65f88757fa63d76b6928d3e3d2 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Sun, 10 Sep 2017 20:13:06 -0600
Subject: [PATCH 16/22] Better docstrings

---
 pvlib/pvsystem.py | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 831dd58628..286c6731f0 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1819,11 +1819,12 @@ def _pwr_optfcn(df, loc):
 def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
              saturation_current, photocurrent):
     '''
-    Computes the device voltage at the given device current using the standard
-     single diode model (SDM) as described in, e.g., Jain and Kapoor 2004 [1].
+    Device voltage at the given device current for the single diode model.
 
-    The solution is per Eq 3 of [1] except when resistance_shunt=numpy.inf, in
-     which case the explict solution for voltage is used.
+    Uses the single diode model (SDM) as described in, e.g.,
+     Jain and Kapoor 2004 [1].
+    The solution is per Eq 3 of [1] except when resistance_shunt=numpy.inf,
+     in which case the explict solution for voltage is used.
     Inputs to this function can include scalars and pandas.Series, but it
      always outputs a float64 numpy.ndarray regardless of input type(s).
     Ideal device parameters are specified by resistance_shunt=np.inf and
@@ -1962,11 +1963,12 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
 def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
              saturation_current, photocurrent):
     '''
-    Computes the device current at the given device voltage using the standard
-     single diode model (SDM) as described in, e.g., Jain and Kapoor 2004 [1].
+    Device current at the given device voltage for the single diode model.
 
-    The solution is per Eq 2 of [1] except when resistance_series=0, in
-     which case the explict solution for current is used.
+    Uses the single diode model (SDM) as described in, e.g.,
+     Jain and Kapoor 2004 [1].
+    The solution is per Eq 2 of [1] except when resistance_series=0,
+     in which case the explict solution for current is used.
     Inputs to this function can include scalars and pandas.Series, but it
      always outputs a float64 numpy.ndarray regardless of input type(s).
     Ideal device parameters are specified by resistance_shunt=np.inf and

From f9c4c829103039879f8d22e49184b5ad47c323fa Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Sun, 10 Sep 2017 20:35:42 -0600
Subject: [PATCH 17/22] Fix parameter ranges in docstrings

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

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 286c6731f0..8082370ad6 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1613,7 +1613,7 @@ def singlediode(photocurrent, saturation_current, resistance_series,
     resistance_shunt : numeric
         Shunt resistance in ohms under desired IV curve conditions.
         Often abbreviated ``Rsh``.
-        0 <= resistance_shunt <= numpy.inf
+        0 < resistance_shunt <= numpy.inf
 
     nNsVth : numeric
         The product of three components. 1) The usual diode ideal factor
@@ -1835,7 +1835,7 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
     resistance_shunt : numeric
         Shunt resistance in ohms under desired IV curve conditions.
         Often abbreviated ``Rsh``.
-        0 <= resistance_shunt <= numpy.inf
+        0 < resistance_shunt <= numpy.inf
 
     resistance_series : numeric
         Series resistance in ohms under desired IV curve conditions.
@@ -1863,7 +1863,7 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
     photocurrent : numeric
         Light-generated current (photocurrent) in amperes under desired
         IV curve conditions. Often abbreviated ``I_L``.
-        0 < photocurrent
+        0 <= photocurrent
 
     Returns
     -------
@@ -1979,7 +1979,7 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
     resistance_shunt : numeric
         Shunt resistance in ohms under desired IV curve conditions.
         Often abbreviated ``Rsh``.
-        0 <= resistance_shunt <= numpy.inf
+        0 < resistance_shunt <= numpy.inf
 
     resistance_series : numeric
         Series resistance in ohms under desired IV curve conditions.

From 26eabdeb4132d2a1146bb28eeb58efb951b93e3f Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Mon, 11 Sep 2017 12:37:47 -0600
Subject: [PATCH 18/22] Add test that overflows lambertw arg

---
 pvlib/test/test_pvsystem.py | 13 ++++++++++++-
 1 file changed, 12 insertions(+), 1 deletion(-)

diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index 70d6a4d25d..53719af471 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -436,7 +436,7 @@ def test_PVSystem_calcparams_desoto(cec_module_params):
     },
     {  # Can handle all python scalar inputs with bigger LambertW arg
       #  1000 W/m^2 on a Canadian Solar 220M with 20 C ambient temp
-      #  github issue 225
+      #  github issue 225 (this appears to be from PR 226 not issue 225)
       'Rsh': 190.,
       'Rs': 1.065,
       'nNsVth': 2.89,
@@ -444,6 +444,17 @@ def test_PVSystem_calcparams_desoto(cec_module_params):
       'I0': 7.05196029e-08,
       'IL': 10.491262,
       'V_expected': 54.303958833791455
+    },
+    {  # Can handle all python scalar inputs with bigger LambertW arg
+      #  1000 W/m^2 on a Canadian Solar 220M with 20 C ambient temp
+      #  github issue 225
+      'Rsh': 381.68,
+      'Rs': 1.065,
+      'nNsVth': 2.681527737715915,
+      'I': 0.,
+      'I0': 1.8739027472625636e-09,
+      'IL': 5.1366949999999996,
+      'V_expected': 58.19323124611128
     }])
 def fixture_v_from_i(request):
     return request.param

From 4d7754d736f5a1b6c9cf75646a71775dcb758d05 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Mon, 11 Sep 2017 13:09:39 -0600
Subject: [PATCH 19/22] Add test for mixed solution types logic

---
 pvlib/test/test_pvsystem.py | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index 53719af471..2259b30623 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -455,6 +455,16 @@ def test_PVSystem_calcparams_desoto(cec_module_params):
       'I0': 1.8739027472625636e-09,
       'IL': 5.1366949999999996,
       'V_expected': 58.19323124611128
+    },
+    {  # Verify mixed solution type indexing logic
+      'Rsh': np.array([np.inf, 190., 381.68]),
+      'Rs': 1.065,
+      'nNsVth': np.array([2.89, 2.89, 2.681527737715915]),
+      'I': 0.,
+      'I0': np.array([7.05196029e-08, 7.05196029e-08, 1.8739027472625636e-09]),
+      'IL': np.array([10.491262, 10.491262, 5.1366949999999996]),
+      'V_expected': np.array([2.89*np.log1p(10.491262/7.05196029e-08),
+                              54.303958833791455, 58.19323124611128])
     }])
 def fixture_v_from_i(request):
     return request.param

From 03ae51e7e15be6a448f22db8584f4d23438323be Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Mon, 11 Sep 2017 21:50:56 -0600
Subject: [PATCH 20/22] Use broadcast_arrays for cleaner code

---
 pvlib/pvsystem.py           | 163 ++++++++++++++++--------------------
 pvlib/test/test_pvsystem.py |   4 +-
 2 files changed, 76 insertions(+), 91 deletions(-)

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 8082370ad6..cddc33a3e2 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1825,10 +1825,11 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
      Jain and Kapoor 2004 [1].
     The solution is per Eq 3 of [1] except when resistance_shunt=numpy.inf,
      in which case the explict solution for voltage is used.
-    Inputs to this function can include scalars and pandas.Series, but it
-     always outputs a float64 numpy.ndarray regardless of input type(s).
     Ideal device parameters are specified by resistance_shunt=np.inf and
      resistance_series=0.
+    Inputs to this function can include scalars and pandas.Series, but it is
+     the caller's responsibility to ensure that the arguments are all float64
+     and within the proper ranges.
 
     Parameters
     ----------
@@ -1880,64 +1881,56 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
     except ImportError:
         raise ImportError('This function requires scipy')
 
-    # asarray turns Series into arrays so that we don't have to worry
-    #  about multidimensional broadcasting failing
-    # Note that lambertw function doesn't support float128
-    Rsh = np.asarray(resistance_shunt, np.float64)
-    Rs = np.asarray(resistance_series, np.float64)
-    nNsVth = np.asarray(nNsVth, np.float64)
-    I = np.asarray(current, np.float64)
-    I0 = np.asarray(saturation_current, np.float64)
-    IL = np.asarray(photocurrent, np.float64)
+    # Record if inputs were all scalar
+    output_is_scalar = all(map(np.isscalar,
+                               [resistance_shunt, resistance_series, nNsVth,
+                                current, saturation_current, photocurrent]))
+
+    # 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
+    Rsh, Rs, a, I, I0, IL = \
+        np.broadcast_arrays(resistance_shunt, resistance_series, nNsVth,
+                            current, saturation_current, photocurrent)
 
     # This transforms any ideal Rsh=np.inf into Gsh=0., which is generally
     #  more numerically stable
     Gsh = 1./Rsh
 
-    # Intitalize output V (including shape) by solving explicit model with
-    #   Gsh=0, multiplying by np.ones_like(Gsh) identity in order to also
-    #   capture shape of Gsh.
-    V = (nNsVth*np.log1p((IL - I)/I0) - I*Rs)*np.ones_like(Gsh)
-
-    # Record if inputs were all scalar
-    output_is_scalar = np.isscalar(V)
+    # Intitalize output V (I might not be float64)
+    V = np.full_like(I, np.nan, dtype=np.float64)
 
-    # Multiply by np.atleast_1d in order to convert scalars to arrays, because
-    #  we need to guarantee the ability to index arrays
-    V = np.atleast_1d(V)
+    # Determine indices where 0 < Gsh requires implicit model solution
+    idx_p = 0. < Gsh
 
-    # Expand Gsh input shape to match output V (if needed)
-    Gsh = Gsh*np.ones_like(V)
+    # Determine indices where 0 = Gsh allows explicit model solution
+    idx_z = 0. == Gsh
 
-    # Determine indices where Gsh>0 requires implicit model solution
-    idx = 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):
-        # Expand remaining inputs to accomodate common indexing
-        Rs = Rs*np.ones_like(V)
-        nNsVth = nNsVth*np.ones_like(V)
-        I = I*np.ones_like(V)
-        I0 = I0*np.ones_like(V)
-        IL = IL*np.ones_like(V)
-
-        # LambertW argument, argW cannot be float128
-        argW = I0[idx] / (Gsh[idx]*nNsVth[idx]) * \
-            np.exp((-I[idx] + IL[idx] + I0[idx]) / (Gsh[idx]*nNsVth[idx]))
+    if np.any(idx_p):
+        # LambertW argument, cannot be float128, may overflow to np.inf
+        argW = I0[idx_p] / (Gsh[idx_p]*a[idx_p]) * \
+            np.exp((-I[idx_p] + IL[idx_p] + I0[idx_p]) /
+                   (Gsh[idx_p]*a[idx_p]))
 
         # lambertw typically returns complex value with zero imaginary part
         lambertwterm = lambertw(argW).real
 
-        # Record indices where LambertW input overflowed output
-        idx_w = np.logical_not(np.isfinite(lambertwterm))
+        # Record indices where lambertw input overflowed output
+        idx_inf = np.logical_not(np.isfinite(lambertwterm))
 
         # Only re-compute LambertW if it overflowed
-        if np.any(idx_w):
+        if np.any(idx_inf):
             # Calculate using log(argW) in case argW is really big
-            logargW = (np.log(I0[idx]) - np.log(Gsh[idx]) -
-                       np.log(nNsVth[idx]) +
-                       (-I[idx] + IL[idx] + I0[idx]) /
-                       (Gsh[idx] * nNsVth[idx]))[idx_w]
+            logargW = (np.log(I0[idx_p]) - np.log(Gsh[idx_p]) -
+                       np.log(a[idx_p]) +
+                       (-I[idx_p] + IL[idx_p] + I0[idx_p]) /
+                       (Gsh[idx_p] * a[idx_p]))[idx_inf]
 
             # Three iterations of Newton-Raphson method to solve
             #  w+log(w)=logargW. The initial guess is w=logargW. Where direct
@@ -1945,14 +1938,14 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
             #  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_w] = w
+                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 - nNsVth*lambertwterm + I0*Rsh
-        # Recasted in terms of Gsh=1/Rsh for better numerical stability.
-        V[idx] = (IL[idx] + I0[idx] - I[idx])/Gsh[idx] - I[idx]*Rs[idx] - \
-            nNsVth[idx]*lambertwterm
+        #  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[idx_p] + I0[idx_p] - I[idx_p])/Gsh[idx_p] - \
+            I[idx_p]*Rs[idx_p] - a[idx_p]*lambertwterm
 
     if output_is_scalar:
         return np.asscalar(V)
@@ -1969,10 +1962,11 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
      Jain and Kapoor 2004 [1].
     The solution is per Eq 2 of [1] except when resistance_series=0,
      in which case the explict solution for current is used.
-    Inputs to this function can include scalars and pandas.Series, but it
-     always outputs a float64 numpy.ndarray regardless of input type(s).
     Ideal device parameters are specified by resistance_shunt=np.inf and
      resistance_series=0.
+    Inputs to this function can include scalars and pandas.Series, but it is
+     the caller's responsibility to ensure that the arguments are all float64
+     and within the proper ranges.
 
     Parameters
     ----------
@@ -2024,51 +2018,42 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
     except ImportError:
         raise ImportError('This function requires scipy')
 
-    # asarray turns Series into arrays so that we don't have to worry
-    #  about multidimensional broadcasting failing
-    # Note that lambertw function doesn't support float128
-    Rsh = np.asarray(resistance_shunt, np.float64)
-    Rs = np.asarray(resistance_series, np.float64)
-    nNsVth = np.asarray(nNsVth, np.float64)
-    V = np.asarray(voltage, np.float64)
-    I0 = np.asarray(saturation_current, np.float64)
-    IL = np.asarray(photocurrent, np.float64)
+    # Record if inputs were all scalar
+    output_is_scalar = all(map(np.isscalar,
+                               [resistance_shunt, resistance_series, nNsVth,
+                                voltage, saturation_current, photocurrent]))
+
+    # 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
+    Rsh, Rs, a, V, I0, IL = \
+        np.broadcast_arrays(resistance_shunt, resistance_series, nNsVth,
+                            voltage, saturation_current, photocurrent)
 
     # This transforms any ideal Rsh=np.inf into Gsh=0., which is generally
     #  more numerically stable
     Gsh = 1./Rsh
 
-    # Intitalize output I (including shape) by solving explicit model with
-    #   Rs=0, multiplying by np.ones_like(Rs) identity in order to also
-    #   capture shape of Rs.
-    I = (IL - I0*np.expm1(V/nNsVth) - Gsh*V)*np.ones_like(Rs)
-
-    # Record if inputs were all scalar
-    output_is_scalar = np.isscalar(I)
+    # Intitalize output I (V might not be float64)
+    I = np.full_like(V, np.nan, dtype=np.float64)
 
-    # Multiply by np.atleast_1d in order to convert scalars to arrays, because
-    #  we need to guarantee the ability to index arrays
-    I = np.atleast_1d(I)
+    # Determine indices where 0 < Rs requires implicit model solution
+    idx_p = 0. < Rs
 
-    # Expand Rs input shape to match output I (if needed)
-    Rs = Rs*np.ones_like(I)
+    # Determine indices where 0 = Rs allows explicit model solution
+    idx_z = 0. == Rs
 
-    # Determine indices where Rs>0 requires implicit model solution
-    idx = 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
-    if np.any(idx):
-        # Expand remaining inputs to accomodate common indexing
-        Gsh = Gsh*np.ones_like(I)
-        nNsVth = nNsVth*np.ones_like(I)
-        V = V*np.ones_like(I)
-        I0 = I0*np.ones_like(I)
-        IL = IL*np.ones_like(I)
-
-        # LambertW argument, argW cannot be float128
-        argW = Rs[idx]*I0[idx]/(nNsVth[idx]*(Rs[idx]*Gsh[idx] + 1.)) * \
-            np.exp((Rs[idx]*(IL[idx] + I0[idx]) + V[idx]) /
-                   (nNsVth[idx]*(Rs[idx]*Gsh[idx] + 1.)))
+    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
         lambertwterm = lambertw(argW).real
@@ -2076,9 +2061,9 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
         # Eqn. 2 in Jain and Kapoor, 2004
         #  I = -V/(Rs + Rsh) - (nNsVth/Rs)*lambertwterm + \
         #      Rsh*(IL + I0)/(Rs + Rsh)
-        # Recasted in terms of Gsh=1/Rsh for better numerical stability.
-        I[idx] = (IL[idx] + I0[idx] - V[idx]*Gsh[idx]) / \
-            (Rs[idx]*Gsh[idx] + 1.) - (nNsVth[idx]/Rs[idx])*lambertwterm
+        # 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 np.asscalar(I)
diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index 2259b30623..c7e7345ca4 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -370,7 +370,7 @@ def test_PVSystem_calcparams_desoto(cec_module_params):
      'I': np.array(3.),
      'I0': np.array(6.e-7),
      'IL': np.array(7.),
-     'V_expected': 7.5049875193450521
+     'V_expected': np.array(7.5049875193450521)
     },
     {  # Can handle all rank-1 singleton array inputs
      'Rsh': np.array([20.]),
@@ -509,7 +509,7 @@ def test_v_from_i(fixture_v_from_i):
       'V': np.array(40.),
       'I0': np.array(6.e-7),
       'IL': np.array(7.),
-      'I_expected': -299.746389916
+      'I_expected': np.array(-299.746389916)
     },
     {  # Can handle all rank-1 singleton array inputs
       'Rsh': np.array([20.]),

From 6c410d27ec1fa497e0d716bc2546d8907f951fd5 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Tue, 12 Sep 2017 02:06:43 -0600
Subject: [PATCH 21/22] One more simplification

---
 pvlib/pvsystem.py | 23 +++++++++++------------
 1 file changed, 11 insertions(+), 12 deletions(-)

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index cddc33a3e2..4cc3d90b4e 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1889,14 +1889,12 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
     # 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
-    Rsh, Rs, a, I, I0, IL = \
-        np.broadcast_arrays(resistance_shunt, resistance_series, nNsVth,
+    # This transforms Gsh=1/Rsh, including ideal Rsh=np.inf into Gsh=0., which
+    #  is generally more numerically stable
+    Gsh, Rs, a, I, I0, IL = \
+        np.broadcast_arrays(1./resistance_shunt, resistance_series, nNsVth,
                             current, saturation_current, photocurrent)
 
-    # This transforms any ideal Rsh=np.inf into Gsh=0., which is generally
-    #  more numerically stable
-    Gsh = 1./Rsh
-
     # Intitalize output V (I might not be float64)
     V = np.full_like(I, np.nan, dtype=np.float64)
 
@@ -1919,6 +1917,7 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
                    (Gsh[idx_p]*a[idx_p]))
 
         # 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
@@ -2026,14 +2025,12 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
     # 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
-    Rsh, Rs, a, V, I0, IL = \
-        np.broadcast_arrays(resistance_shunt, resistance_series, nNsVth,
+    # This transforms Gsh=1/Rsh, including ideal Rsh=np.inf into Gsh=0., which
+    #  is generally more numerically stable
+    Gsh, Rs, a, V, I0, IL = \
+        np.broadcast_arrays(1./resistance_shunt, resistance_series, nNsVth,
                             voltage, saturation_current, photocurrent)
 
-    # This transforms any ideal Rsh=np.inf into Gsh=0., which is generally
-    #  more numerically stable
-    Gsh = 1./Rsh
-
     # Intitalize output I (V might not be float64)
     I = np.full_like(V, np.nan, dtype=np.float64)
 
@@ -2049,6 +2046,7 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
             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.)) * \
@@ -2056,6 +2054,7 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
                    (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

From 19af023021494881a8dc7aa24dfe76f6a2b2a31b Mon Sep 17 00:00:00 2001
From: Mark Campanelli <campy@instreamdesigns.com>
Date: Tue, 12 Sep 2017 07:21:35 -0600
Subject: [PATCH 22/22] Better use of broadcast_arrays

---
 pvlib/pvsystem.py | 19 +++++++++++--------
 1 file changed, 11 insertions(+), 8 deletions(-)

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 4cc3d90b4e..b19cce54f6 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1886,13 +1886,15 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
                                [resistance_shunt, resistance_series, nNsVth,
                                 current, saturation_current, photocurrent]))
 
+    # 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
-    # This transforms Gsh=1/Rsh, including ideal Rsh=np.inf into Gsh=0., which
-    #  is generally more numerically stable
     Gsh, Rs, a, I, I0, IL = \
-        np.broadcast_arrays(1./resistance_shunt, resistance_series, nNsVth,
+        np.broadcast_arrays(conductance_shunt, resistance_series, nNsVth,
                             current, saturation_current, photocurrent)
 
     # Intitalize output V (I might not be float64)
@@ -2022,13 +2024,15 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
                                [resistance_shunt, resistance_series, nNsVth,
                                 voltage, saturation_current, photocurrent]))
 
+    # 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
-    # This transforms Gsh=1/Rsh, including ideal Rsh=np.inf into Gsh=0., which
-    #  is generally more numerically stable
     Gsh, Rs, a, V, I0, IL = \
-        np.broadcast_arrays(1./resistance_shunt, resistance_series, nNsVth,
+        np.broadcast_arrays(conductance_shunt, resistance_series, nNsVth,
                             voltage, saturation_current, photocurrent)
 
     # Intitalize output I (V might not be float64)
@@ -2058,8 +2062,7 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
         lambertwterm = lambertw(argW).real
 
         # Eqn. 2 in Jain and Kapoor, 2004
-        #  I = -V/(Rs + Rsh) - (nNsVth/Rs)*lambertwterm + \
-        #      Rsh*(IL + I0)/(Rs + Rsh)
+        #  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