Implement iam.schlick, iam.schlick_diffuse

docs/sphinx/source/reference/pv_modeling.rst

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

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

@@ -10,6 +10,9 @@ Deprecations
 Enhancements
 ~~~~~~~~~~~~
 * Multiple code style issues fixed that were reported by LGTM analysis. (:issue:`1275`, :pull:`1559`)
+* Added :py:func:`pvlib.iam.fedis`, a model for calculating direct and diffuse IAM
+  coefficients. (:pull:`1562`)
+
 
 Bug fixes
 ~~~~~~~~~
@@ -34,4 +37,5 @@ Requirements
 
 Contributors
 ~~~~~~~~~~~~
-* Christian Orner (:ghuser:`chrisorner`)
+* Christian Orner (:ghuser:`chrisorner`)
+* Yu Xie (:ghuser:`xieyupku`)

pvlib/iam.py

@@ -748,3 +748,108 @@ def marion_integrate(function, surface_tilt, region, num=None):
         Fd = pd.Series(Fd, surface_tilt.index)
 
     return Fd
+
+
+def fedis(aoi, surface_tilt, n=1.5, n_ref=1.4585):
+    """
+    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 an analytical solution of diffuse transmission
+    based on the rigorous integration of an alternate form of the
+    Fresnel equations. The approach leads to a simple yet accurate
+    relative transmittance model that reconciles the solar energy
+    sensed by pyranometers and PV panels.
+
+    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 cover.  The default value of 1.5
+        was used for an IMT reference cell in [1]_. [unitless]
+
+    n_ref : float, default 1.4585
+        Refractive index of the pyranometer cover. The default value
+        was used for a fused silica dome over a CMP22 in [1]_.
+
+    Returns
+    -------
+    iam : dict
+        IAM values for each type of irradiance:
+
+            * 'direct': radiation from the solar disc
+            * 'sky': radiation from the sky dome (zenith <= 90)
+            * 'ground': radiation reflected from the ground (zenith >= 90)
+
+    Notes
+    -----
+    This implementation corrects a typo in the reference 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`
+    """
+
+    # 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 >= 90, 90 - zeroang, surface_tilt)
+
+    # and for aoi:
+    aoi = np.where(aoi <= 0, zeroang, aoi)
+    # similar for AOI > 90
+    aoi = np.where(aoi >= 90, 90 - zeroang, aoi)
+
+    # angle between module normal and refracted ray:
+    theta_0tp = asind(sind(aoi) / n)  # Eq 3c
+
+    # reflectance of direct radiation on PV cover:
+    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
+
+    # reflectance on pyranometer cover:
+    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
+
+    # 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:
+    cug = 40 * w / (21 * (1 - cosB)) - (1 + cosB) / (1 - cosB) * cuk  # Eq 6
+
+    # handle tilt=0 case correctly:
+    cug = np.where(surface_tilt == zeroang, 0, cug)
+    out = {
+        'direct': cd,
+        'sky': cuk,
+        'ground': cug,
+    }
+
+    return out

pvlib/tests/test_iam.py

@@ -322,3 +322,51 @@ def test_marion_integrate_invalid():
 
     with pytest.raises(ValueError):
         _iam.marion_integrate(_iam.ashrae, 0, 'bad', 180)
+
+
+def test_fedis_defaults():
+    # expected values generated with code from model authors
+    surface_tilt = np.array([0, 30, 60, 90])
+    aoi = np.array([0, 30, 60, 90])
+
+    expected = {
+        'direct': np.array([0.994592623, 0.99301513, 0.943633517, 0.0]),
+        'sky': np.array([0.9410498, 0.95449355, 0.95357626, 0.9410498]),
+        'ground': np.array([0.0, 0.75380253, 0.9034704, 0.9410498]),
+    }
+    # numpy arrays
+    actual = _iam.fedis(aoi, surface_tilt)
+    assert_allclose(expected['direct'], actual['direct'], atol=1e-6)
+    assert_allclose(expected['sky'], actual['sky'])
+    assert_allclose(expected['ground'], actual['ground'], atol=1e-6)
+
+    # scalars
+    for i in range(len(aoi)):
+        actual = _iam.fedis(aoi[i], surface_tilt[i])
+        assert_allclose(expected['direct'][i], actual['direct'], atol=1e-6)
+        assert_allclose(expected['sky'][i], actual['sky'])
+        assert_allclose(expected['ground'][i], actual['ground'], atol=1e-6)
+
+
+def test_fedis_kwargs():
+    # expected values generated with code from model authors
+    surface_tilt = np.array([0, 30, 60, 90])
+    aoi = np.array([0, 30, 60, 90])
+
+    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)
+
+    # scalars
+    for i in range(len(aoi)):
+        actual = _iam.fedis(aoi[i], surface_tilt[i], n=1.7, n_ref=1.3)
+        assert_allclose(expected['direct'][i], actual['direct'], atol=1e-6)
+        assert_allclose(expected['sky'][i], actual['sky'])
+        assert_allclose(expected['ground'][i], actual['ground'], atol=1e-6)