Implement IEC 61853 IAM calculations for diffuse irradiance (#793)

* Functional new function.

* Undo

* Working function; partial docstring; no tests.

* Complete docstring; simplify a bit; start tests.

* Getting better all the time!

* Adjust some comments and update whatsnew.

* Improvements as per review.

* More changes as requested.

* Fun with stickler.

* Docstring corrections.

* Remove all controversial code.

Author: adriesse

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
 

docs/sphinx/source/whatsnew/v0.7.0.rst

@@ -114,6 +114,8 @@ Enhancements
 * Add :py:func:`~pvlib.ivtools.fit_sdm_desoto`, a method to fit the De Soto single
   diode model to the typical specifications given in manufacturers datasheets.
 * 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
 ~~~~~~~~~

pvlib/iam.py

@@ -12,7 +12,6 @@
 import pandas as pd
 from pvlib.tools import cosd, sind, tand, asind
 
-
 # a dict of required parameter names for each IAM model
 # keys are the function names for the IAM models
 IAM_MODEL_PARAMS = {
@@ -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
@@ -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,6 +274,111 @@ def martin_ruiz(aoi, a_r=0.16):
     return iam
 
 
+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
+    ----------
+    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)
+        surface_tilt must be in the range [0, 180]
+
+    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 greater than zero.
+
+    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 : 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].
+
+    Returns
+    -------
+    iam_sky : numeric
+        The incident angle modifier for sky diffuse
+
+    iam_ground : numeric
+        The incident angle modifier for ground-reflected diffuse
+
+    Notes
+    -----
+    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
+    ----------
+    [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.
+
+    [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
+    iam.sapm
+    '''
+    # Contributed by Anton Driesse (@adriesse), PV Performance Labs. Oct. 2019
+
+    if isinstance(surface_tilt, pd.Series):
+        out_index = surface_tilt.index
+    else:
+        out_index = None
+
+    surface_tilt = np.asanyarray(surface_tilt)
+
+    # avoid undefined results for horizontal or upside-down surfaces
+    zeroang = 1e-06
+
+    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
+
+    beta = np.radians(surface_tilt)
+
+    from numpy import pi, sin, cos, exp
+
+    # 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
+
+    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

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,11 +100,53 @@ def test_martin_ruiz():
     iam = _iam.martin_ruiz(aoi, a_r)
     assert_series_equal(iam, expected)
 
-    # check exception clause
-    with pytest.raises(RuntimeError):
+
+def test_martin_ruiz_exception():
+
+    with pytest.raises(ValueError):
         _iam.martin_ruiz(0.0, a_r=0.0)
 
 
+def test_martin_ruiz_diffuse():
+
+    surface_tilt = 30.
+    a_r = 0.16
+    expected = (0.9549735, 0.7944426)
+
+    # 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)
+
+    a_r = 0.18
+    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.9407678,
+                    np.nan, np.nan]
+
+    # check various inputs as list
+    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(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
+    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(surface_tilt, a_r)
+    assert_series_equal(iam[0], expected_sky)
+    assert_series_equal(iam[1], expected_gnd)
+
+
 @requires_scipy
 def test_iam_interp():