Repository: pvlib/pvlib-python
Pull request #752 (2019-07-11 14:53:08)

Implement IEC 61853 IAM calculations
 - [x] Closes  #751 
 - [x] I am familiar with the [contributing guidelines](http://pvlib-python.readthedocs.io/en/latest/contributing.html).
 - [x] Fully tested. Added and/or modified tests to ensure correct behavior for all reasonable inputs. Tests (usually) must pass on the TravisCI and Appveyor testing services.
 - [x] Updates entries to `docs/sphinx/source/api.rst` for API changes.
 - [x] Adds description and name entries in the appropriate `docs/sphinx/source/whatsnew` file for all changes.
 - [x] Code quality and style is sufficient. Passes LGTM and SticklerCI checks.
 - [x] New code is fully documented. Includes sphinx/numpydoc compliant docstrings and comments in the code where necessary.
 - [x] Pull request is nearly complete and ready for detailed review.

Brief description of the problem and proposed solution (if not already fully described in the issue linked to above):
The aforementioned models are not yet in PVLIB and would be useful to the community.

adriesse: Ok, I think there's enough code already to generate some discussion.  Perhaps @jsstein will be interested in following it as well?
cwhanse: Suggest that we separate implementing the Faiman model from the IAM model work, #682 anticipates creating `celltemp.py` to contain related functions.
adriesse: Sure, what's the timeline on `celltemp.py`?
cwhanse: > Sure, what's the timeline on celltemp.py?

This month, next on my list. #682 help is welcomed.
adriesse: What kind of help were you thinking of?
cwhanse: > What kind of help were you thinking of?

Don't know yet. We need to refresh the PR in #682 with master and see where it's at, been a while since I started it.  
adriesse: I guess I can just work on the code where it is and move it after your PR is merged, but before my PR is merged.
cwhanse: Personally I'm not a fan of the greek symbol, I don't know how to type it and the IDE I frequently use (Spyder) doesn't have a Insert Symbol function.

Feel free to edit the docstring from MATLAB, for style and clarity.  If you'd like I'll suggest some changes in comments.
adriesse: Not too many comments coming... do I need to invite people?  @wholmgren , @markcampanelli , @mikofski , @jsstein.  The work is not done, but I would like to receive early nudges to keep me on the straight and narrow.
mikofski: Sorry, yes gentle nudges sometimes useful. It's not for lack of interest, just busy, too many competing priorities. I'll try to take a quick look, for code and style, but Cliff or Daniel Riley probably best for content since this is ported from MATLAB right?
adriesse: Thanks Mark--I ported only comments from MATLAB, not code.
adriesse: Refresh.
cwhanse: > Refresh.

@adriesse What action are you expecting from the maintainers and reviewers? Reviewing the comments above, I think we are saying 1) the docstring could use some editing, 2) let's not use the Unicode symbol as the variable, 3) perhaps not use the annotation. 
adriesse: The discussion didn't seem like it was finished, and I thought other people might want to say something.  Anyway, I think it's done now.
wholmgren: I prefer descriptive arguments like “aoi” instead of “theta” but up for discussion in this case. 
adriesse: @wholmgren It says that only those with write access to this repository can merge pull requests.
cwhanse: > @adriesse are you able to resolve the merge conflict in the whatsnew file using your git client or the web interface? See here for web interface instructions, though I'm not sure if that's a maintainer-only feature. I can do it if you'd prefer.

I took care of the conflict in the whatsnew file.

Patch:
diff --git a/docs/sphinx/source/whatsnew/v0.7.0.rst b/docs/sphinx/source/whatsnew/v0.7.0.rst
index a5251d9..505bd00 100644
--- a/docs/sphinx/source/whatsnew/v0.7.0.rst
+++ b/docs/sphinx/source/whatsnew/v0.7.0.rst
@@ -10,10 +10,16 @@ compatibility notes.
 **Python 2.7 support ended on June 1, 2019.** (:issue:`501`)
 **Minimum numpy version is now 1.10.4. Minimum pandas version is now 0.18.1.**
 
+Enhancements
+~~~~~~~~~~~~
+* Created two new incidence angle modifier functions: :py:func:`pvlib.pvsystem.iam_martin_ruiz`
+  and :py:func:`pvlib.pvsystem.iam_interp`. (:issue:`751`)
+
 Bug fixes
 ~~~~~~~~~
 * Fix handling of keyword arguments in `forecasts.get_processed_data`.
   (:issue:`745`)
+* Fix output as Series feature in :py:func:`pvlib.pvsystem.ashraeiam`.
 
 Testing
 ~~~~~~~
@@ -26,4 +32,5 @@ Contributors
 * Mark Campanellli (:ghuser:`markcampanelli`)
 * Will Holmgren (:ghuser:`wholmgren`)
 * Oscar Dowson (:ghuser:`odow`)
+* Anton Driesse (:ghuser:`adriesse`)
 * Alexander Morgan (:ghuser:`alexandermorgan`)
diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 35701bd..9000b9f 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -17,7 +17,7 @@ import numpy as np
 import pandas as pd
 
 from pvlib import atmosphere, irradiance, tools, singlediode as _singlediode
-from pvlib.tools import _build_kwargs
+from pvlib.tools import _build_kwargs, cosd
 from pvlib.location import Location
 
 
@@ -946,7 +946,7 @@ def ashraeiam(aoi, b=0.05):
     iam = np.where(aoi_gte_90, 0, iam)
     iam = np.maximum(0, iam)
 
-    if isinstance(iam, pd.Series):
+    if isinstance(aoi, pd.Series):
         iam = pd.Series(iam, index=aoi.index)
 
     return iam
@@ -1064,6 +1064,161 @@ def physicaliam(aoi, n=1.526, K=4., L=0.002):
     return iam
 
 
+def iam_martin_ruiz(aoi, a_r=0.16):
+    '''
+    Determine the incidence angle modifier (iam) using the Martin
+    and Ruiz incident angle model.
+
+    Parameters
+    ----------
+    aoi : numeric, degrees
+        The angle of incidence between the module normal vector and the
+        sun-beam vector in degrees. Theta must be a numeric scalar or vector.
+        iam is 0 where |aoi| > 90.
+
+    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 or vector (same length as aoi).
+
+    Returns
+    -------
+    iam : numeric
+        The incident angle modifier(s)
+
+    Notes
+    -----
+    iam_martin_ruiz calculates the incidence angle modifier (iam)
+    as described by Martin and Ruiz in [1]. The information
+    required is the incident angle (aoi) and the angular losses
+    coefficient (a_r). Please note that [1] has a corrigendum which makes
+    the document much simpler to understand.
+
+    The incident angle modifier is defined as
+    [1-exp(-cos(aoi/ar))] / [1-exp(-1/ar)], which is
+    presented as AL(alpha) = 1 - IAM in equation 4 of [1]. Thus IAM is
+    equal to 1 at aoi = 0, and equal to 0 at aoi = 90.  This equation is only
+    valid for -90 <= aoi <= 90, therefore iam must be constrained to 0.0
+    beyond this range.
+
+    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
+    --------
+    physicaliam
+    ashraeiam
+    iam_interp
+    '''
+    # Contributed by Anton Driesse (@adriesse), PV Performance Labs. July, 2019
+
+    aoi_input = aoi
+
+    aoi = np.asanyarray(aoi)
+    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.")
+
+    with np.errstate(invalid='ignore'):
+        iam = (1 - np.exp(-cosd(aoi) / a_r)) / (1 - np.exp(-1 / a_r))
+        iam = np.where(np.abs(aoi) >= 90.0, 0.0, iam)
+
+    if isinstance(aoi_input, pd.Series):
+        iam = pd.Series(iam, index=aoi_input.index)
+
+    return iam
+
+
+def iam_interp(aoi, theta_ref, iam_ref, method='linear', normalize=True):
+    '''
+    Determine the incidence angle modifier (iam) by interpolating a set of
+    reference values, which are usually measured values.
+
+    Parameters
+    ----------
+    aoi : numeric, degrees
+        The angle of incidence between the module normal vector and the
+        sun-beam vector in degrees.
+
+    theta_ref : numeric, degrees
+        Vector of angles at which the iam is known.
+
+    iam_ref : numeric, unitless
+        iam values for each angle in theta_ref.
+
+    method : str, default 'linear'
+        Specifies the interpolation method.
+        Useful options are: 'linear', 'quadratic','cubic'.
+        See scipy.interpolate.interp1d for more options.
+
+    normalize : boolean
+        When true, the interpolated values are divided by the interpolated
+        value at zero degrees.  This ensures that the iam at normal
+        incidence is equal to 1.0.
+
+    Returns
+    -------
+    iam : numeric
+        The incident angle modifier(s)
+
+    Notes:
+    ------
+    theta_ref must have two or more points and may span any range of angles.
+    Typically there will be a dozen or more points in the range 0-90 degrees.
+    iam beyond the range of theta_ref are extrapolated, but constrained to be
+    non-negative.
+
+    The sign of aoi is ignored; only the magnitude is used.
+
+    See Also
+    --------
+    physicaliam
+    ashraeiam
+    iam_martin_ruiz
+    '''
+    # Contributed by Anton Driesse (@adriesse), PV Performance Labs. July, 2019
+
+    from scipy.interpolate import interp1d
+
+    # Scipy doesn't give the clearest feedback, so check number of points here.
+    MIN_REF_VALS = {'linear': 2, 'quadratic': 3, 'cubic': 4, 1: 2, 2: 3, 3: 4}
+
+    if len(theta_ref) < MIN_REF_VALS.get(method, 2):
+        raise ValueError("Too few reference points defined "
+                         "for interpolation method '%s'." % method)
+
+    if np.any(np.less(iam_ref, 0)):
+        raise ValueError("Negative value(s) found in 'iam_ref'. "
+                         "This is not physically possible.")
+
+    interpolator = interp1d(theta_ref, iam_ref, kind=method,
+                            fill_value='extrapolate')
+    aoi_input = aoi
+
+    aoi = np.asanyarray(aoi)
+    aoi = np.abs(aoi)
+    iam = interpolator(aoi)
+    iam = np.clip(iam, 0, None)
+
+    if normalize:
+        iam /= interpolator(0)
+
+    if isinstance(aoi_input, pd.Series):
+        iam = pd.Series(iam, index=aoi_input.index)
+
+    return iam
+
+
 def calcparams_desoto(effective_irradiance, temp_cell,
                       alpha_sc, a_ref, I_L_ref, I_o_ref, R_sh_ref, R_s,
                       EgRef=1.121, dEgdT=-0.0002677,


Test patch:
diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index 17f70e9..7de8930 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -148,6 +148,88 @@ def test_PVSystem_physicaliam(mocker):
     assert iam < 1.
 
 
+def test_iam_martin_ruiz():
+
+    aoi = 45.
+    a_r = 0.16
+    expected = 0.98986965
+
+    # will fail of default values change
+    iam = pvsystem.iam_martin_ruiz(aoi)
+    assert_allclose(iam, expected)
+    # will fail of parameter names change
+    iam = pvsystem.iam_martin_ruiz(aoi=aoi, a_r=a_r)
+    assert_allclose(iam, expected)
+
+    a_r = 0.18
+    aoi = [-100, -60, 0, 60, 100, np.nan, np.inf]
+    expected = [0.0, 0.9414631, 1.0, 0.9414631, 0.0, np.nan, 0.0]
+
+    # check out of range of inputs as list
+    iam = pvsystem.iam_martin_ruiz(aoi, a_r)
+    assert_allclose(iam, expected, equal_nan=True)
+
+    # check out of range of inputs as array
+    iam = pvsystem.iam_martin_ruiz(np.array(aoi), a_r)
+    assert_allclose(iam, expected, equal_nan=True)
+
+    # check out of range of inputs as Series
+    aoi = pd.Series(aoi)
+    expected = pd.Series(expected)
+    iam = pvsystem.iam_martin_ruiz(aoi, a_r)
+    assert_series_equal(iam, expected)
+
+    # check exception clause
+    with pytest.raises(RuntimeError):
+        pvsystem.iam_martin_ruiz(0.0, a_r=0.0)
+
+
+@requires_scipy
+def test_iam_interp():
+
+    aoi_meas = [0.0, 45.0, 65.0, 75.0]
+    iam_meas = [1.0,  0.9,  0.8,  0.6]
+
+    # simple default linear method
+    aoi = 55.0
+    expected = 0.85
+    iam = pvsystem.iam_interp(aoi, aoi_meas, iam_meas)
+    assert_allclose(iam, expected)
+
+    # simple non-default method
+    aoi = 55.0
+    expected = 0.8878062
+    iam = pvsystem.iam_interp(aoi, aoi_meas, iam_meas, method='cubic')
+    assert_allclose(iam, expected)
+
+    # check with all reference values
+    aoi = aoi_meas
+    expected = iam_meas
+    iam = pvsystem.iam_interp(aoi, aoi_meas, iam_meas)
+    assert_allclose(iam, expected)
+
+    # check normalization and Series
+    aoi = pd.Series(aoi)
+    expected = pd.Series(expected)
+    iam_mult = np.multiply(0.9, iam_meas)
+    iam = pvsystem.iam_interp(aoi, aoi_meas, iam_mult, normalize=True)
+    assert_series_equal(iam, expected)
+
+    # check beyond reference values
+    aoi = [-45, 0, 45, 85, 90, 95, 100, 105, 110]
+    expected = [0.9, 1.0, 0.9, 0.4, 0.3, 0.2, 0.1, 0.0, 0.0]
+    iam = pvsystem.iam_interp(aoi, aoi_meas, iam_meas)
+    assert_allclose(iam, expected)
+
+    # check exception clause
+    with pytest.raises(ValueError):
+        pvsystem.iam_interp(0.0, [0], [1])
+
+    # check exception clause
+    with pytest.raises(ValueError):
+        pvsystem.iam_interp(0.0, [0, 90], [1, -1])
+
+
 # if this completes successfully we'll be able to do more tests below.
 @pytest.fixture(scope="session")
 def sam_data():