Add marion_diffuse_tracking, to be used instead of marion_diffuse for vector surface_tilt inputs

docs/sphinx/source/reference/pv_modeling/iam.rst

@@ -14,6 +14,7 @@ Incident angle modifiers
    iam.sapm
    iam.interp
    iam.marion_diffuse
+   iam.marion_diffuse_tracking
    iam.marion_integrate
    iam.schlick
    iam.schlick_diffuse

pvlib/iam.py

@@ -12,6 +12,7 @@
 import pandas as pd
 import functools
 from scipy.optimize import minimize
+from scipy.interpolate import PchipInterpolator
 from pvlib.tools import cosd, sind, acosd
 
 # a dict of required parameter names for each IAM model
@@ -644,6 +645,105 @@ def marion_diffuse(model, surface_tilt, **kwargs):
     return iam
 
 
+def marion_diffuse_tracking(model, surface_tilt, resolution=0.5, **kwargs):
+    """
+    Determine diffuse irradiance incidence angle modifiers using Marion's
+    method of integrating over solid angle. This function is designed for
+    trackers, where ``surface_tilt`` is a vector, to avoid the computational
+    burden of calling ``marion_integrate`` for each tilt angle.
+    Instead, the IAM function is integrated once for a range of tilt angles,
+    and then interpolated to the requested tilt angles. For fixed-tilt systems,
+    use :py:func:`marion_diffuse`.
+
+    Parameters
+    ----------
+    model : str
+        The IAM function to evaluate across solid angle. Must be one of
+        `'ashrae', 'physical', 'martin_ruiz', 'sapm', 'schlick'`.
+
+    surface_tilt : numeric
+        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).
+
+    resolution : float, default 0.5
+        The resolution in degrees of the tilt angle vector used for
+        the integration.
+
+    **kwargs
+        Extra parameters passed to the IAM function.
+
+    Returns
+    -------
+    iam : dict
+        IAM values for each type of diffuse irradiance:
+
+        * 'sky': radiation from the sky dome (zenith <= 90)
+        * 'horizon': radiation from the region of the sky near the horizon
+          (89.5 <= zenith <= 90)
+        * 'ground': radiation reflected from the ground (zenith >= 90)
+
+        See [1]_ for a detailed description of each class.
+
+    See Also
+    --------
+    pvlib.iam.marion_diffuse
+    pvlib.iam.marion_integrate
+
+    References
+    ----------
+    .. [1] B. Marion "Numerical method for angle-of-incidence correction
+       factors for diffuse radiation incident photovoltaic modules",
+       Solar Energy, Volume 147, Pages 344-348. 2017.
+       :doi:`10.1016/j.solener.2017.03.027`
+    """
+
+    models = {
+        'physical': physical,
+        'ashrae': ashrae,
+        'sapm': sapm,
+        'martin_ruiz': martin_ruiz,
+        'schlick': schlick,
+    }
+
+    try:
+        iam_model = models[model]
+    except KeyError:
+        raise ValueError('model must be one of: ' + str(list(models.keys())))
+
+    full_range = (np.asarray(surface_tilt) > 90).any()
+
+    iam_function = functools.partial(iam_model, **kwargs)
+    iam = {}
+    for region in ['sky', 'horizon', 'ground']:
+        interpolator = _get_marion_interpolator(iam_function,
+                                                region, resolution, full_range)
+        iam_values = interpolator(surface_tilt)
+        if isinstance(surface_tilt, pd.Series):
+            iam_values = pd.Series(iam_values, index=surface_tilt.index)
+        elif isinstance(surface_tilt, list):
+            iam_values = iam_values.tolist()
+        iam[region] = iam_values
+
+    return iam
+
+
+@functools.cache
+def _get_marion_interpolator(iam_function, region, resolution,
+                             full_range=False):
+    """
+    Cached interpolator for the Marion integration of an IAM function over
+    solid angle. Helper function for :py:func:`marion_efficient` function
+    to avoid repeated calculations leading to excessive memory use.
+    """
+    if full_range:
+        tilt = np.arange(0, 180.5, resolution)
+    else:
+        tilt = np.arange(0, 90.5, resolution)
+    iam = marion_integrate(iam_function, tilt, region)
+    return PchipInterpolator(tilt, iam)
+
+
 def marion_integrate(function, surface_tilt, region, num=None):
     """
     Integrate an incidence angle modifier (IAM) function over solid angle

tests/test_iam.py

@@ -289,6 +289,60 @@ def test_marion_diffuse_invalid():
         _iam.marion_diffuse('not_a_model', 20)
 
 
+def test_marion_diffuse_tracking_list():
+    # tilt angles under 90
+    expected = {
+        'sky': [0.9523612, 0.95960858, 0.96198112],
+        'horizon': [0.0, 0.83290704, 0.89872877],
+        'ground': [0.0, 0.71982356, 0.81863602]
+    }
+    tilt = [0, 20, 30]
+    actual = _iam.marion_diffuse_tracking('ashrae', tilt)
+    for key, value in expected.items():
+        assert_allclose(actual[key], value)
+    # tilt angles above 90
+    expected = {
+        'sky': [0.9523612, 0.95960858, 0.94530815],
+        'horizon': [0.0, 0.83290704, 0.97141949],
+        'ground': [0.0, 0.71982356, 0.95735835]
+    }
+    tilt = [0, 20, 100]
+    actual = _iam.marion_diffuse_tracking('ashrae', tilt)
+    for key, values in expected.items():
+        assert_allclose(actual[key], values)
+
+
+def test_marion_diffuse_tracking_array():
+    expected = {
+        'sky': np.array([0.9523612, 0.95960858, 0.96198112]),
+        'horizon': np.array([0.0, 0.83290704, 0.89872877]),
+        'ground': np.array([0.0, 0.71982356, 0.81863602])
+    }
+    tilt = np.array([0, 20, 30])
+    actual = _iam.marion_diffuse_tracking('ashrae', tilt)
+    for key, value in expected.items():
+        assert_allclose(actual[key], value)
+
+
+def test_marion_diffuse_tracking_series():
+    expected = {
+        'sky': [0.9523612, 0.95960858, 0.96198112],
+        'horizon': [0.0, 0.83290704, 0.89872877],
+        'ground': [0.0, 0.71982356, 0.81863602]
+    }
+    idx = pd.date_range('2019-01-01', periods=3, freq='h')
+    tilt = pd.Series([0, 20, 30], index=idx)
+    actual = _iam.marion_diffuse_tracking('ashrae', tilt)
+    for key, values in expected.items():
+        expected_series = pd.Series(values, index=idx)
+        assert_series_equal(actual[key], expected_series)
+
+
+def test_marion_diffuse_tracking_invalid():
+    with pytest.raises(ValueError):
+        _iam.marion_diffuse_tracking('not_a_model', 20)
+
+
 @pytest.mark.parametrize('region,N,expected', [
     ('sky', 180, 0.9596085829811408),
     ('horizon', 1800, 0.8329070417832541),