From dc1be1463cc950cc8e6c495dba51c9110b7504f4 Mon Sep 17 00:00:00 2001
From: Anton Driesse <anton.driesse@pvperformancelabs.com>
Date: Fri, 18 Oct 2019 22:32:26 +0200
Subject: [PATCH 01/11] Functional new function.

---
 pvlib/pvsystem.py | 73 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 73 insertions(+)

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index bebe2f1263..86786beb85 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1202,6 +1202,79 @@ def iam_martin_ruiz(aoi, a_r=0.16):
     return iam
 
 
+def iam_martin_ruiz_diffuse(slope, a_r=0.16, c1=0.4244, c2=None):
+    '''
+    Determine the incidence angle modifier (iam) for diffuse irradiance using
+    the Martin and Ruiz incident angle model.
+
+    Parameters
+    ----------
+    slope : numeric, degrees
+        The slope or tilt angle of of the module, where 0 is horizontal.
+
+    a_r : numeric
+        The angular losses coefficient described in equation 3 of [1].
+        This is an empirical dimensionless parameter. Values of a_r are
+        generally on the order of 0.08 to 0.25 for flat-plate PV modules.
+        a_r must be a positive numeric scalar.
+
+    Returns
+    -------
+    iam_sky : numeric
+        The incident angle modifier(s) for sky diffuse
+
+    iam_gnd : numeric
+        The incident angle modifier(s) for ground reflected diffuse
+
+    Notes
+    -----
+    Sky and ground modifiers are complementary, and equal for vertical surfaces
+    '''
+    # Contributed by Anton Driesse (@adriesse), PV Performance Labs. Oct., 2019
+
+    if isinstance(slope, pd.Series):
+        out_index = slope.index
+    else:
+        out_index = None
+
+    slope = np.asanyarray(slope)
+
+    # undo possible surface rotations
+    slope = (slope + 180) % 360 - 180
+
+    # flip backward tilts forward
+    slope = np.abs(slope)
+
+    # avoid undefined results for horizontal or upside-down surfaces
+    slope = np.clip(slope, 0.001, 180 - 0.001)
+
+    if np.any(np.less_equal(a_r, 0)):
+        raise RuntimeError("The parameter 'a_r' cannot be zero or negative.")
+
+    if c2 is None:
+        # See IEC 61853-3; valid for 0.16 <= a_r <= 0.18
+        c2 = 0.5 * a_r - 0.154
+        if np.any(np.less(a_r, 0.15)) or np.any(np.greater(a_r, 0.19)):
+            warnings.warn('A calculated value for c2 lies outside the range'
+                      'reported in the literature. See docstring for details.')
+
+    from numpy import pi, radians, sin, cos, exp
+
+    beta = radians(slope)
+
+    trig_term_sky = sin(beta) + (pi - beta - sin(beta)) / (1 + cos(beta))
+    trig_term_gnd = sin(beta) +      (beta - sin(beta)) / (1 - cos(beta))
+
+    iam_sky = 1 - exp(-(c1 + c2 * trig_term_sky ) * trig_term_sky / a_r)
+    iam_gnd = 1 - exp(-(c1 + c2 * trig_term_gnd ) * trig_term_gnd / a_r)
+
+    if out_index is not None:
+        iam_sky = pd.Series(iam_sky, index=out_index, name='iam_sky')
+        iam_gnd = pd.Series(iam_gnd, index=out_index, name='iam_ground')
+
+    return iam_sky, iam_gnd
+
+
 def iam_interp(aoi, theta_ref, iam_ref, method='linear', normalize=True):
     '''
     Determine the incidence angle modifier (iam) by interpolating a set of

From ae3eea9e7591f05b7df9a133c0c1dd2c34f9c43d Mon Sep 17 00:00:00 2001
From: Anton Driesse <anton.driesse@pvperformancelabs.com>
Date: Fri, 18 Oct 2019 22:39:36 +0200
Subject: [PATCH 02/11] Undo

---
 pvlib/pvsystem.py | 73 -----------------------------------------------
 1 file changed, 73 deletions(-)

diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 86786beb85..bebe2f1263 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1202,79 +1202,6 @@ def iam_martin_ruiz(aoi, a_r=0.16):
     return iam
 
 
-def iam_martin_ruiz_diffuse(slope, a_r=0.16, c1=0.4244, c2=None):
-    '''
-    Determine the incidence angle modifier (iam) for diffuse irradiance using
-    the Martin and Ruiz incident angle model.
-
-    Parameters
-    ----------
-    slope : numeric, degrees
-        The slope or tilt angle of of the module, where 0 is horizontal.
-
-    a_r : numeric
-        The angular losses coefficient described in equation 3 of [1].
-        This is an empirical dimensionless parameter. Values of a_r are
-        generally on the order of 0.08 to 0.25 for flat-plate PV modules.
-        a_r must be a positive numeric scalar.
-
-    Returns
-    -------
-    iam_sky : numeric
-        The incident angle modifier(s) for sky diffuse
-
-    iam_gnd : numeric
-        The incident angle modifier(s) for ground reflected diffuse
-
-    Notes
-    -----
-    Sky and ground modifiers are complementary, and equal for vertical surfaces
-    '''
-    # Contributed by Anton Driesse (@adriesse), PV Performance Labs. Oct., 2019
-
-    if isinstance(slope, pd.Series):
-        out_index = slope.index
-    else:
-        out_index = None
-
-    slope = np.asanyarray(slope)
-
-    # undo possible surface rotations
-    slope = (slope + 180) % 360 - 180
-
-    # flip backward tilts forward
-    slope = np.abs(slope)
-
-    # avoid undefined results for horizontal or upside-down surfaces
-    slope = np.clip(slope, 0.001, 180 - 0.001)
-
-    if np.any(np.less_equal(a_r, 0)):
-        raise RuntimeError("The parameter 'a_r' cannot be zero or negative.")
-
-    if c2 is None:
-        # See IEC 61853-3; valid for 0.16 <= a_r <= 0.18
-        c2 = 0.5 * a_r - 0.154
-        if np.any(np.less(a_r, 0.15)) or np.any(np.greater(a_r, 0.19)):
-            warnings.warn('A calculated value for c2 lies outside the range'
-                      'reported in the literature. See docstring for details.')
-
-    from numpy import pi, radians, sin, cos, exp
-
-    beta = radians(slope)
-
-    trig_term_sky = sin(beta) + (pi - beta - sin(beta)) / (1 + cos(beta))
-    trig_term_gnd = sin(beta) +      (beta - sin(beta)) / (1 - cos(beta))
-
-    iam_sky = 1 - exp(-(c1 + c2 * trig_term_sky ) * trig_term_sky / a_r)
-    iam_gnd = 1 - exp(-(c1 + c2 * trig_term_gnd ) * trig_term_gnd / a_r)
-
-    if out_index is not None:
-        iam_sky = pd.Series(iam_sky, index=out_index, name='iam_sky')
-        iam_gnd = pd.Series(iam_gnd, index=out_index, name='iam_ground')
-
-    return iam_sky, iam_gnd
-
-
 def iam_interp(aoi, theta_ref, iam_ref, method='linear', normalize=True):
     '''
     Determine the incidence angle modifier (iam) by interpolating a set of

From fa319f1c7eb0292ec4eb2a2bf00d913da0acb36e Mon Sep 17 00:00:00 2001
From: Anton Driesse <anton.driesse@pvperformancelabs.com>
Date: Sat, 19 Oct 2019 15:59:28 +0200
Subject: [PATCH 03/11] Working function; partial docstring; no tests.

---
 pvlib/iam.py | 98 +++++++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 97 insertions(+), 1 deletion(-)

diff --git a/pvlib/iam.py b/pvlib/iam.py
index 3d869f298d..6350a9e7c7 100644
--- a/pvlib/iam.py
+++ b/pvlib/iam.py
@@ -11,7 +11,7 @@
 import numpy as np
 import pandas as pd
 from pvlib.tools import cosd, sind, tand, asind
-
+import warnings
 
 # a dict of required parameter names for each IAM model
 # keys are the function names for the IAM models
@@ -274,6 +274,102 @@ def martin_ruiz(aoi, a_r=0.16):
     return iam
 
 
+def martin_ruiz_diffuse(slope, a_r=0.16, c1=0.4244, c2=None):
+    '''
+    Determine the incidence angle modifiers (iam) for diffuse irradiance using
+    the Martin and Ruiz incident angle model.
+
+    Parameters
+    ----------
+    slope : numeric, degrees
+        The slope or tilt angle of of the module, where 0 is horizontal.
+
+    a_r : numeric
+        The angular losses coefficient described in equation 3 of [1].
+        This is an empirical dimensionless parameter. Values of a_r are
+        generally on the order of 0.08 to 0.25 for flat-plate PV modules.
+        a_r must be a positive numeric scalar.
+
+    c1 :
+
+    c2 :
+
+
+    Returns
+    -------
+    iam_sky : numeric
+        The incident angle modifier(s) for sky diffuse
+
+    iam_gnd : numeric
+        The incident angle modifier(s) for ground reflected diffuse
+
+    Notes
+    -----
+    Sky and ground modifiers are complementary, and equal for vertical surfaces
+
+    References
+    ----------
+    [1] N. Martin and J. M. Ruiz, "Calculation of the PV modules angular
+    losses under field conditions by means of an analytical model", Solar
+    Energy Materials & Solar Cells, vol. 70, pp. 25-38, 2001.
+
+    [2] N. Martin and J. M. Ruiz, "Corrigendum to 'Calculation of the PV
+    modules angular losses under field conditions by means of an
+    analytical model'", Solar Energy Materials & Solar Cells, vol. 110,
+    pp. 154, 2013.
+
+    See Also
+    --------
+    iam.physical
+    iam.ashrae
+    iam.interp
+    iam.sapm
+    '''
+    # Contributed by Anton Driesse (@adriesse), PV Performance Labs. Oct. 2019
+
+    if isinstance(slope, pd.Series):
+        out_index = slope.index
+    else:
+        out_index = None
+
+    slope = np.asanyarray(slope)
+
+    # undo possible surface rotations
+    slope = (slope + 180) % 360 - 180
+
+    # flip backward tilts forward
+    slope = np.abs(slope)
+
+    # avoid undefined results for horizontal or upside-down surfaces
+    slope = np.clip(slope, 0.001, 180 - 0.001)
+
+    if np.any(np.less_equal(a_r, 0)):
+        raise RuntimeError("The parameter 'a_r' cannot be zero or negative.")
+
+    if c2 is None:
+        # See IEC 61853-3; valid for 0.16 <= a_r <= 0.18; to be explained above
+        c2 = 0.5 * a_r - 0.154
+        if np.any(np.less(a_r, 0.15)) or np.any(np.greater(a_r, 0.19)):
+            warnings.warn('A calculated value for c2 lies outside the range'
+                      'reported in the literature. See docstring for details.')
+
+    from numpy import pi, radians, sin, cos, exp
+
+    beta = radians(slope)
+
+    trig_term_sky = sin(beta) + (pi - beta - sin(beta)) / (1 + cos(beta))
+    trig_term_gnd = sin(beta) +      (beta - sin(beta)) / (1 - cos(beta))
+
+    iam_sky = 1 - exp(-(c1 + c2 * trig_term_sky ) * trig_term_sky / a_r)
+    iam_gnd = 1 - exp(-(c1 + c2 * trig_term_gnd ) * trig_term_gnd / a_r)
+
+    if out_index is not None:
+        iam_sky = pd.Series(iam_sky, index=out_index, name='iam_sky')
+        iam_gnd = pd.Series(iam_gnd, index=out_index, name='iam_ground')
+
+    return iam_sky, iam_gnd
+
+
 def interp(aoi, theta_ref, iam_ref, method='linear', normalize=True):
     r'''
     Determine the incidence angle modifier (IAM) by interpolating a set of

From d65d953afb014d25a7a6aa13ee10f0875e6f69f3 Mon Sep 17 00:00:00 2001
From: Anton Driesse <anton.driesse@pvperformancelabs.com>
Date: Sat, 19 Oct 2019 21:44:24 +0200
Subject: [PATCH 04/11] Complete docstring; simplify a bit; start tests.

---
 pvlib/iam.py           | 31 +++++++++++++++++--------------
 pvlib/test/test_iam.py | 14 ++++++++++++++
 2 files changed, 31 insertions(+), 14 deletions(-)

diff --git a/pvlib/iam.py b/pvlib/iam.py
index 6350a9e7c7..c070ddc872 100644
--- a/pvlib/iam.py
+++ b/pvlib/iam.py
@@ -11,7 +11,6 @@
 import numpy as np
 import pandas as pd
 from pvlib.tools import cosd, sind, tand, asind
-import warnings
 
 # a dict of required parameter names for each IAM model
 # keys are the function names for the IAM models
@@ -220,8 +219,8 @@ def martin_ruiz(aoi, a_r=0.16):
     -----
     `martin_ruiz` calculates the incidence angle modifier (IAM) as described in
     [1]. The information required is the incident angle (AOI) and the angular
-    losses coefficient (a_r). Note that [1] has a corrigendum [2] which makes
-    the document much simpler to understand.
+    losses coefficient (a_r). Note that [1] has a corrigendum [2] which
+    clarifies a mix-up of 'alpha's and 'a's in the former.
 
     The incident angle modifier is defined as
 
@@ -249,6 +248,7 @@ def martin_ruiz(aoi, a_r=0.16):
 
     See Also
     --------
+    iam.martin_ruiz_diffuse
     iam.physical
     iam.ashrae
     iam.interp
@@ -290,10 +290,11 @@ def martin_ruiz_diffuse(slope, a_r=0.16, c1=0.4244, c2=None):
         generally on the order of 0.08 to 0.25 for flat-plate PV modules.
         a_r must be a positive numeric scalar.
 
-    c1 :
-
-    c2 :
-
+    c1, c2 : numeric scalar
+        Two fitting parameters for the expressions that approximate the
+        integral of diffuse irradiance coming from different directions.
+        c1 is contant at 4 / 3 / pi (0.4244) and c2 varies with a_r.  The
+        calculation of c2 from a_r is according to [3].
 
     Returns
     -------
@@ -305,7 +306,9 @@ def martin_ruiz_diffuse(slope, a_r=0.16, c1=0.4244, c2=None):
 
     Notes
     -----
-    Sky and ground modifiers are complementary, and equal for vertical surfaces
+    Sky and ground modifiers are complementary, so iam_sky for slope = 30 is
+    equal to iam_gnd for slope = 180 - 30.  For vertical surfaces, slope = 90,
+    the two factors are equal.
 
     References
     ----------
@@ -318,8 +321,12 @@ def martin_ruiz_diffuse(slope, a_r=0.16, c1=0.4244, c2=None):
     analytical model'", Solar Energy Materials & Solar Cells, vol. 110,
     pp. 154, 2013.
 
+    [3] "IEC 61853-3 Photovoltaic (PV) module performance testing and energy
+    rating - Part 3: Energy rating of PV modules". IEC, Geneva, 2018.
+
     See Also
     --------
+    iam.martin_ruiz
     iam.physical
     iam.ashrae
     iam.interp
@@ -347,11 +354,7 @@ def martin_ruiz_diffuse(slope, a_r=0.16, c1=0.4244, c2=None):
         raise RuntimeError("The parameter 'a_r' cannot be zero or negative.")
 
     if c2 is None:
-        # See IEC 61853-3; valid for 0.16 <= a_r <= 0.18; to be explained above
         c2 = 0.5 * a_r - 0.154
-        if np.any(np.less(a_r, 0.15)) or np.any(np.greater(a_r, 0.19)):
-            warnings.warn('A calculated value for c2 lies outside the range'
-                      'reported in the literature. See docstring for details.')
 
     from numpy import pi, radians, sin, cos, exp
 
@@ -360,8 +363,8 @@ def martin_ruiz_diffuse(slope, a_r=0.16, c1=0.4244, c2=None):
     trig_term_sky = sin(beta) + (pi - beta - sin(beta)) / (1 + cos(beta))
     trig_term_gnd = sin(beta) +      (beta - sin(beta)) / (1 - cos(beta))
 
-    iam_sky = 1 - exp(-(c1 + c2 * trig_term_sky ) * trig_term_sky / a_r)
-    iam_gnd = 1 - exp(-(c1 + c2 * trig_term_gnd ) * trig_term_gnd / a_r)
+    iam_sky = 1 - exp(-(c1 + c2 * trig_term_sky) * trig_term_sky / a_r)
+    iam_gnd = 1 - exp(-(c1 + c2 * trig_term_gnd) * trig_term_gnd / a_r)
 
     if out_index is not None:
         iam_sky = pd.Series(iam_sky, index=out_index, name='iam_sky')
diff --git a/pvlib/test/test_iam.py b/pvlib/test/test_iam.py
index 0c51b005a1..b50b321ff2 100644
--- a/pvlib/test/test_iam.py
+++ b/pvlib/test/test_iam.py
@@ -104,6 +104,20 @@ def test_martin_ruiz():
         _iam.martin_ruiz(0.0, a_r=0.0)
 
 
+def test_martin_ruiz_diffuse():
+
+    slope = 30.
+    a_r = 0.16
+    expected = (0.9549735, 0.7944426)
+
+    # will fail if default values change
+    iam = _iam.martin_ruiz_diffuse(slope)
+    assert_allclose(iam, expected)
+    # will fail if parameter names change
+    iam = _iam.martin_ruiz_diffuse(slope=slope, a_r=a_r)
+    assert_allclose(iam, expected)
+
+
 @requires_scipy
 def test_iam_interp():
 

From 143aecde1d03f80343a40b88912ab32fee243b93 Mon Sep 17 00:00:00 2001
From: Anton Driesse <anton.driesse@pvperformancelabs.com>
Date: Sun, 20 Oct 2019 23:36:34 +0200
Subject: [PATCH 05/11] Getting better all the time!

---
 pvlib/iam.py           | 26 ++++++++++++++++----------
 pvlib/test/test_iam.py | 28 ++++++++++++++++++++++++++++
 2 files changed, 44 insertions(+), 10 deletions(-)

diff --git a/pvlib/iam.py b/pvlib/iam.py
index c070ddc872..0dba21fe77 100644
--- a/pvlib/iam.py
+++ b/pvlib/iam.py
@@ -341,14 +341,16 @@ def martin_ruiz_diffuse(slope, a_r=0.16, c1=0.4244, c2=None):
 
     slope = np.asanyarray(slope)
 
-    # undo possible surface rotations
-    slope = (slope + 180) % 360 - 180
+    with np.errstate(invalid='ignore'):
+        # undo possible surface rotations
+        slope = (slope + 180) % 360 - 180
 
-    # flip backward tilts forward
-    slope = np.abs(slope)
+        # flip backward tilts forward
+        slope = np.abs(slope)
 
-    # avoid undefined results for horizontal or upside-down surfaces
-    slope = np.clip(slope, 0.001, 180 - 0.001)
+        # avoid undefined results for horizontal or upside-down surfaces
+        small_angle = 1e-06
+        slope = np.clip(slope, small_angle, 180 - small_angle)
 
     if np.any(np.less_equal(a_r, 0)):
         raise RuntimeError("The parameter 'a_r' cannot be zero or negative.")
@@ -360,11 +362,15 @@ def martin_ruiz_diffuse(slope, a_r=0.16, c1=0.4244, c2=None):
 
     beta = radians(slope)
 
-    trig_term_sky = sin(beta) + (pi - beta - sin(beta)) / (1 + cos(beta))
-    trig_term_gnd = sin(beta) +      (beta - sin(beta)) / (1 - cos(beta))
+    with np.errstate(invalid='ignore'):
+        # because sin(pi) isn't exactly zero
+        sin_beta = np.where(slope < 90, sin(beta), sin(pi - beta))
+
+        trig_term_sky = sin_beta + (pi - beta - sin_beta) / (1 + cos(beta))
+        trig_term_gnd = sin_beta +      (beta - sin_beta) / (1 - cos(beta))
 
-    iam_sky = 1 - exp(-(c1 + c2 * trig_term_sky) * trig_term_sky / a_r)
-    iam_gnd = 1 - exp(-(c1 + c2 * trig_term_gnd) * trig_term_gnd / a_r)
+        iam_sky = 1 - exp(-(c1 + c2 * trig_term_sky) * trig_term_sky / a_r)
+        iam_gnd = 1 - exp(-(c1 + c2 * trig_term_gnd) * trig_term_gnd / a_r)
 
     if out_index is not None:
         iam_sky = pd.Series(iam_sky, index=out_index, name='iam_sky')
diff --git a/pvlib/test/test_iam.py b/pvlib/test/test_iam.py
index b50b321ff2..ec15214b9b 100644
--- a/pvlib/test/test_iam.py
+++ b/pvlib/test/test_iam.py
@@ -117,6 +117,34 @@ def test_martin_ruiz_diffuse():
     iam = _iam.martin_ruiz_diffuse(slope=slope, a_r=a_r)
     assert_allclose(iam, expected)
 
+    a_r = 0.18
+    slope = [0, 30, 90, 120, -30, np.nan, np.inf]
+    expected_sky = [0.9407678, 0.9452250, 0.9407678, 0.9055541, 0.9452250,
+                    np.nan, np.nan]
+    expected_gnd = [0.0000000, 0.7610849, 0.9407678, 0.9483508, 0.7610849,
+                    np.nan, np.nan]
+
+    # check various inputs as list
+    iam = _iam.martin_ruiz_diffuse(slope, a_r)
+    assert_allclose(iam[0], expected_sky, atol=1e-7, equal_nan=True)
+    assert_allclose(iam[1], expected_gnd, atol=1e-7, equal_nan=True)
+
+    # check various inputs as array
+    iam = _iam.martin_ruiz_diffuse(np.array(slope), a_r)
+    assert_allclose(iam[0], expected_sky, atol=1e-7, equal_nan=True)
+    assert_allclose(iam[1], expected_gnd, atol=1e-7, equal_nan=True)
+
+    # check various inputs as Series
+    slope = pd.Series(slope)
+    expected_sky = pd.Series(expected_sky, name='iam_sky')
+    expected_gnd = pd.Series(expected_gnd, name='iam_ground')
+    iam = _iam.martin_ruiz_diffuse(slope, a_r)
+    assert_series_equal(iam[0], expected_sky)
+    assert_series_equal(iam[1], expected_gnd)
+
+    # check exception clause
+    with pytest.raises(RuntimeError):
+        _iam.martin_ruiz_diffuse(0.0, a_r=0.0)
 
 @requires_scipy
 def test_iam_interp():

From 1624484d9d8e6769b8ca9895272ba226a6e79665 Mon Sep 17 00:00:00 2001
From: Anton Driesse <anton.driesse@pvperformancelabs.com>
Date: Mon, 21 Oct 2019 13:53:33 +0200
Subject: [PATCH 06/11] Adjust some comments and update whatsnew.

---
 docs/sphinx/source/whatsnew/v0.7.0.rst |  2 ++
 pvlib/iam.py                           | 16 ++++++++--------
 2 files changed, 10 insertions(+), 8 deletions(-)

diff --git a/docs/sphinx/source/whatsnew/v0.7.0.rst b/docs/sphinx/source/whatsnew/v0.7.0.rst
index 5a90f11e64..21860ac7aa 100644
--- a/docs/sphinx/source/whatsnew/v0.7.0.rst
+++ b/docs/sphinx/source/whatsnew/v0.7.0.rst
@@ -112,6 +112,8 @@ Enhancements
 * Add :py:func:`~pvlib.ivtools.fit_sdm_cec_sam`, a wrapper for the CEC single
   diode model fitting function '6parsolve' from NREL's System Advisor Model.
 * Add `timeout` to :py:func:`pvlib.iotools.get_psm3`.
+* Created one new incidence angle modifier (IAM) function for diffuse irradiance:
+  :py:func:`pvlib.iam.martin_ruiz_diffuse`. (:issue:`751`)
 
 Bug fixes
 ~~~~~~~~~
diff --git a/pvlib/iam.py b/pvlib/iam.py
index 0dba21fe77..9224bb5e7a 100644
--- a/pvlib/iam.py
+++ b/pvlib/iam.py
@@ -276,13 +276,13 @@ def martin_ruiz(aoi, a_r=0.16):
 
 def martin_ruiz_diffuse(slope, a_r=0.16, c1=0.4244, c2=None):
     '''
-    Determine the incidence angle modifiers (iam) for diffuse irradiance using
-    the Martin and Ruiz incident angle model.
+    Determine the incidence angle modifiers (iam) for diffuse sky and
+    ground-reflected irradiance using the Martin and Ruiz incident angle model.
 
     Parameters
     ----------
     slope : numeric, degrees
-        The slope or tilt angle of of the module, where 0 is horizontal.
+        The slope or tilt angle of the module, where 0 is horizontal.
 
     a_r : numeric
         The angular losses coefficient described in equation 3 of [1].
@@ -301,14 +301,14 @@ def martin_ruiz_diffuse(slope, a_r=0.16, c1=0.4244, c2=None):
     iam_sky : numeric
         The incident angle modifier(s) for sky diffuse
 
-    iam_gnd : numeric
-        The incident angle modifier(s) for ground reflected diffuse
+    iam_ground : numeric
+        The incident angle modifier(s) for ground-reflected diffuse
 
     Notes
     -----
-    Sky and ground modifiers are complementary, so iam_sky for slope = 30 is
-    equal to iam_gnd for slope = 180 - 30.  For vertical surfaces, slope = 90,
-    the two factors are equal.
+    Sky and ground modifiers are complementary: iam_sky for slope = 30 is
+    equal to iam_ground for slope = 180 - 30.  For vertical surfaces,
+    slope = 90, the two factors are equal.
 
     References
     ----------

From c5f907c15d94263caea4399d819a75ca1d65c18e Mon Sep 17 00:00:00 2001
From: Anton Driesse <anton.driesse@pvperformancelabs.com>
Date: Mon, 21 Oct 2019 21:35:59 +0200
Subject: [PATCH 07/11] Improvements as per review.

---
 docs/sphinx/source/api.rst |  1 +
 pvlib/iam.py               | 41 ++++++++++++++++++++------------------
 pvlib/test/test_iam.py     | 20 +++++++++----------
 3 files changed, 33 insertions(+), 29 deletions(-)

diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst
index c37bedf7df..bcd8359b61 100644
--- a/docs/sphinx/source/api.rst
+++ b/docs/sphinx/source/api.rst
@@ -220,6 +220,7 @@ Incident angle modifiers
    iam.physical
    iam.ashrae
    iam.martin_ruiz
+   iam.martin_ruiz_diffuse
    iam.sapm
    iam.interp
 
diff --git a/pvlib/iam.py b/pvlib/iam.py
index 9224bb5e7a..cd9d052850 100644
--- a/pvlib/iam.py
+++ b/pvlib/iam.py
@@ -262,7 +262,7 @@ def martin_ruiz(aoi, a_r=0.16):
     a_r = np.asanyarray(a_r)
 
     if np.any(np.less_equal(a_r, 0)):
-        raise RuntimeError("The parameter 'a_r' cannot be zero or negative.")
+        raise ValueError("The parameter 'a_r' cannot be zero or negative.")
 
     with np.errstate(invalid='ignore'):
         iam = (1 - np.exp(-cosd(aoi) / a_r)) / (1 - np.exp(-1 / a_r))
@@ -274,15 +274,17 @@ def martin_ruiz(aoi, a_r=0.16):
     return iam
 
 
-def martin_ruiz_diffuse(slope, a_r=0.16, c1=0.4244, c2=None):
+def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None):
     '''
     Determine the incidence angle modifiers (iam) for diffuse sky and
     ground-reflected irradiance using the Martin and Ruiz incident angle model.
 
     Parameters
     ----------
-    slope : numeric, degrees
-        The slope or tilt angle of the module, where 0 is horizontal.
+    surface_tilt: float or array-like, default 0
+        Surface tilt angles in decimal degrees.
+        The tilt angle is defined as degrees from horizontal
+        (e.g. surface facing up = 0, surface facing horizon = 90)
 
     a_r : numeric
         The angular losses coefficient described in equation 3 of [1].
@@ -293,22 +295,22 @@ def martin_ruiz_diffuse(slope, a_r=0.16, c1=0.4244, c2=None):
     c1, c2 : numeric scalar
         Two fitting parameters for the expressions that approximate the
         integral of diffuse irradiance coming from different directions.
-        c1 is contant at 4 / 3 / pi (0.4244) and c2 varies with a_r.  The
+        c1 is constant at 4 / 3 / pi (0.4244) and c2 varies with a_r.  The
         calculation of c2 from a_r is according to [3].
 
     Returns
     -------
     iam_sky : numeric
-        The incident angle modifier(s) for sky diffuse
+        The incident angle modifier for sky diffuse
 
     iam_ground : numeric
-        The incident angle modifier(s) for ground-reflected diffuse
+        The incident angle modifier for ground-reflected diffuse
 
     Notes
     -----
-    Sky and ground modifiers are complementary: iam_sky for slope = 30 is
-    equal to iam_ground for slope = 180 - 30.  For vertical surfaces,
-    slope = 90, the two factors are equal.
+    Sky and ground modifiers are complementary: iam_sky for tilt = 30 is
+    equal to iam_ground for tilt = 180 - 30.  For vertical surfaces,
+    tilt = 90, the two factors are equal.
 
     References
     ----------
@@ -334,37 +336,38 @@ def martin_ruiz_diffuse(slope, a_r=0.16, c1=0.4244, c2=None):
     '''
     # Contributed by Anton Driesse (@adriesse), PV Performance Labs. Oct. 2019
 
-    if isinstance(slope, pd.Series):
-        out_index = slope.index
+    if isinstance(surface_tilt, pd.Series):
+        out_index = surface_tilt.index
     else:
         out_index = None
 
-    slope = np.asanyarray(slope)
+    surface_tilt = np.asanyarray(surface_tilt)
 
     with np.errstate(invalid='ignore'):
         # undo possible surface rotations
-        slope = (slope + 180) % 360 - 180
+        surface_tilt = (surface_tilt + 180) % 360 - 180
 
         # flip backward tilts forward
-        slope = np.abs(slope)
+        surface_tilt = np.abs(surface_tilt)
 
         # avoid undefined results for horizontal or upside-down surfaces
         small_angle = 1e-06
-        slope = np.clip(slope, small_angle, 180 - small_angle)
+        surface_tilt = np.clip(surface_tilt, small_angle, 180 - small_angle)
 
     if np.any(np.less_equal(a_r, 0)):
-        raise RuntimeError("The parameter 'a_r' cannot be zero or negative.")
+        raise ValueError("The parameter 'a_r' cannot be zero or negative.")
 
     if c2 is None:
+        # This equation is from [3] Sect. 7.2
         c2 = 0.5 * a_r - 0.154
 
     from numpy import pi, radians, sin, cos, exp
 
-    beta = radians(slope)
+    beta = radians(surface_tilt)
 
     with np.errstate(invalid='ignore'):
         # because sin(pi) isn't exactly zero
-        sin_beta = np.where(slope < 90, sin(beta), sin(pi - beta))
+        sin_beta = np.where(surface_tilt < 90, sin(beta), sin(pi - beta))
 
         trig_term_sky = sin_beta + (pi - beta - sin_beta) / (1 + cos(beta))
         trig_term_gnd = sin_beta +      (beta - sin_beta) / (1 - cos(beta))
diff --git a/pvlib/test/test_iam.py b/pvlib/test/test_iam.py
index ec15214b9b..ef142c92e3 100644
--- a/pvlib/test/test_iam.py
+++ b/pvlib/test/test_iam.py
@@ -100,50 +100,50 @@ def test_martin_ruiz():
     assert_series_equal(iam, expected)
 
     # check exception clause
-    with pytest.raises(RuntimeError):
+    with pytest.raises(ValueError):
         _iam.martin_ruiz(0.0, a_r=0.0)
 
 
 def test_martin_ruiz_diffuse():
 
-    slope = 30.
+    surface_tilt = 30.
     a_r = 0.16
     expected = (0.9549735, 0.7944426)
 
     # will fail if default values change
-    iam = _iam.martin_ruiz_diffuse(slope)
+    iam = _iam.martin_ruiz_diffuse(surface_tilt)
     assert_allclose(iam, expected)
     # will fail if parameter names change
-    iam = _iam.martin_ruiz_diffuse(slope=slope, a_r=a_r)
+    iam = _iam.martin_ruiz_diffuse(surface_tilt=surface_tilt, a_r=a_r)
     assert_allclose(iam, expected)
 
     a_r = 0.18
-    slope = [0, 30, 90, 120, -30, np.nan, np.inf]
+    surface_tilt = [0, 30, 90, 120, -30, np.nan, np.inf]
     expected_sky = [0.9407678, 0.9452250, 0.9407678, 0.9055541, 0.9452250,
                     np.nan, np.nan]
     expected_gnd = [0.0000000, 0.7610849, 0.9407678, 0.9483508, 0.7610849,
                     np.nan, np.nan]
 
     # check various inputs as list
-    iam = _iam.martin_ruiz_diffuse(slope, a_r)
+    iam = _iam.martin_ruiz_diffuse(surface_tilt, a_r)
     assert_allclose(iam[0], expected_sky, atol=1e-7, equal_nan=True)
     assert_allclose(iam[1], expected_gnd, atol=1e-7, equal_nan=True)
 
     # check various inputs as array
-    iam = _iam.martin_ruiz_diffuse(np.array(slope), a_r)
+    iam = _iam.martin_ruiz_diffuse(np.array(surface_tilt), a_r)
     assert_allclose(iam[0], expected_sky, atol=1e-7, equal_nan=True)
     assert_allclose(iam[1], expected_gnd, atol=1e-7, equal_nan=True)
 
     # check various inputs as Series
-    slope = pd.Series(slope)
+    surface_tilt = pd.Series(surface_tilt)
     expected_sky = pd.Series(expected_sky, name='iam_sky')
     expected_gnd = pd.Series(expected_gnd, name='iam_ground')
-    iam = _iam.martin_ruiz_diffuse(slope, a_r)
+    iam = _iam.martin_ruiz_diffuse(surface_tilt, a_r)
     assert_series_equal(iam[0], expected_sky)
     assert_series_equal(iam[1], expected_gnd)
 
     # check exception clause
-    with pytest.raises(RuntimeError):
+    with pytest.raises(ValueError):
         _iam.martin_ruiz_diffuse(0.0, a_r=0.0)
 
 @requires_scipy

From a53b2538b5cdae1152765877111247ad481dff9b Mon Sep 17 00:00:00 2001
From: Anton Driesse <anton.driesse@pvperformancelabs.com>
Date: Tue, 22 Oct 2019 14:05:09 +0200
Subject: [PATCH 08/11] More changes as requested.

---
 pvlib/iam.py           | 29 +++++++++++++++++------------
 pvlib/test/test_iam.py | 11 +++++++++--
 2 files changed, 26 insertions(+), 14 deletions(-)

diff --git a/pvlib/iam.py b/pvlib/iam.py
index cd9d052850..e4e084858d 100644
--- a/pvlib/iam.py
+++ b/pvlib/iam.py
@@ -290,13 +290,18 @@ def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None):
         The angular losses coefficient described in equation 3 of [1].
         This is an empirical dimensionless parameter. Values of a_r are
         generally on the order of 0.08 to 0.25 for flat-plate PV modules.
-        a_r must be a positive numeric scalar.
+        a_r must be greater than zero.
 
-    c1, c2 : numeric scalar
-        Two fitting parameters for the expressions that approximate the
+    c1 : numeric scalar
+        First fitting parameter for the expressions that approximates the
         integral of diffuse irradiance coming from different directions.
-        c1 is constant at 4 / 3 / pi (0.4244) and c2 varies with a_r.  The
-        calculation of c2 from a_r is according to [3].
+        c1 is given as the constant 4 / 3 / pi (0.4244) in [1].
+
+    c2 : numeric scalar
+        Second fitting parameter for the expressions that approximates the
+        integral of diffuse irradiance coming from different directions.
+        If c2 is None, it will be calculated according to the linear
+        relationship given in [3].
 
     Returns
     -------
@@ -343,16 +348,16 @@ def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None):
 
     surface_tilt = np.asanyarray(surface_tilt)
 
+    # undo possible surface rotations
     with np.errstate(invalid='ignore'):
-        # undo possible surface rotations
         surface_tilt = (surface_tilt + 180) % 360 - 180
 
-        # flip backward tilts forward
-        surface_tilt = np.abs(surface_tilt)
+    # flip backward tilts forward
+    surface_tilt = np.abs(surface_tilt)
 
-        # avoid undefined results for horizontal or upside-down surfaces
-        small_angle = 1e-06
-        surface_tilt = np.clip(surface_tilt, small_angle, 180 - small_angle)
+    # avoid undefined results for horizontal or upside-down surfaces
+    small_angle = 1e-06
+    surface_tilt = np.clip(surface_tilt, small_angle, 180 - small_angle)
 
     if np.any(np.less_equal(a_r, 0)):
         raise ValueError("The parameter 'a_r' cannot be zero or negative.")
@@ -370,7 +375,7 @@ def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None):
         sin_beta = np.where(surface_tilt < 90, sin(beta), sin(pi - beta))
 
         trig_term_sky = sin_beta + (pi - beta - sin_beta) / (1 + cos(beta))
-        trig_term_gnd = sin_beta +      (beta - sin_beta) / (1 - cos(beta))
+        trig_term_gnd = sin_beta +      (beta - sin_beta) / (1 - cos(beta)) # noqa: E222
 
         iam_sky = 1 - exp(-(c1 + c2 * trig_term_sky) * trig_term_sky / a_r)
         iam_gnd = 1 - exp(-(c1 + c2 * trig_term_gnd) * trig_term_gnd / a_r)
diff --git a/pvlib/test/test_iam.py b/pvlib/test/test_iam.py
index ef142c92e3..29152c4b7f 100644
--- a/pvlib/test/test_iam.py
+++ b/pvlib/test/test_iam.py
@@ -77,6 +77,7 @@ def test_martin_ruiz():
     # will fail if default values change
     iam = _iam.martin_ruiz(aoi)
     assert_allclose(iam, expected)
+
     # will fail if parameter names change
     iam = _iam.martin_ruiz(aoi=aoi, a_r=a_r)
     assert_allclose(iam, expected)
@@ -99,7 +100,9 @@ def test_martin_ruiz():
     iam = _iam.martin_ruiz(aoi, a_r)
     assert_series_equal(iam, expected)
 
-    # check exception clause
+
+def test_martin_ruiz_exception():
+
     with pytest.raises(ValueError):
         _iam.martin_ruiz(0.0, a_r=0.0)
 
@@ -113,6 +116,7 @@ def test_martin_ruiz_diffuse():
     # will fail if default values change
     iam = _iam.martin_ruiz_diffuse(surface_tilt)
     assert_allclose(iam, expected)
+
     # will fail if parameter names change
     iam = _iam.martin_ruiz_diffuse(surface_tilt=surface_tilt, a_r=a_r)
     assert_allclose(iam, expected)
@@ -142,10 +146,13 @@ def test_martin_ruiz_diffuse():
     assert_series_equal(iam[0], expected_sky)
     assert_series_equal(iam[1], expected_gnd)
 
-    # check exception clause
+
+def test_martin_ruiz_diffuse_exception():
+
     with pytest.raises(ValueError):
         _iam.martin_ruiz_diffuse(0.0, a_r=0.0)
 
+
 @requires_scipy
 def test_iam_interp():
 

From 496c7db979c92a893d6f711178be33cc8501f9d1 Mon Sep 17 00:00:00 2001
From: Anton Driesse <anton.driesse@pvperformancelabs.com>
Date: Tue, 22 Oct 2019 14:08:08 +0200
Subject: [PATCH 09/11] Fun with stickler.

---
 pvlib/iam.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pvlib/iam.py b/pvlib/iam.py
index e4e084858d..9d694ab6e5 100644
--- a/pvlib/iam.py
+++ b/pvlib/iam.py
@@ -375,7 +375,7 @@ def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None):
         sin_beta = np.where(surface_tilt < 90, sin(beta), sin(pi - beta))
 
         trig_term_sky = sin_beta + (pi - beta - sin_beta) / (1 + cos(beta))
-        trig_term_gnd = sin_beta +      (beta - sin_beta) / (1 - cos(beta)) # noqa: E222
+        trig_term_gnd = sin_beta +      (beta - sin_beta) / (1 - cos(beta)) # noqa: E222 E261 E501
 
         iam_sky = 1 - exp(-(c1 + c2 * trig_term_sky) * trig_term_sky / a_r)
         iam_gnd = 1 - exp(-(c1 + c2 * trig_term_gnd) * trig_term_gnd / a_r)

From df428901355938f02f1d094b638a6497b638c7ec Mon Sep 17 00:00:00 2001
From: Anton Driesse <anton.driesse@pvperformancelabs.com>
Date: Tue, 29 Oct 2019 13:15:23 +0100
Subject: [PATCH 10/11] Docstring corrections.

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

diff --git a/pvlib/iam.py b/pvlib/iam.py
index 9d694ab6e5..a261724b95 100644
--- a/pvlib/iam.py
+++ b/pvlib/iam.py
@@ -292,13 +292,13 @@ def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None):
         generally on the order of 0.08 to 0.25 for flat-plate PV modules.
         a_r must be greater than zero.
 
-    c1 : numeric scalar
-        First fitting parameter for the expressions that approximates the
+    c1 : float
+        First fitting parameter for the expressions that approximate the
         integral of diffuse irradiance coming from different directions.
         c1 is given as the constant 4 / 3 / pi (0.4244) in [1].
 
-    c2 : numeric scalar
-        Second fitting parameter for the expressions that approximates the
+    c2 : float
+        Second fitting parameter for the expressions that approximate the
         integral of diffuse irradiance coming from different directions.
         If c2 is None, it will be calculated according to the linear
         relationship given in [3].

From 83ffa694c0e2f9e22d131d58f9f41c07bf26f82e Mon Sep 17 00:00:00 2001
From: Anton Driesse <anton.driesse@pvperformancelabs.com>
Date: Fri, 1 Nov 2019 11:08:49 +0100
Subject: [PATCH 11/11] Remove all controversial code.

---
 pvlib/iam.py           | 32 ++++++++++++--------------------
 pvlib/test/test_iam.py | 12 +++---------
 2 files changed, 15 insertions(+), 29 deletions(-)

diff --git a/pvlib/iam.py b/pvlib/iam.py
index a261724b95..f2f968fac3 100644
--- a/pvlib/iam.py
+++ b/pvlib/iam.py
@@ -285,6 +285,7 @@ def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None):
         Surface tilt angles in decimal degrees.
         The tilt angle is defined as degrees from horizontal
         (e.g. surface facing up = 0, surface facing horizon = 90)
+        surface_tilt must be in the range [0, 180]
 
     a_r : numeric
         The angular losses coefficient described in equation 3 of [1].
@@ -348,37 +349,28 @@ def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None):
 
     surface_tilt = np.asanyarray(surface_tilt)
 
-    # undo possible surface rotations
-    with np.errstate(invalid='ignore'):
-        surface_tilt = (surface_tilt + 180) % 360 - 180
-
-    # flip backward tilts forward
-    surface_tilt = np.abs(surface_tilt)
-
     # avoid undefined results for horizontal or upside-down surfaces
-    small_angle = 1e-06
-    surface_tilt = np.clip(surface_tilt, small_angle, 180 - small_angle)
+    zeroang = 1e-06
 
-    if np.any(np.less_equal(a_r, 0)):
-        raise ValueError("The parameter 'a_r' cannot be zero or negative.")
+    surface_tilt = np.where(surface_tilt == 0, zeroang, surface_tilt)
+    surface_tilt = np.where(surface_tilt == 180, 180 - zeroang, surface_tilt)
 
     if c2 is None:
         # This equation is from [3] Sect. 7.2
         c2 = 0.5 * a_r - 0.154
 
-    from numpy import pi, radians, sin, cos, exp
+    beta = np.radians(surface_tilt)
 
-    beta = radians(surface_tilt)
+    from numpy import pi, sin, cos, exp
 
-    with np.errstate(invalid='ignore'):
-        # because sin(pi) isn't exactly zero
-        sin_beta = np.where(surface_tilt < 90, sin(beta), sin(pi - beta))
+    # because sin(pi) isn't exactly zero
+    sin_beta = np.where(surface_tilt < 90, sin(beta), sin(pi - beta))
 
-        trig_term_sky = sin_beta + (pi - beta - sin_beta) / (1 + cos(beta))
-        trig_term_gnd = sin_beta +      (beta - sin_beta) / (1 - cos(beta)) # noqa: E222 E261 E501
+    trig_term_sky = sin_beta + (pi - beta - sin_beta) / (1 + cos(beta))
+    trig_term_gnd = sin_beta +      (beta - sin_beta) / (1 - cos(beta)) # noqa: E222 E261 E501
 
-        iam_sky = 1 - exp(-(c1 + c2 * trig_term_sky) * trig_term_sky / a_r)
-        iam_gnd = 1 - exp(-(c1 + c2 * trig_term_gnd) * trig_term_gnd / a_r)
+    iam_sky = 1 - exp(-(c1 + c2 * trig_term_sky) * trig_term_sky / a_r)
+    iam_gnd = 1 - exp(-(c1 + c2 * trig_term_gnd) * trig_term_gnd / a_r)
 
     if out_index is not None:
         iam_sky = pd.Series(iam_sky, index=out_index, name='iam_sky')
diff --git a/pvlib/test/test_iam.py b/pvlib/test/test_iam.py
index 29152c4b7f..8e04ee2975 100644
--- a/pvlib/test/test_iam.py
+++ b/pvlib/test/test_iam.py
@@ -122,10 +122,10 @@ def test_martin_ruiz_diffuse():
     assert_allclose(iam, expected)
 
     a_r = 0.18
-    surface_tilt = [0, 30, 90, 120, -30, np.nan, np.inf]
-    expected_sky = [0.9407678, 0.9452250, 0.9407678, 0.9055541, 0.9452250,
+    surface_tilt = [0, 30, 90, 120, 180, np.nan, np.inf]
+    expected_sky = [0.9407678, 0.9452250, 0.9407678, 0.9055541, 0.0000000,
                     np.nan, np.nan]
-    expected_gnd = [0.0000000, 0.7610849, 0.9407678, 0.9483508, 0.7610849,
+    expected_gnd = [0.0000000, 0.7610849, 0.9407678, 0.9483508, 0.9407678,
                     np.nan, np.nan]
 
     # check various inputs as list
@@ -147,12 +147,6 @@ def test_martin_ruiz_diffuse():
     assert_series_equal(iam[1], expected_gnd)
 
 
-def test_martin_ruiz_diffuse_exception():
-
-    with pytest.raises(ValueError):
-        _iam.martin_ruiz_diffuse(0.0, a_r=0.0)
-
-
 @requires_scipy
 def test_iam_interp():