From 70c0493236c65276ac18d0e77dc42827501e43c2 Mon Sep 17 00:00:00 2001 From: Karel De Brabandere Date: Sun, 4 Dec 2022 21:11:59 +0100 Subject: [PATCH 1/7] [ENH] extend iam.physical with optional n_ar parameter for AR coating --- pvlib/iam.py | 111 +++++++++++++++++++++++----------------- pvlib/tests/test_iam.py | 18 ++++++- 2 files changed, 81 insertions(+), 48 deletions(-) diff --git a/pvlib/iam.py b/pvlib/iam.py index dfad91b6ff..ebac44b88e 100644 --- a/pvlib/iam.py +++ b/pvlib/iam.py @@ -11,7 +11,7 @@ import numpy as np import pandas as pd import functools -from pvlib.tools import cosd, sind, tand, asind +from pvlib.tools import cosd, sind # a dict of required parameter names for each IAM model # keys are the function names for the IAM models @@ -91,21 +91,20 @@ def ashrae(aoi, b=0.05): return iam -def physical(aoi, n=1.526, K=4., L=0.002): +def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None): r""" Determine the incidence angle modifier using refractive index ``n``, extinction coefficient ``K``, and glazing thickness ``L``. ``iam.physical`` calculates the incidence angle modifier as described in - [1]_, Section 3. The calculation is based on a physical model of absorbtion + [1]_, Section 3. The calculation is based on a physical model of absorption and transmission through a transparent cover. Parameters ---------- aoi : numeric The angle of incidence between the module normal vector and the - sun-beam vector in degrees. Angles of 0 are replaced with 1e-06 - to ensure non-nan results. Angles of nan will result in nan. + sun-beam vector in degrees. Angles of nan will result in nan. n : numeric, default 1.526 The effective index of refraction (unitless). Reference [1]_ @@ -121,6 +120,11 @@ def physical(aoi, n=1.526, K=4., L=0.002): indicates that 0.002 meters (2 mm) is reasonable for most glass-covered PV panels. + n_ar : numeric, default None + The effective index of refraction of the anti-reflective coating + (unitless). + If n_ar is None, no anti-reflective coating is applied. + Returns ------- iam : numeric @@ -149,48 +153,61 @@ def physical(aoi, n=1.526, K=4., L=0.002): pvlib.iam.interp pvlib.iam.sapm """ - zeroang = 1e-06 - - # hold a new reference to the input aoi object since we're going to - # overwrite the aoi reference below, but we'll need it for the - # series check at the end of the function - aoi_input = aoi - - aoi = np.where(aoi == 0, zeroang, aoi) - - # angle of reflection - thetar_deg = asind(1.0 / n * (sind(aoi))) - - # reflectance and transmittance for normal incidence light - rho_zero = ((1-n) / (1+n)) ** 2 - tau_zero = np.exp(-K*L) - - # reflectance for parallel and perpendicular polarized light - rho_para = (tand(thetar_deg - aoi) / tand(thetar_deg + aoi)) ** 2 - rho_perp = (sind(thetar_deg - aoi) / sind(thetar_deg + aoi)) ** 2 - - # transmittance for non-normal light - tau = np.exp(-K * L / cosd(thetar_deg)) - - # iam is ratio of non-normal to normal incidence transmitted light - # after deducting the reflected portion of each - iam = ((1 - (rho_para + rho_perp) / 2) / (1 - rho_zero) * tau / tau_zero) - - with np.errstate(invalid='ignore'): - # angles near zero produce nan, but iam is defined as one - small_angle = 1e-06 - iam = np.where(np.abs(aoi) < small_angle, 1.0, iam) - - # angles at 90 degrees can produce tiny negative values, - # which should be zero. this is a result of calculation precision - # rather than the physical model - iam = np.where(iam < 0, 0, iam) - - # for light coming from behind the plane, none can enter the module - iam = np.where(aoi > 90, 0, iam) - - if isinstance(aoi_input, pd.Series): - iam = pd.Series(iam, index=aoi_input.index) + n1, n3 = 1, n + if n_ar in [None, 1]: + # no AR coating + n2 = n + else: + n2 = n_ar + + # incidence angle + costheta = np.maximum(0, cosd(aoi)) # always >= 0 + sintheta = np.sqrt(1 - costheta**2) # always >= 0 + n1costheta1 = n1 * costheta + n2costheta1 = n2 * costheta + + # refraction angle of first interface + sintheta = n1 / n2 * sintheta + costheta = np.sqrt(1 - sintheta**2) + n1costheta2 = n1 * costheta + n2costheta2 = n2 * costheta + + # reflectance of s-, p-polarized, and normal light by the first interface + rho12_s = ((n1costheta1 - n2costheta2) / (n1costheta1 + n2costheta2)) ** 2 + rho12_p = ((n1costheta2 - n2costheta1) / (n1costheta2 + n2costheta1)) ** 2 + rho12_0 = ((n1 - n2) / (n1 + n2)) ** 2 + + # transmittance of s-, p-polarized, and normal light through the first interface + tau_s = 1 - rho12_s + tau_p = 1 - rho12_p + tau_0 = 1 - rho12_0 + + if n3 != n2: # AR coated glass + n3costheta2 = n3 * costheta + # refraction angle of second interface + sintheta = n2 / n3 * sintheta + costheta = np.sqrt(1 - sintheta**2) + n2costheta3 = n2 * costheta + n3costheta3 = n3 * costheta + + # reflectance of s-, p-polarized, and normal light by the second interface + rho23_s = ((n2costheta2 - n3costheta3) / (n2costheta2 + n3costheta3)) ** 2 + rho23_p = ((n2costheta3 - n3costheta2) / (n2costheta3 + n3costheta2)) ** 2 + rho23_0 = ((n2 - n3) / (n2 + n3)) ** 2 + + # transmittance through the coating, including internal reflections + # 1 + rho23 * rho12 + (rho23 * rho12)**2 + ... = 1 / (1 - rho23 * rho12) + tau_s *= (1 - rho23_s) / (1 - rho23_s * rho12_s) + tau_p *= (1 - rho23_p) / (1 - rho23_p * rho12_p) + tau_0 *= (1 - rho23_0) / (1 - rho23_0 * rho12_0) + + # transmittance of s-, p-polarized, and normal light after absorption in the glass + tau_s *= np.exp(-K * L / costheta) + tau_p *= np.exp(-K * L / costheta) + tau_0 *= np.exp(-K * L) + + # incidence angle modifier + iam = (tau_s + tau_p) / 2 / tau_0 return iam diff --git a/pvlib/tests/test_iam.py b/pvlib/tests/test_iam.py index df4d9ee877..eba1c66cb0 100644 --- a/pvlib/tests/test_iam.py +++ b/pvlib/tests/test_iam.py @@ -42,7 +42,7 @@ def test_physical(): expected = np.array([0, 0.8893998, 0.98797788, 0.99926198, 1, 0.99926198, 0.98797788, 0.8893998, 0, np.nan]) iam = _iam.physical(aoi, 1.526, 0.002, 4) - assert_allclose(iam, expected, equal_nan=True) + assert_allclose(iam, expected, atol=1e-7, equal_nan=True) # GitHub issue 397 aoi = pd.Series(aoi) @@ -51,6 +51,22 @@ def test_physical(): assert_series_equal(iam, expected) +def test_physical_ar(): + aoi = np.array([0, 22.5, 45, 67.5, 90, 100, np.nan]) + expected = np.array([1, 0.99944171, 0.9917463, 0.91506158, 0, 0, np.nan]) + iam = _iam.physical(aoi, n_ar=1.29) + assert_allclose(iam, expected, atol=1e-7, equal_nan=True) + + +def test_physical_noar(): + aoi = np.array([0, 22.5, 45, 67.5, 90, 100, np.nan]) + expected = _iam.physical(aoi) + iam0 = _iam.physical(aoi, n_ar=1) + iam1 = _iam.physical(aoi, n_ar=1.526) + assert_allclose(iam0, expected, equal_nan=True) + assert_allclose(iam1, expected, equal_nan=True) + + def test_physical_scalar(): aoi = -45. iam = _iam.physical(aoi, 1.526, 0.002, 4) From e54ff70fd4d0723eaba370010ba47b891e801651 Mon Sep 17 00:00:00 2001 From: Karel De Brabandere Date: Wed, 1 Feb 2023 20:25:49 +0100 Subject: [PATCH 2/7] [CLN] resolve stickler-ci failures; [ENH] ensure inputs could be vectors --- pvlib/iam.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/pvlib/iam.py b/pvlib/iam.py index ebac44b88e..d546ccc9c3 100644 --- a/pvlib/iam.py +++ b/pvlib/iam.py @@ -154,7 +154,7 @@ def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None): pvlib.iam.sapm """ n1, n3 = 1, n - if n_ar in [None, 1]: + if n_ar is None or np.allclose(n_ar, n1): # no AR coating n2 = n else: @@ -177,12 +177,12 @@ def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None): rho12_p = ((n1costheta2 - n2costheta1) / (n1costheta2 + n2costheta1)) ** 2 rho12_0 = ((n1 - n2) / (n1 + n2)) ** 2 - # transmittance of s-, p-polarized, and normal light through the first interface + # transmittance through the first interface tau_s = 1 - rho12_s tau_p = 1 - rho12_p tau_0 = 1 - rho12_0 - if n3 != n2: # AR coated glass + if not np.allclose(n3, n2): # AR coated glass n3costheta2 = n3 * costheta # refraction angle of second interface sintheta = n2 / n3 * sintheta @@ -190,18 +190,22 @@ def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None): n2costheta3 = n2 * costheta n3costheta3 = n3 * costheta - # reflectance of s-, p-polarized, and normal light by the second interface - rho23_s = ((n2costheta2 - n3costheta3) / (n2costheta2 + n3costheta3)) ** 2 - rho23_p = ((n2costheta3 - n3costheta2) / (n2costheta3 + n3costheta2)) ** 2 + # reflectance by the second interface + rho23_s = ( + (n2costheta2 - n3costheta3) / (n2costheta2 + n3costheta3) + ) ** 2 + rho23_p = ( + (n2costheta3 - n3costheta2) / (n2costheta3 + n3costheta2) + ) ** 2 rho23_0 = ((n2 - n3) / (n2 + n3)) ** 2 # transmittance through the coating, including internal reflections - # 1 + rho23 * rho12 + (rho23 * rho12)**2 + ... = 1 / (1 - rho23 * rho12) + # 1 + rho23*rho12 + (rho23*rho12)^2 + ... = 1/(1 - rho23*rho12) tau_s *= (1 - rho23_s) / (1 - rho23_s * rho12_s) tau_p *= (1 - rho23_p) / (1 - rho23_p * rho12_p) tau_0 *= (1 - rho23_0) / (1 - rho23_0 * rho12_0) - # transmittance of s-, p-polarized, and normal light after absorption in the glass + # transmittance after absorption in the glass tau_s *= np.exp(-K * L / costheta) tau_p *= np.exp(-K * L / costheta) tau_0 *= np.exp(-K * L) From f3a82dd3cd6aa43e326ae7fdde0d51bdf041a6cb Mon Sep 17 00:00:00 2001 From: Karel De Brabandere Date: Thu, 2 Feb 2023 08:35:50 +0100 Subject: [PATCH 3/7] [DOC] iam.physical: update docstring and add whatsnew entry --- docs/sphinx/source/whatsnew/v0.9.5.rst | 2 ++ pvlib/iam.py | 6 ++++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/docs/sphinx/source/whatsnew/v0.9.5.rst b/docs/sphinx/source/whatsnew/v0.9.5.rst index 816f41dfa7..c940adcd54 100644 --- a/docs/sphinx/source/whatsnew/v0.9.5.rst +++ b/docs/sphinx/source/whatsnew/v0.9.5.rst @@ -17,6 +17,8 @@ Deprecations Enhancements ~~~~~~~~~~~~ +* Added optional ``n_ar`` parameter to :py:func:`pvlib.iam.physical` to + support an anti-reflective coating. (:issue:`1501`, :pull:`1616`) Bug fixes ~~~~~~~~~ diff --git a/pvlib/iam.py b/pvlib/iam.py index d546ccc9c3..30fdbfccdb 100644 --- a/pvlib/iam.py +++ b/pvlib/iam.py @@ -94,10 +94,12 @@ def ashrae(aoi, b=0.05): def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None): r""" Determine the incidence angle modifier using refractive index ``n``, - extinction coefficient ``K``, and glazing thickness ``L``. + extinction coefficient ``K``, glazing thickness ``L`` and refractive + index ``n_ar`` of an eventual anti-reflective coating. ``iam.physical`` calculates the incidence angle modifier as described in - [1]_, Section 3. The calculation is based on a physical model of absorption + [1]_, Section 3, with additional support of an anti-reflective coating. + The calculation is based on a physical model of reflections, absorption, and transmission through a transparent cover. Parameters From c0ba45a7f6642ae7e3ad47dbb33029a7e227ada6 Mon Sep 17 00:00:00 2001 From: Karel De Brabandere Date: Thu, 2 Feb 2023 08:52:58 +0100 Subject: [PATCH 4/7] [CLN] minor non-functional change to improve readability --- pvlib/iam.py | 144 +++++++++++++++++++++++++++++---------------------- 1 file changed, 82 insertions(+), 62 deletions(-) diff --git a/pvlib/iam.py b/pvlib/iam.py index 30fdbfccdb..027da7681c 100644 --- a/pvlib/iam.py +++ b/pvlib/iam.py @@ -16,11 +16,11 @@ # a dict of required parameter names for each IAM model # keys are the function names for the IAM models _IAM_MODEL_PARAMS = { - 'ashrae': {'b'}, - 'physical': {'n', 'K', 'L'}, - 'martin_ruiz': {'a_r'}, - 'sapm': {'B0', 'B1', 'B2', 'B3', 'B4', 'B5'}, - 'interp': set() + "ashrae": {"b"}, + "physical": {"n", "K", "L"}, + "martin_ruiz": {"a_r"}, + "sapm": {"B0", "B1", "B2", "B3", "B4", "B5"}, + "interp": set(), } @@ -80,7 +80,7 @@ def ashrae(aoi, b=0.05): """ iam = 1 - b * (1 / np.cos(np.radians(aoi)) - 1) - aoi_gte_90 = np.full_like(aoi, False, dtype='bool') + aoi_gte_90 = np.full_like(aoi, False, dtype="bool") np.greater_equal(np.abs(aoi), 90, where=~np.isnan(aoi), out=aoi_gte_90) iam = np.where(aoi_gte_90, 0, iam) iam = np.maximum(0, iam) @@ -155,12 +155,14 @@ def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None): pvlib.iam.interp pvlib.iam.sapm """ - n1, n3 = 1, n + n1 = 1 if n_ar is None or np.allclose(n_ar, n1): # no AR coating - n2 = n + n2 = n3 = n else: + # AR coating n2 = n_ar + n3 = n # incidence angle costheta = np.maximum(0, cosd(aoi)) # always >= 0 @@ -219,7 +221,7 @@ def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None): def martin_ruiz(aoi, a_r=0.16): - r''' + r""" Determine the incidence angle modifier (IAM) using the Martin and Ruiz incident angle model. @@ -277,7 +279,7 @@ def martin_ruiz(aoi, a_r=0.16): pvlib.iam.ashrae pvlib.iam.interp pvlib.iam.sapm - ''' + """ # Contributed by Anton Driesse (@adriesse), PV Performance Labs. July, 2019 aoi_input = aoi @@ -288,7 +290,7 @@ def martin_ruiz(aoi, a_r=0.16): if np.any(np.less_equal(a_r, 0)): raise ValueError("The parameter 'a_r' cannot be zero or negative.") - with np.errstate(invalid='ignore'): + 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) @@ -299,7 +301,7 @@ def martin_ruiz(aoi, a_r=0.16): def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None): - ''' + """ Determine the incidence angle modifiers (iam) for diffuse sky and ground-reflected irradiance using the Martin and Ruiz incident angle model. @@ -363,7 +365,7 @@ def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None): pvlib.iam.ashrae pvlib.iam.interp pvlib.iam.sapm - ''' + """ # Contributed by Anton Driesse (@adriesse), PV Performance Labs. Oct. 2019 if isinstance(surface_tilt, pd.Series): @@ -389,25 +391,27 @@ def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None): cos = np.cos # avoid RuntimeWarnings for <, sin, and cos with nan - with np.errstate(invalid='ignore'): + with np.errstate(invalid="ignore"): # because sin(pi) isn't exactly zero sin_beta = np.where(surface_tilt < 90, sin(beta), sin(pi - beta)) trig_term_sky = sin_beta + (pi - beta - sin_beta) / (1 + cos(beta)) - trig_term_gnd = sin_beta + (beta - sin_beta) / (1 - cos(beta)) # noqa: E222 E261 E501 + trig_term_gnd = sin_beta + (beta - sin_beta) / ( + 1 - cos(beta) + ) # noqa: E222 E261 E501 iam_sky = 1 - np.exp(-(c1 + c2 * trig_term_sky) * trig_term_sky / a_r) iam_gnd = 1 - np.exp(-(c1 + c2 * trig_term_gnd) * trig_term_gnd / a_r) if out_index is not None: - iam_sky = pd.Series(iam_sky, index=out_index, name='iam_sky') - iam_gnd = pd.Series(iam_gnd, index=out_index, name='iam_ground') + iam_sky = pd.Series(iam_sky, index=out_index, name="iam_sky") + iam_gnd = pd.Series(iam_gnd, index=out_index, name="iam_ground") return iam_sky, iam_gnd -def interp(aoi, theta_ref, iam_ref, method='linear', normalize=True): - r''' +def interp(aoi, theta_ref, iam_ref, method="linear", normalize=True): + r""" Determine the incidence angle modifier (IAM) by interpolating a set of reference values, which are usually measured values. @@ -453,24 +457,29 @@ def interp(aoi, theta_ref, iam_ref, method='linear', normalize=True): pvlib.iam.ashrae pvlib.iam.martin_ruiz pvlib.iam.sapm - ''' + """ # 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} + 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) + 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') + 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) @@ -538,13 +547,19 @@ def sapm(aoi, module, upper=None): pvlib.iam.interp """ - aoi_coeff = [module['B5'], module['B4'], module['B3'], module['B2'], - module['B1'], module['B0']] + aoi_coeff = [ + module["B5"], + module["B4"], + module["B3"], + module["B2"], + module["B1"], + module["B0"], + ] iam = np.polyval(aoi_coeff, aoi) iam = np.clip(iam, 0, upper) # nan tolerant masking - aoi_lt_0 = np.full_like(aoi, False, dtype='bool') + aoi_lt_0 = np.full_like(aoi, False, dtype="bool") np.less(aoi, 0, where=~np.isnan(aoi), out=aoi_lt_0) iam = np.where(aoi_lt_0, 0, iam) @@ -610,21 +625,21 @@ def marion_diffuse(model, surface_tilt, **kwargs): """ models = { - 'physical': physical, - 'ashrae': ashrae, - 'sapm': sapm, - 'martin_ruiz': martin_ruiz, - 'schlick': schlick, + "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()))) + raise ValueError("model must be one of: " + str(list(models.keys()))) iam_function = functools.partial(iam_model, **kwargs) iam = {} - for region in ['sky', 'horizon', 'ground']: + for region in ["sky", "horizon", "ground"]: iam[region] = marion_integrate(iam_function, surface_tilt, region) return iam @@ -695,34 +710,34 @@ def marion_integrate(function, surface_tilt, region, num=None): """ if num is None: - if region in ['sky', 'ground']: + if region in ["sky", "ground"]: num = 180 - elif region == 'horizon': + elif region == "horizon": num = 1800 else: - raise ValueError(f'Invalid region: {region}') + raise ValueError(f"Invalid region: {region}") beta = np.radians(surface_tilt) if isinstance(beta, pd.Series): # convert Series to np array for broadcasting later beta = beta.values - ai = np.pi/num # angular increment + ai = np.pi / num # angular increment phi_range = np.linspace(0, np.pi, num, endpoint=False) - psi_range = np.linspace(0, 2*np.pi, 2*num, endpoint=False) + psi_range = np.linspace(0, 2 * np.pi, 2 * num, endpoint=False) # the pseudocode in [1] do these checks at the end, but it's # faster to do this criteria check up front instead of later. - if region == 'sky': - mask = phi_range + ai <= np.pi/2 - elif region == 'horizon': - lo = 89.5 * np.pi/180 - hi = np.pi/2 + if region == "sky": + mask = phi_range + ai <= np.pi / 2 + elif region == "horizon": + lo = 89.5 * np.pi / 180 + hi = np.pi / 2 mask = (lo <= phi_range) & (phi_range + ai <= hi) - elif region == 'ground': - mask = (phi_range >= np.pi/2) + elif region == "ground": + mask = phi_range >= np.pi / 2 else: - raise ValueError(f'Invalid region: {region}') + raise ValueError(f"Invalid region: {region}") phi_range = phi_range[mask] # fast Cartesian product of phi and psi @@ -733,8 +748,8 @@ def marion_integrate(function, surface_tilt, region, num=None): psi_1 = angles[:, [1]] phi_2 = phi_1 + ai # psi_2 = psi_1 + ai # not needed - phi_avg = phi_1 + 0.5*ai - psi_avg = psi_1 + 0.5*ai + phi_avg = phi_1 + 0.5 * ai + psi_avg = psi_1 + 0.5 * ai term_1 = np.cos(beta) * np.cos(phi_avg) # The AOI formula includes a term based on the difference between # panel azimuth and the photon azimuth, but because we assume each class @@ -751,18 +766,18 @@ def marion_integrate(function, surface_tilt, region, num=None): dAs = ai * (np.cos(phi_1) - np.cos(phi_2)) cosaoi_dAs = cosaoi * dAs # apply the final AOI check, zeroing out non-passing points - mask = aoi < np.pi/2 + mask = aoi < np.pi / 2 cosaoi_dAs = np.where(mask, cosaoi_dAs, 0) numerator = np.sum(function(np.degrees(aoi)) * cosaoi_dAs, axis=0) denominator = np.sum(cosaoi_dAs, axis=0) - with np.errstate(invalid='ignore'): + with np.errstate(invalid="ignore"): # in some cases, no points pass the criteria # (e.g. region='ground', surface_tilt=0), so we override the division # by zero to set Fd=0. Also, preserve nans in beta. - Fd = np.where((denominator != 0) | ~np.isfinite(beta), - numerator / denominator, - 0) + Fd = np.where( + (denominator != 0) | ~np.isfinite(beta), numerator / denominator, 0 + ) # preserve input type if np.isscalar(surface_tilt): @@ -871,13 +886,18 @@ def schlick_diffuse(surface_tilt): 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 + (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] + 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) From ab07267612c9531df6bdd376ce22e39086564f17 Mon Sep 17 00:00:00 2001 From: Karel De Brabandere Date: Thu, 2 Feb 2023 12:51:37 +0100 Subject: [PATCH 5/7] revert latest commit --- pvlib/iam.py | 144 ++++++++++++++++++++++----------------------------- 1 file changed, 62 insertions(+), 82 deletions(-) diff --git a/pvlib/iam.py b/pvlib/iam.py index 027da7681c..30fdbfccdb 100644 --- a/pvlib/iam.py +++ b/pvlib/iam.py @@ -16,11 +16,11 @@ # a dict of required parameter names for each IAM model # keys are the function names for the IAM models _IAM_MODEL_PARAMS = { - "ashrae": {"b"}, - "physical": {"n", "K", "L"}, - "martin_ruiz": {"a_r"}, - "sapm": {"B0", "B1", "B2", "B3", "B4", "B5"}, - "interp": set(), + 'ashrae': {'b'}, + 'physical': {'n', 'K', 'L'}, + 'martin_ruiz': {'a_r'}, + 'sapm': {'B0', 'B1', 'B2', 'B3', 'B4', 'B5'}, + 'interp': set() } @@ -80,7 +80,7 @@ def ashrae(aoi, b=0.05): """ iam = 1 - b * (1 / np.cos(np.radians(aoi)) - 1) - aoi_gte_90 = np.full_like(aoi, False, dtype="bool") + aoi_gte_90 = np.full_like(aoi, False, dtype='bool') np.greater_equal(np.abs(aoi), 90, where=~np.isnan(aoi), out=aoi_gte_90) iam = np.where(aoi_gte_90, 0, iam) iam = np.maximum(0, iam) @@ -155,14 +155,12 @@ def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None): pvlib.iam.interp pvlib.iam.sapm """ - n1 = 1 + n1, n3 = 1, n if n_ar is None or np.allclose(n_ar, n1): # no AR coating - n2 = n3 = n + n2 = n else: - # AR coating n2 = n_ar - n3 = n # incidence angle costheta = np.maximum(0, cosd(aoi)) # always >= 0 @@ -221,7 +219,7 @@ def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None): def martin_ruiz(aoi, a_r=0.16): - r""" + r''' Determine the incidence angle modifier (IAM) using the Martin and Ruiz incident angle model. @@ -279,7 +277,7 @@ def martin_ruiz(aoi, a_r=0.16): pvlib.iam.ashrae pvlib.iam.interp pvlib.iam.sapm - """ + ''' # Contributed by Anton Driesse (@adriesse), PV Performance Labs. July, 2019 aoi_input = aoi @@ -290,7 +288,7 @@ def martin_ruiz(aoi, a_r=0.16): if np.any(np.less_equal(a_r, 0)): raise ValueError("The parameter 'a_r' cannot be zero or negative.") - with np.errstate(invalid="ignore"): + 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) @@ -301,7 +299,7 @@ def martin_ruiz(aoi, a_r=0.16): def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None): - """ + ''' Determine the incidence angle modifiers (iam) for diffuse sky and ground-reflected irradiance using the Martin and Ruiz incident angle model. @@ -365,7 +363,7 @@ def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None): pvlib.iam.ashrae pvlib.iam.interp pvlib.iam.sapm - """ + ''' # Contributed by Anton Driesse (@adriesse), PV Performance Labs. Oct. 2019 if isinstance(surface_tilt, pd.Series): @@ -391,27 +389,25 @@ def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None): cos = np.cos # avoid RuntimeWarnings for <, sin, and cos with nan - with np.errstate(invalid="ignore"): + with np.errstate(invalid='ignore'): # because sin(pi) isn't exactly zero sin_beta = np.where(surface_tilt < 90, sin(beta), sin(pi - beta)) trig_term_sky = sin_beta + (pi - beta - sin_beta) / (1 + cos(beta)) - trig_term_gnd = sin_beta + (beta - sin_beta) / ( - 1 - cos(beta) - ) # noqa: E222 E261 E501 + trig_term_gnd = sin_beta + (beta - sin_beta) / (1 - cos(beta)) # noqa: E222 E261 E501 iam_sky = 1 - np.exp(-(c1 + c2 * trig_term_sky) * trig_term_sky / a_r) iam_gnd = 1 - np.exp(-(c1 + c2 * trig_term_gnd) * trig_term_gnd / a_r) if out_index is not None: - iam_sky = pd.Series(iam_sky, index=out_index, name="iam_sky") - iam_gnd = pd.Series(iam_gnd, index=out_index, name="iam_ground") + iam_sky = pd.Series(iam_sky, index=out_index, name='iam_sky') + iam_gnd = pd.Series(iam_gnd, index=out_index, name='iam_ground') return iam_sky, iam_gnd -def interp(aoi, theta_ref, iam_ref, method="linear", normalize=True): - r""" +def interp(aoi, theta_ref, iam_ref, method='linear', normalize=True): + r''' Determine the incidence angle modifier (IAM) by interpolating a set of reference values, which are usually measured values. @@ -457,29 +453,24 @@ def interp(aoi, theta_ref, iam_ref, method="linear", normalize=True): pvlib.iam.ashrae pvlib.iam.martin_ruiz pvlib.iam.sapm - """ + ''' # 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} + 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 - ) + 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" - ) + 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) @@ -547,19 +538,13 @@ def sapm(aoi, module, upper=None): pvlib.iam.interp """ - aoi_coeff = [ - module["B5"], - module["B4"], - module["B3"], - module["B2"], - module["B1"], - module["B0"], - ] + aoi_coeff = [module['B5'], module['B4'], module['B3'], module['B2'], + module['B1'], module['B0']] iam = np.polyval(aoi_coeff, aoi) iam = np.clip(iam, 0, upper) # nan tolerant masking - aoi_lt_0 = np.full_like(aoi, False, dtype="bool") + aoi_lt_0 = np.full_like(aoi, False, dtype='bool') np.less(aoi, 0, where=~np.isnan(aoi), out=aoi_lt_0) iam = np.where(aoi_lt_0, 0, iam) @@ -625,21 +610,21 @@ def marion_diffuse(model, surface_tilt, **kwargs): """ models = { - "physical": physical, - "ashrae": ashrae, - "sapm": sapm, - "martin_ruiz": martin_ruiz, - "schlick": schlick, + '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()))) + raise ValueError('model must be one of: ' + str(list(models.keys()))) iam_function = functools.partial(iam_model, **kwargs) iam = {} - for region in ["sky", "horizon", "ground"]: + for region in ['sky', 'horizon', 'ground']: iam[region] = marion_integrate(iam_function, surface_tilt, region) return iam @@ -710,34 +695,34 @@ def marion_integrate(function, surface_tilt, region, num=None): """ if num is None: - if region in ["sky", "ground"]: + if region in ['sky', 'ground']: num = 180 - elif region == "horizon": + elif region == 'horizon': num = 1800 else: - raise ValueError(f"Invalid region: {region}") + raise ValueError(f'Invalid region: {region}') beta = np.radians(surface_tilt) if isinstance(beta, pd.Series): # convert Series to np array for broadcasting later beta = beta.values - ai = np.pi / num # angular increment + ai = np.pi/num # angular increment phi_range = np.linspace(0, np.pi, num, endpoint=False) - psi_range = np.linspace(0, 2 * np.pi, 2 * num, endpoint=False) + psi_range = np.linspace(0, 2*np.pi, 2*num, endpoint=False) # the pseudocode in [1] do these checks at the end, but it's # faster to do this criteria check up front instead of later. - if region == "sky": - mask = phi_range + ai <= np.pi / 2 - elif region == "horizon": - lo = 89.5 * np.pi / 180 - hi = np.pi / 2 + if region == 'sky': + mask = phi_range + ai <= np.pi/2 + elif region == 'horizon': + lo = 89.5 * np.pi/180 + hi = np.pi/2 mask = (lo <= phi_range) & (phi_range + ai <= hi) - elif region == "ground": - mask = phi_range >= np.pi / 2 + elif region == 'ground': + mask = (phi_range >= np.pi/2) else: - raise ValueError(f"Invalid region: {region}") + raise ValueError(f'Invalid region: {region}') phi_range = phi_range[mask] # fast Cartesian product of phi and psi @@ -748,8 +733,8 @@ def marion_integrate(function, surface_tilt, region, num=None): psi_1 = angles[:, [1]] phi_2 = phi_1 + ai # psi_2 = psi_1 + ai # not needed - phi_avg = phi_1 + 0.5 * ai - psi_avg = psi_1 + 0.5 * ai + phi_avg = phi_1 + 0.5*ai + psi_avg = psi_1 + 0.5*ai term_1 = np.cos(beta) * np.cos(phi_avg) # The AOI formula includes a term based on the difference between # panel azimuth and the photon azimuth, but because we assume each class @@ -766,18 +751,18 @@ def marion_integrate(function, surface_tilt, region, num=None): dAs = ai * (np.cos(phi_1) - np.cos(phi_2)) cosaoi_dAs = cosaoi * dAs # apply the final AOI check, zeroing out non-passing points - mask = aoi < np.pi / 2 + mask = aoi < np.pi/2 cosaoi_dAs = np.where(mask, cosaoi_dAs, 0) numerator = np.sum(function(np.degrees(aoi)) * cosaoi_dAs, axis=0) denominator = np.sum(cosaoi_dAs, axis=0) - with np.errstate(invalid="ignore"): + with np.errstate(invalid='ignore'): # in some cases, no points pass the criteria # (e.g. region='ground', surface_tilt=0), so we override the division # by zero to set Fd=0. Also, preserve nans in beta. - Fd = np.where( - (denominator != 0) | ~np.isfinite(beta), numerator / denominator, 0 - ) + Fd = np.where((denominator != 0) | ~np.isfinite(beta), + numerator / denominator, + 0) # preserve input type if np.isscalar(surface_tilt): @@ -886,18 +871,13 @@ def schlick_diffuse(surface_tilt): 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 + (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] + 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) From f2957638a1dcf896890c0a7d75ba0ded18aba949 Mon Sep 17 00:00:00 2001 From: Karel De Brabandere Date: Thu, 2 Feb 2023 13:26:51 +0100 Subject: [PATCH 6/7] apply suggested changes --- pvlib/iam.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pvlib/iam.py b/pvlib/iam.py index 30fdbfccdb..98e4290962 100644 --- a/pvlib/iam.py +++ b/pvlib/iam.py @@ -95,7 +95,7 @@ def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None): r""" Determine the incidence angle modifier using refractive index ``n``, extinction coefficient ``K``, glazing thickness ``L`` and refractive - index ``n_ar`` of an eventual anti-reflective coating. + index ``n_ar`` of an optional anti-reflective coating. ``iam.physical`` calculates the incidence angle modifier as described in [1]_, Section 3, with additional support of an anti-reflective coating. @@ -124,8 +124,7 @@ def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None): n_ar : numeric, default None The effective index of refraction of the anti-reflective coating - (unitless). - If n_ar is None, no anti-reflective coating is applied. + (unitless). If n_ar is None, no anti-reflective coating is applied. Returns ------- From 14cd046e1f10f6a9a293d682541c96b1c5c10b89 Mon Sep 17 00:00:00 2001 From: Karel De Brabandere Date: Sat, 4 Feb 2023 14:55:54 +0100 Subject: [PATCH 7/7] update docstring, add contributors --- docs/sphinx/source/whatsnew/v0.9.5.rst | 4 ++++ pvlib/iam.py | 7 ++++--- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/docs/sphinx/source/whatsnew/v0.9.5.rst b/docs/sphinx/source/whatsnew/v0.9.5.rst index c940adcd54..e4cec140ec 100644 --- a/docs/sphinx/source/whatsnew/v0.9.5.rst +++ b/docs/sphinx/source/whatsnew/v0.9.5.rst @@ -48,3 +48,7 @@ Contributors * Kevin Anderson (:ghuser:`kanderso-nrel`) * Will Holmgren (:ghuser:`wholmgren`) * Pratham Chauhan (:ghuser:`ooprathamm`) +* Karel De Brabandere (:ghuser:`kdebrab`) +* Mark Mikofski (:ghuser:`mikofski`) +* Anton Driesse (:ghuser:`adriesse`) +* Adam R. Jensen (:ghuser:`AdamRJensen`) diff --git a/pvlib/iam.py b/pvlib/iam.py index 98e4290962..3eaa6b4c8e 100644 --- a/pvlib/iam.py +++ b/pvlib/iam.py @@ -122,9 +122,10 @@ def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None): indicates that 0.002 meters (2 mm) is reasonable for most glass-covered PV panels. - n_ar : numeric, default None - The effective index of refraction of the anti-reflective coating - (unitless). If n_ar is None, no anti-reflective coating is applied. + n_ar : numeric, optional + The effective index of refraction of the anti-reflective (AR) coating + (unitless). If n_ar is None (default), no AR coating is applied. + A typical value for the effective index of an AR coating is 1.29. Returns -------