Implement iam.schlick, iam.schlick_diffuse

docs/sphinx/source/reference/pv_modeling.rst

@@ -28,6 +28,8 @@ Incident angle modifiers
    iam.interp
    iam.marion_diffuse
    iam.marion_integrate
+   iam.fedis
+   iam.schlick
 
 PV temperature models
 ---------------------

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

@@ -10,6 +10,10 @@ Deprecations
 Enhancements
 ~~~~~~~~~~~~
 * Multiple code style issues fixed that were reported by LGTM analysis. (:issue:`1275`, :pull:`1559`)
+* Added a simple IAM model :py:func:`pvlib.iam.schlick` and made it available
+  in :py:func:`~pvlib.iam.marion_diffuse` via ``model='schlick'`` (:pull:`1562`, :issue:`1564`)
+* Added a direct and diffuse IAM model :py:func:`pvlib.iam.fedis` (:pull:`1562`)
+
 
 Bug fixes
 ~~~~~~~~~
@@ -34,4 +38,5 @@ Requirements
 
 Contributors
 ~~~~~~~~~~~~
-* Christian Orner (:ghuser:`chrisorner`)
+* Christian Orner (:ghuser:`chrisorner`)
+* Yu Xie (:ghuser:`xieyupku`)

pvlib/iam.py

@@ -541,7 +541,7 @@ def marion_diffuse(model, surface_tilt, **kwargs):
     ----------
     model : str
         The IAM function to evaluate across solid angle. Must be one of
-        `'ashrae', 'physical', 'martin_ruiz', 'sapm'`.
+        `'ashrae', 'physical', 'martin_ruiz', 'sapm', 'schlick'`.
 
     surface_tilt : numeric
         Surface tilt angles in decimal degrees.
@@ -592,6 +592,7 @@ def marion_diffuse(model, surface_tilt, **kwargs):
         'ashrae': ashrae,
         'sapm': sapm,
         'martin_ruiz': martin_ruiz,
+        'schlick': schlick,
     }
 
     try:
@@ -748,3 +749,177 @@ def marion_integrate(function, surface_tilt, region, num=None):
         Fd = pd.Series(Fd, surface_tilt.index)
 
     return Fd
+
+
+def schlick(aoi):
+    """
+    Determine incidence angle modifier (IAM) using the Schlick approximation
+    to the Fresnel equations.
+
+    The Schlick approximation was proposed in [1]_ as a computationally
+    efficient alternative to computing the Fresnel factor in computer
+    graphics contexts.  This implementation is a normalized form of the
+    equation in [1]_ so that it can be used as a PV IAM model.
+    Unlike other IAM models, this model has no ability to describe
+    different reflection profiles.
+
+    In PV contexts, the Schlick approximation has been used as an analytically
+    integrable alternative to the Fresnel equations for estimating IAM
+    for diffuse irradiance [2]_.
+
+    Parameters
+    ----------
+    aoi : numeric
+        The angle of incidence (AOI) between the module normal vector and the
+        sun-beam vector in degrees. Angles of nan will result in nan.
+
+    Returns
+    -------
+    iam : numeric
+        The incident angle modifier
+
+    References
+    ----------
+    .. [1] Schlick, C. An inexpensive BRDF model for physically-based
+       rendering. Computer graphics forum 13 (1994).
+
+    .. [2] Xie, Y., M. Sengupta, A. Habte, A. Andreas, "The 'Fresnel Equations'
+       for Diffuse radiation on Inclined photovoltaic Surfaces (FEDIS)",
+       Renewable and Sustainable Energy Reviews, vol. 161, 112362. June 2022.
+       :doi:`10.1016/j.rser.2022.112362`
+
+    See Also
+    --------
+    pvlib.iam.fedis
+    """
+    iam = 1 - (1 - cosd(aoi)) ** 5
+    iam = np.where(np.abs(aoi) >= 90.0, 0.0, iam)
+
+    # preserve input type
+    if np.isscalar(aoi):
+        iam = iam.item()
+    elif isinstance(aoi, pd.Series):
+        iam = pd.Series(iam, aoi.index)
+
+    return iam
+
+
+def fedis(aoi, surface_tilt, n=1.5, n_ref=None):
+    """
+    Determine the incidence angle modifiers (IAM) for direct, diffuse sky,
+    and ground-reflected radiation using the FEDIS transmittance model.
+
+    The "Fresnel Equations" for Diffuse radiation on Inclined photovoltaic
+    Surfaces (FEDIS) [1]_ is the result of analytical integration
+    the Schlick approximation [2]_ to the Fresnel equations.
+
+    Parameters
+    ----------
+    aoi : numeric
+        Angle of incidence. [degrees]
+
+    surface_tilt : numeric
+        Surface tilt angle measured from horizontal (e.g. surface facing
+        up = 0, surface facing horizon = 90). [degrees]
+
+    n : float, default 1.5
+        Refractive index of the PV front surface material.  The default value
+        of 1.5 was used for an IMT reference cell in [1]_. [unitless]
+
+    n_ref : float, optional
+        Reference refractive index. In [1]_ this was set to 1.4585 for
+        was used for a fused silica dome over a CMP22, but in conventional
+        PV applications it is appropriate to set this to the same value as
+        ``n`` (the default behavior).
+
+    Returns
+    -------
+    iam : dict
+        IAM values for each type of irradiance:
+
+            * 'direct': irradiance from the solar disc
+            * 'sky': irradiance radiation from the sky dome
+            * 'ground': irradiance reflected from the ground
+
+    Notes
+    -----
+    This implementation corrects a typo in [1]_ regarding the sign
+    of the last polynomial term in Equation 5.
+
+    References
+    ----------
+    .. [1] Xie, Y., M. Sengupta, A. Habte, A. Andreas, "The 'Fresnel Equations'
+       for Diffuse radiation on Inclined photovoltaic Surfaces (FEDIS)",
+       Renewable and Sustainable Energy Reviews, vol. 161, 112362. June 2022.
+       :doi:`10.1016/j.rser.2022.112362`
+
+    .. [2] Schlick, C. An inexpensive BRDF model for physically-based
+       rendering. Computer graphics forum 13 (1994).
+
+    See Also
+    --------
+    pvlib.iam.schlick
+    """
+    if n_ref is None:
+        n_ref = n
+
+    # angle between module normal and refracted ray:
+    theta_0tp = asind(sind(aoi) / n)  # Eq 3c
+
+    # reflectance of direct radiation on PV cover:
+    with np.errstate(invalid='ignore'):
+        sin_term = sind(aoi - theta_0tp)**2 / sind(aoi + theta_0tp)**2 / 2
+        tan_term = tand(aoi - theta_0tp)**2 / tand(aoi + theta_0tp)**2 / 2
+
+    rd = sin_term + tan_term  # Eq 3b
+    rd = np.where(np.abs(aoi) < 1e-6, ((n-1.0)/(n+1.0))**2.0, rd)
+
+    # reflectance for normal incidence with reference refractive index:
+    r0 = ((n_ref-1.0)/(n_ref+1.0))**2.0  # Eq 3e
+
+    # relative transmittance of direct radiation by PV cover:
+    cd = (1 - rd) / (1 - r0)  # Eq 3a
+    cd = np.where(np.abs(aoi) > 90, 0, cd)
+
+    # weighting function
+    term1 = n*(n_ref+1)**2 / (n_ref*(n+1)**2)
+    # note: the last coefficient here differs in sign from the reference
+    polycoeffs = [2.77526e-09, 3.74953, -5.18727, 3.41186, -1.08794, 0.136060]
+    term2 = np.polynomial.polynomial.polyval(n, polycoeffs)
+    w = term1 * term2  # Eq 5
+
+    # relative transmittance of sky diffuse radiation by PV cover:
+    cosB = cosd(surface_tilt)
+    sinB = sind(surface_tilt)
+    cuk = (2*w / (np.pi * (1 + cosB))) * (
+        (30/7)*np.pi - (160/21)*np.radians(surface_tilt) - (10/3)*np.pi*cosB
+        + (160/21)*cosB*sinB - (5/3)*np.pi*cosB*sinB**2 + (20/7)*cosB*sinB**3
+        - (5/16)*np.pi*cosB*sinB**4 + (16/105)*cosB*sinB**5
+    )  # Eq 4
+
+    # relative transmittance of ground-reflected radiation by PV cover:
+    with np.errstate(divide='ignore', invalid='ignore'):  # Eq 6
+        cug = 40 * w / (21 * (1 - cosB)) - (1 + cosB) / (1 - cosB) * cuk
+
+    cug = np.where(surface_tilt < 1e-6, 0, cug)
+
+    # respect input types:
+    if np.isscalar(aoi):
+        cd = cd.item()
+    elif isinstance(aoi, pd.Series):
+        cd = pd.Series(cd, aoi.index)
+
+    if np.isscalar(surface_tilt):
+        cuk = cuk.item()
+        cug = cug.item()
+    elif isinstance(surface_tilt, pd.Series):
+        cuk = pd.Series(cuk, surface_tilt.index)
+        cug = pd.Series(cug, surface_tilt.index)
+
+    out = {
+        'direct': cd,
+        'sky': cuk,
+        'ground': cug,
+    }
+
+    return out

pvlib/tests/test_iam.py

@@ -322,3 +322,89 @@ def test_marion_integrate_invalid():
 
     with pytest.raises(ValueError):
         _iam.marion_integrate(_iam.ashrae, 0, 'bad', 180)
+
+
+def test_schlick():
+    idx = pd.date_range('2019-01-01', freq='h', periods=9)
+    aoi = pd.Series([-180, -135, -90, -45, 0, 45, 90, 135, 180], idx)
+    expected = pd.Series([0, 0, 0, 0.99784451, 1.0, 0.99784451, 0, 0, 0], idx)
+
+    # scalars
+    for aoi_scalar, expected_scalar in zip(aoi, expected):
+        actual = _iam.schlick(aoi_scalar)
+        assert_allclose(expected_scalar, actual)
+
+    # numpy arrays
+    actual = _iam.schlick(aoi.values)
+    assert_allclose(expected.values, actual)
+
+    # pandas Series
+    actual = _iam.schlick(aoi)
+    assert_series_equal(expected, actual)
+
+
+def test_fedis_defaults_direct():
+    idx = pd.date_range('2019-01-01', freq='h', periods=9)
+    aoi = pd.Series([-180, -135, -90, -45, 0, 45, 90, 135, 180], idx)
+    expected = pd.Series([0, 0, 0, 0.989333426, 1, 0.989333426, 0, 0, 0], idx)
+
+    # numpy arrays
+    actual = _iam.fedis(aoi.values, surface_tilt=0)
+    assert_allclose(expected, expected, atol=1e-6)
+
+    # scalars
+    for i in range(len(aoi)):
+        actual = _iam.fedis(aoi[i], surface_tilt=0)
+        assert_allclose(expected[i], actual['direct'], atol=1e-6)
+
+    # pandas Series
+    actual = _iam.fedis(aoi, surface_tilt=0)
+    assert_series_equal(expected, actual['direct'])
+
+
+def test_fedis_defaults_diffuse():
+    idx = pd.date_range('2019-01-01', freq='h', periods=3)
+    surface_tilt = pd.Series([0, 20, 90], idx)
+
+    # expected values generated with code from model authors:
+    # https://github.com/NREL/FEDIS/commit/7ae7186caa39aa85848163a39dac46df56fb9819  # noqa: E501
+    expected = {
+        'sky': pd.Series([0.946166074, 0.956218435, 0.946166074], idx),
+        'ground': pd.Series([0.0, 0.62284759, 0.946166074], idx),
+    }
+
+    # numpy arrays
+    actual = _iam.fedis(aoi=0, surface_tilt=surface_tilt.values)
+    assert_allclose(expected['sky'], actual['sky'])
+    assert_allclose(expected['ground'], actual['ground'])
+
+    # scalars
+    for i in range(len(surface_tilt)):
+        actual = _iam.fedis(aoi=0, surface_tilt=surface_tilt[i])
+        assert_allclose(expected['sky'][i], actual['sky'])
+        assert_allclose(expected['ground'][i], actual['ground'])
+
+    # pandas Series
+    actual = _iam.fedis(aoi=0, surface_tilt=surface_tilt)
+    assert_series_equal(expected['sky'], actual['sky'])
+    assert_series_equal(expected['ground'], actual['ground'])
+
+
+def test_fedis_kwargs():
+    # custom n and n_ref
+
+    surface_tilt = np.array([0, 30, 60, 90])
+    aoi = np.array([0, 30, 60, 90])
+
+    # expected values generated with code from model authors:
+    # https://github.com/NREL/FEDIS/commit/7ae7186caa39aa85848163a39dac46df56fb9819  # noqa: E501
+    expected = {
+        'direct': np.array([0.948928986, 0.947089588, 0.894889901, 0.0]),
+        'sky': np.array([0.89536659, 0.90815772, 0.90728496, 0.89536659]),
+        'ground': np.array([0.0,  0.717209232, 0.859611482, 0.895366593]),
+    }
+    # numpy arrays
+    actual = _iam.fedis(aoi, surface_tilt, n=1.7, n_ref=1.3)
+    assert_allclose(expected['direct'], actual['direct'], atol=1e-6)
+    assert_allclose(expected['sky'], actual['sky'])
+    assert_allclose(expected['ground'], actual['ground'], atol=1e-6)