Implement iam.schlick, iam.schlick_diffuse

docs/sphinx/source/reference/pv_modeling.rst

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

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

@@ -10,6 +10,11 @@ Deprecations
 Enhancements
 ~~~~~~~~~~~~
 * Multiple code style issues fixed that were reported by LGTM analysis. (:issue:`1275`, :pull:`1559`)
+* Added two direct IAM models :py:func:`pvlib.iam.schlick` and :py:func:`pvlib.iam.fedis`,
+  both of which can be used with :py:func:`~pvlib.iam.marion_diffuse` (:pull:`1562`, :issue:`1564`)
+* Added two diffuse IAM models :py:func:`pvlib.iam.schlick_diffuse` and 
+  :py:func:`pvlib.iam.fedis_diffuse` (:pull:`1562`)
+
 
 Bug fixes
 ~~~~~~~~~
@@ -34,4 +39,8 @@ Requirements
 
 Contributors
 ~~~~~~~~~~~~
-* Christian Orner (:ghuser:`chrisorner`)
+* Christian Orner (:ghuser:`chrisorner`)
+* Yu Xie (:ghuser:`xieyupku`)
+* Anton Driesse (:ghuser:`adriesse`)
+* Cliff Hansen (:ghuser:`cwhanse`)
+* Kevin Anderson (:ghuser:`kanderso-nrel`)

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', 'fedis'`.
 
     surface_tilt : numeric
         Surface tilt angles in decimal degrees.
@@ -592,6 +592,8 @@ def marion_diffuse(model, surface_tilt, **kwargs):
         'ashrae': ashrae,
         'sapm': sapm,
         'martin_ruiz': martin_ruiz,
+        'schlick': schlick,
+        'fedis': fedis,
     }
 
     try:
@@ -748,3 +750,255 @@ 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. Angles of nan will result in nan. [degrees]
+
+    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
+    pvlib.iam.schlick_diffuse
+    """
+    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 schlick_diffuse(surface_tilt):
+    """
+    Determine the incidence angle modifiers (IAM) for diffuse sky and
+    ground-reflected irradiance on a tilted surface using the Schlick
+    incident angle model.
+
+    The diffuse iam values are calculated using an analytical integration
+    of the Schlick equation [1]_ over the portion of an isotropic sky and
+    isotropic foreground that is visible from the tilted surface [2]_.
+
+    Parameters
+    ----------
+    surface_tilt : numeric
+        Surface tilt angle measured from horizontal (e.g. surface facing
+        up = 0, surface facing horizon = 90). [degrees]
+
+    Returns
+    -------
+    iam_sky : numeric
+        The incident angle modifier for sky diffuse
+
+    iam_ground : numeric
+        The incident angle modifier for ground-reflected diffuse
+
+    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.schlick
+    pvlib.iam.fedis_diffuse
+    """
+    # these calculations are as in [2]_, but with the refractive index
+    # weighting coefficient w set to 1.0 (so it is omitted)
+
+    # relative transmittance of sky diffuse radiation by PV cover:
+    cosB = cosd(surface_tilt)
+    sinB = sind(surface_tilt)
+    cuk = (2 / (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 in [2]
+
+    # relative transmittance of ground-reflected radiation by PV cover:
+    with np.errstate(divide='ignore', invalid='ignore'):  # Eq 6 in [2]
+        cug = 40 / (21 * (1 - cosB)) - (1 + cosB) / (1 - cosB) * cuk
+
+    cug = np.where(surface_tilt < 1e-6, 0, cug)
+
+    # respect input types:
+    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)
+
+    return cuk, cug
+
+
+def fedis(aoi, n=1.5, n_ref=None):
+    """
+    Determine the incidence angle modifier (IAM) for direct irradiance
+    using the FEDIS transmittance model.
+
+    Note that FEDIS [1]_ calculates direct IAM using the Fresnel equations
+    without extinction, thus in the default case of ``n_ref=None``, this is
+    equivalent to calling :py:func:`physical` with ``K=0``.
+
+    Parameters
+    ----------
+    aoi : numeric
+        Angle of incidence. [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. If None (default), set equal to ``n``.
+        In [1]_ ``n_ref`` was set to 1.4585 for
+        a fused silica dome over a CMP22, but in conventional
+        PV applications it is appropriate to set ``n_ref=n`` (the default
+        behavior).
+
+    Returns
+    -------
+    iam : numeric
+        The incident angle modifier
+
+    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`
+
+    See Also
+    --------
+    pvlib.iam.physical
+    """
+    iam_physical = physical(aoi, n, K=0, L=0)
+
+    if n_ref is None:
+        return iam_physical
+
+    else:
+        # reflectance for normal incidence with refractive index n:
+        rd0 = ((n - 1) / (n + 1)) ** 2
+
+        # reflectance for normal incidence with refractive index n_ref:
+        r0 = ((n_ref - 1) / (n_ref + 1)) ** 2
+
+        # weighting function
+        w = (1 - rd0) / (1 - r0)
+
+        return w * iam_physical
+
+
+def fedis_diffuse(surface_tilt, n=1.5, n_ref=None):
+    """
+    Determine the incidence angle modifiers (IAM) for diffuse sky and
+    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
+    ----------
+    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. If None (default), set equal to ``n``.
+        In [1]_ ``n_ref`` was set to 1.4585 for
+        a fused silica dome over a CMP22, but in conventional
+        PV applications it is appropriate to set ``n_ref=n`` (the default
+        behavior).
+
+    Returns
+    -------
+    iam_sky : numeric
+        The incident angle modifier for sky diffuse
+
+    iam_ground : numeric
+        The incident angle modifier for ground-reflected diffuse
+
+    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_diffuse
+    """
+    if n_ref is None:
+        n_ref = n
+
+    cuk, cug = schlick_diffuse(surface_tilt)
+
+    # 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
+
+    cuk = cuk * w
+    cug = cug * w
+
+    return cuk, cug