r"""
The ``iam`` module contains functions that implement models for the incidence
angle modifier (IAM). The IAM quantifies the fraction of direct irradiance on
a module's front surface that is transmitted through the module materials to
the cells. Stated differently, the quantity 1 - IAM is the fraction of direct
irradiance that is reflected away or absorbed by the module's front materials.
IAM is typically a function of the angle of incidence (AOI) of the direct
irradiance to the module's surface.
"""
import numpy as np
import pandas as pd
import functools
from scipy.optimize import minimize
from pvlib.tools import cosd, sind, acosd
# 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': {'theta_ref', 'iam_ref'}
}
[docs]
def ashrae(aoi, b=0.05):
r"""
Determine the incidence angle modifier using the ASHRAE transmission
model.
The ASHRAE (American Society of Heating, Refrigeration, and Air
Conditioning Engineers) transmission model is developed in
[1]_, and in [2]_. The model has been used in software such as PVSyst [3]_.
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.
b : float, default 0.05
A parameter to adjust the incidence angle modifier as a function of
angle of incidence. Typical values are on the order of 0.05 [3].
Returns
-------
iam : numeric
The incident angle modifier (IAM). Returns zero for all abs(aoi) >= 90
and for all ``iam`` values that would be less than 0.
Notes
-----
The incidence angle modifier is calculated as
.. math::
IAM = 1 - b (\sec(aoi) - 1)
As AOI approaches 90 degrees, the model yields negative values for IAM;
negative IAM values are set to zero in this implementation.
References
----------
.. [1] Souka A.F., Safwat H.H., "Determination of the optimum
orientations for the double exposure flat-plate collector and its
reflections". Solar Energy vol .10, pp 170-174. 1966.
.. [2] ASHRAE standard 93-77
.. [3] PVsyst 7 Help.
https://www.pvsyst.com/help/index.html?iam_loss.htm retrieved on
January 30, 2024
See Also
--------
pvlib.iam.physical
pvlib.iam.martin_ruiz
pvlib.iam.interp
"""
iam = 1 - b * (1 / np.cos(np.radians(aoi)) - 1)
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)
if isinstance(aoi, pd.Series):
iam = pd.Series(iam, index=aoi.index)
return iam
[docs]
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 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.
The calculation is based on a physical model of reflections, 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 nan will result in nan.
n : numeric, default 1.526
The effective index of refraction (unitless). Reference [1]_
indicates that a value of 1.526 is acceptable for glass.
K : numeric, default 4.0
The glazing extinction coefficient in units of 1/meters.
Reference [1] indicates that a value of 4 is reasonable for
"water white" glass.
L : numeric, default 0.002
The glazing thickness in units of meters. Reference [1]_
indicates that 0.002 meters (2 mm) is reasonable for most
glass-covered PV panels.
n_ar : numeric, optional
The effective index of refraction of the anti-reflective (AR) coating
(unitless). If ``n_ar`` is not supplied, no AR coating is applied.
A typical value for the effective index of an AR coating is 1.29.
Returns
-------
iam : numeric
The incident angle modifier
Notes
-----
The pvlib python authors believe that Eqn. 14 in [1]_ is
incorrect, which presents :math:`\theta_{r} = \arcsin(n \sin(AOI))`.
Here, :math:`\theta_{r} = \arcsin(1/n \times \sin(AOI))`
References
----------
.. [1] W. De Soto et al., "Improvement and validation of a model for
photovoltaic array performance", Solar Energy, vol 80, pp. 78-88,
2006.
.. [2] Duffie, John A. & Beckman, William A.. (2006). Solar Engineering
of Thermal Processes, third edition. [Books24x7 version] Available
from http://common.books24x7.com/toc.aspx?bookid=17160.
See Also
--------
pvlib.iam.martin_ruiz
pvlib.iam.ashrae
pvlib.iam.interp
pvlib.iam.sapm
"""
n1, n3 = 1, n
if n_ar is None or np.allclose(n_ar, n1):
# 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
with np.errstate(divide='ignore', invalid='ignore'):
rho12_s = \
((n1costheta1 - n2costheta2) / (n1costheta1 + n2costheta2)) ** 2
rho12_p = \
((n1costheta2 - n2costheta1) / (n1costheta2 + n2costheta1)) ** 2
rho12_0 = ((n1 - n2) / (n1 + n2)) ** 2
# transmittance through the first interface
tau_s = 1 - rho12_s
tau_p = 1 - rho12_p
tau_0 = 1 - rho12_0
if not np.allclose(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 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 after absorption in the glass
with np.errstate(divide='ignore', invalid='ignore'):
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
# for light coming from behind the plane, none can enter the module
# when n2 > 1, this is already the case
if np.isclose(n2, 1).any():
iam = np.where(aoi >= 90, 0, iam)
if isinstance(aoi, pd.Series):
iam = pd.Series(iam, index=aoi.index)
return iam
[docs]
def martin_ruiz(aoi, a_r=0.16):
r'''
Determine the incidence angle modifier (IAM) using the Martin
and Ruiz incident angle model.
Parameters
----------
aoi : numeric, degrees
The angle of incidence between the module normal vector and the
sun-beam vector in degrees.
a_r : numeric
The angular losses coefficient described in equation 3 of [1]_.
This is an empirical dimensionless parameter. Values of ``a_r`` are
generally on the order of 0.08 to 0.25 for flat-plate PV modules.
Returns
-------
iam : numeric
The incident angle modifier(s)
Notes
-----
`martin_ruiz` calculates the incidence angle modifier (IAM) as described in
[1]_. The information required is the incident angle (AOI) and the angular
losses coefficient (a_r). Note that [1]_ has a corrigendum [2]_ which
clarifies a mix-up of 'alpha's and 'a's in the former.
The incident angle modifier is defined as
.. math::
IAM = \frac{1 - \exp(-\frac{\cos(aoi)}{a_r})}
{1 - \exp(\frac{-1}{a_r})}
which is presented as :math:`AL(\alpha) = 1 - IAM` in equation 4 of [1]_,
with :math:`\alpha` representing the angle of incidence AOI. Thus IAM = 1
at AOI = 0°, and IAM = 0 at AOI = 90°. This equation is only valid for
0° <= aoi <= 90°, therefore `iam` is constrained to 0.0 outside this
interval.
References
----------
.. [1] N. Martin and J. M. Ruiz, "Calculation of the PV modules angular
losses under field conditions by means of an analytical model", Solar
Energy Materials & Solar Cells, vol. 70, pp. 25-38, 2001.
.. [2] N. Martin and J. M. Ruiz, "Corrigendum to 'Calculation of the PV
modules angular losses under field conditions by means of an
analytical model'", Solar Energy Materials & Solar Cells, vol. 110,
pp. 154, 2013.
See Also
--------
pvlib.iam.martin_ruiz_diffuse
pvlib.iam.physical
pvlib.iam.ashrae
pvlib.iam.interp
pvlib.iam.sapm
'''
# Contributed by Anton Driesse (@adriesse), PV Performance Labs. July, 2019
aoi_input = aoi
aoi = np.asanyarray(aoi)
a_r = np.asanyarray(a_r)
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'):
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)
if isinstance(aoi_input, pd.Series):
iam = pd.Series(iam, index=aoi_input.index)
return iam
[docs]
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.
Parameters
----------
surface_tilt: float or array-like, default 0
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)
surface_tilt must be in the range [0, 180]
a_r : numeric
The angular losses coefficient described in equation 3 of [1]_.
This is an empirical dimensionless parameter. Values of a_r are
generally on the order of 0.08 to 0.25 for flat-plate PV modules.
a_r must be greater than zero.
c1 : float
First fitting parameter for the expressions that approximate the
integral of diffuse irradiance coming from different directions.
c1 is given as the constant 4 / 3 / pi (0.4244) in [1]_.
c2 : float
Second fitting parameter for the expressions that approximate the
integral of diffuse irradiance coming from different directions.
If c2 is not specified, it will be calculated according to the linear
relationship given in [3]_.
Returns
-------
iam_sky : numeric
The incident angle modifier for sky diffuse
iam_ground : numeric
The incident angle modifier for ground-reflected diffuse
Notes
-----
Sky and ground modifiers are complementary: iam_sky for tilt = 30 is
equal to iam_ground for tilt = 180 - 30. For vertical surfaces,
tilt = 90, the two factors are equal.
References
----------
.. [1] N. Martin and J. M. Ruiz, "Calculation of the PV modules angular
losses under field conditions by means of an analytical model", Solar
Energy Materials & Solar Cells, vol. 70, pp. 25-38, 2001.
.. [2] N. Martin and J. M. Ruiz, "Corrigendum to 'Calculation of the PV
modules angular losses under field conditions by means of an
analytical model'", Solar Energy Materials & Solar Cells, vol. 110,
pp. 154, 2013.
.. [3] "IEC 61853-3 Photovoltaic (PV) module performance testing and energy
rating - Part 3: Energy rating of PV modules". IEC, Geneva, 2018.
See Also
--------
pvlib.iam.martin_ruiz
pvlib.iam.physical
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):
out_index = surface_tilt.index
else:
out_index = None
surface_tilt = np.asanyarray(surface_tilt)
# 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 == 180, 180 - zeroang, surface_tilt)
if c2 is None:
# This equation is from [3] Sect. 7.2
c2 = 0.5 * a_r - 0.154
beta = np.radians(surface_tilt)
sin = np.sin
pi = np.pi
cos = np.cos
# avoid RuntimeWarnings for <, sin, and cos with nan
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
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')
return iam_sky, iam_gnd
[docs]
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.
Parameters
----------
aoi : numeric
The angle of incidence between the module normal vector and the
sun-beam vector [degrees].
theta_ref : numeric
Vector of angles at which the IAM is known [degrees].
iam_ref : numeric
IAM values for each angle in ``theta_ref`` [unitless].
method : str, default 'linear'
Specifies the interpolation method.
Useful options are: 'linear', 'quadratic', 'cubic'.
See scipy.interpolate.interp1d for more options.
normalize : boolean, default True
When true, the interpolated values are divided by the interpolated
value at zero degrees. This ensures that ``iam=1.0`` at normal
incidence.
Returns
-------
iam : numeric
The incident angle modifier(s) [unitless]
Notes
-----
``theta_ref`` must have two or more points and may span any range of
angles. Typically there will be a dozen or more points in the range 0-90
degrees. Beyond the range of ``theta_ref``, IAM values are extrapolated,
but constrained to be non-negative.
The sign of ``aoi`` is ignored; only the magnitude is used.
See Also
--------
pvlib.iam.physical
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}
if len(theta_ref) < MIN_REF_VALS.get(method, 2):
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')
aoi_input = aoi
aoi = np.asanyarray(aoi)
aoi = np.abs(aoi)
iam = interpolator(aoi)
iam = np.clip(iam, 0, None)
if normalize:
iam /= interpolator(0)
if isinstance(aoi_input, pd.Series):
iam = pd.Series(iam, index=aoi_input.index)
return iam
[docs]
def sapm(aoi, module, upper=None):
r"""
Determine the incidence angle modifier (IAM) using the SAPM model.
Parameters
----------
aoi : numeric
Angle of incidence in degrees. Negative input angles will return
zeros.
module : dict-like
A dict or Series with the SAPM IAM model parameters.
See the :py:func:`sapm` notes section for more details.
upper : float, optional
Upper limit on the results.
Returns
-------
iam : numeric
The SAPM angle of incidence loss coefficient, termed F2 in [1]_.
Notes
-----
The SAPM [1]_ traditionally does not define an upper limit on the AOI
loss function and values slightly exceeding 1 may exist for moderate
angles of incidence (15-40 degrees). However, users may consider
imposing an upper limit of 1.
References
----------
.. [1] King, D. et al, 2004, "Sandia Photovoltaic Array Performance
Model", SAND Report 3535, Sandia National Laboratories, Albuquerque,
NM.
.. [2] B.H. King et al, "Procedure to Determine Coefficients for the
Sandia Array Performance Model (SAPM)," SAND2016-5284, Sandia
National Laboratories (2016).
.. [3] B.H. King et al, "Recent Advancements in Outdoor Measurement
Techniques for Angle of Incidence Effects," 42nd IEEE PVSC (2015).
:doi:`10.1109/PVSC.2015.7355849`
See Also
--------
pvlib.iam.physical
pvlib.iam.ashrae
pvlib.iam.martin_ruiz
pvlib.iam.interp
"""
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')
np.less(aoi, 0, where=~np.isnan(aoi), out=aoi_lt_0)
iam = np.where(aoi_lt_0, 0, iam)
if isinstance(aoi, pd.Series):
iam = pd.Series(iam, aoi.index)
return iam
[docs]
def marion_diffuse(model, surface_tilt, **kwargs):
"""
Determine diffuse irradiance incidence angle modifiers using Marion's
method of integrating over solid angle.
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).
**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_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`
Examples
--------
>>> marion_diffuse('physical', surface_tilt=20)
{'sky': 0.9539178294437575,
'horizon': 0.7652650139134007,
'ground': 0.6387140117795903}
>>> marion_diffuse('ashrae', [20, 30], b=0.04)
{'sky': array([0.96748999, 0.96938408]),
'horizon': array([0.86478428, 0.91825792]),
'ground': array([0.77004435, 0.8522436 ])}
"""
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())))
iam_function = functools.partial(iam_model, **kwargs)
iam = {}
for region in ['sky', 'horizon', 'ground']:
iam[region] = marion_integrate(iam_function, surface_tilt, region)
return iam
[docs]
def marion_integrate(function, surface_tilt, region, num=None):
"""
Integrate an incidence angle modifier (IAM) function over solid angle
to determine a diffuse irradiance correction factor using Marion's method.
This lower-level function actually performs the IAM integration for the
specified solid angle region.
Parameters
----------
function : callable(aoi)
The IAM function to evaluate across solid angle. The function must
be vectorized and take only one parameter, the angle of incidence in
degrees.
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).
region : {'sky', 'horizon', 'ground'}
The region to integrate over. Must be one of:
* '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.
num : int, optional
The number of increments in the zenith integration.
If not specified, N will follow the values used in [1]_:
* 'sky' or 'ground': num = 180
* 'horizon': num = 1800
Returns
-------
iam : numeric
AOI diffuse correction factor for the specified region.
See Also
--------
pvlib.iam.marion_diffuse
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`
Examples
--------
>>> marion_integrate(pvlib.iam.ashrae, 20, 'sky')
0.9596085829811408
>>> from functools import partial
>>> f = partial(pvlib.iam.physical, n=1.3)
>>> marion_integrate(f, [20, 30], 'sky')
array([0.96225034, 0.9653219 ])
"""
if num is None:
if region in ['sky', 'ground']:
num = 180
elif region == 'horizon':
num = 1800
else:
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
phi_range = np.linspace(0, np.pi, 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
mask = (lo <= phi_range) & (phi_range + ai <= hi)
elif region == 'ground':
mask = (phi_range >= np.pi/2)
else:
raise ValueError(f'Invalid region: {region}')
phi_range = phi_range[mask]
# fast Cartesian product of phi and psi
angles = np.array(np.meshgrid(phi_range, psi_range)).T.reshape(-1, 2)
# index with single-element lists to maintain 2nd dimension so that
# these angle arrays broadcast across the beta array
phi_1 = angles[:, [0]]
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
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
# of diffuse irradiance is isotropic and we are integrating over all
# angles, it doesn't matter what panel azimuth we choose (i.e., the
# system is rotationally invariant). So we choose gamma to be zero so
# that we can omit it from the cos(psi_avg) term.
# Marion's paper mentions this in the Section 3 pseudocode:
# "set gamma to pi (or any value between 0 and 2pi)"
term_2 = np.sin(beta) * np.sin(phi_avg) * np.cos(psi_avg)
cosaoi = term_1 + term_2
aoi = np.arccos(cosaoi)
# simplify Eq 8, (psi_2 - psi_1) is always ai
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
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'):
# 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)
# preserve input type
if np.isscalar(surface_tilt):
Fd = Fd.item()
elif isinstance(surface_tilt, pd.Series):
Fd = pd.Series(Fd, surface_tilt.index)
return Fd
[docs]
def schlick(aoi):
"""
Determine incidence angle modifier (IAM) for direct irradiance 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]_ (see :py:func:`schlick_diffuse`).
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.
See Also
--------
pvlib.iam.schlick_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`
"""
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
[docs]
def schlick_diffuse(surface_tilt):
r"""
Determine the incidence angle modifiers (IAM) for diffuse sky and
ground-reflected irradiance on a tilted surface using the Schlick
incident angle model.
The Schlick equation (or "Schlick's approximation") [1]_ is an
approximation to the Fresnel reflection factor which can be recast as
a simple photovoltaic IAM model like so:
.. math::
IAM = 1 - (1 - \cos(aoi))^5
Unlike the Fresnel reflection factor itself, Schlick's approximation can
be integrated analytically to derive a closed-form equation for diffuse
IAM factors for the portions of the sky and ground visible
from a tilted surface if isotropic distributions are assumed.
This function implements the integration of the
Schlick approximation provided by Xie et al. [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.
See Also
--------
pvlib.iam.schlick
Notes
-----
The analytical integration of the Schlick approximation was derived
as part of the FEDIS diffuse IAM model [2]_. Compared with the model
implemented in this function, the FEDIS model includes an additional term
to account for reflection off a pyranometer's glass dome. Because that
reflection should already be accounted for in the instrument's calibration,
the pvlib authors believe it is inappropriate to account for pyranometer
reflection again in an IAM model. Thus, this function omits that term and
implements only the integrated Schlick approximation.
Note also that the output of this function (which is an exact integration)
can be compared with the output of :py:func:`marion_diffuse` which
numerically integrates the Schlick approximation:
.. code::
>>> pvlib.iam.marion_diffuse('schlick', surface_tilt=20)
{'sky': 0.9625000227247358,
'horizon': 0.7688174948510073,
'ground': 0.6267861879241405}
>>> pvlib.iam.schlick_diffuse(surface_tilt=20)
(0.9624993421569652, 0.6269387554469255)
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`
"""
# 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 _get_model(model_name):
# check that model is implemented
model_dict = {'ashrae': ashrae, 'martin_ruiz': martin_ruiz,
'physical': physical}
try:
model = model_dict[model_name]
except KeyError:
raise NotImplementedError(f"The {model_name} model has not been "
"implemented")
return model
def _check_params(model_name, params):
# check that the parameters passed in with the model
# belong to the model
exp_params = _IAM_MODEL_PARAMS[model_name]
if set(params.keys()) != exp_params:
raise ValueError(f"The {model_name} model was expecting to be passed "
"{', '.join(list(exp_params))}, but "
"was handed {', '.join(list(params.keys()))}")
def _sin_weight(aoi):
return 1 - sind(aoi)
def _residual(aoi, source_iam, target, target_params,
weight=_sin_weight):
# computes a sum of weighted differences between the source model
# and target model, using the provided weight function
weights = weight(aoi)
# if aoi contains values outside of interval (0, 90), annihilate
# the associated weights (we don't want IAM values from AOI outside
# of (0, 90) to affect the fit; this is a possible issue when using
# `iam.fit`, but not an issue when using `iam.convert`, since in
# that case aoi is defined internally)
weights = weights * np.logical_and(aoi >= 0, aoi <= 90).astype(int)
diff = np.abs(source_iam - np.nan_to_num(target(aoi, *target_params)))
return np.sum(diff * weights)
def _get_ashrae_intercept(b):
# find x-intercept of ashrae model
return acosd(b / (1 + b))
def _ashrae_to_physical(aoi, ashrae_iam, weight, fix_n, b):
if fix_n:
# the ashrae model has an x-intercept less than 90
# we solve for this intercept, and fix n so that the physical
# model will have the same x-intercept
intercept = _get_ashrae_intercept(b)
n = sind(intercept)
# with n fixed, we will optimize for L (recall that K and L always
# appear in the physical model as a product, so it is enough to
# optimize for just L, and to fix K=4)
# we will pass n to the optimizer to simplify things later on,
# but because we are setting (n, n) as the bounds, the optimizer
# will leave n fixed
bounds = [(1e-6, 0.08), (n, n)]
guess = [0.002, n]
else:
# we don't fix n, so physical won't have same x-intercept as ashrae
# the fit will be worse, but the parameters returned for the physical
# model will be more realistic
bounds = [(1e-6, 0.08), (0.8, 2)] # L, n
guess = [0.002, 1.0]
def residual_function(target_params):
L, n = target_params
return _residual(aoi, ashrae_iam, physical, [n, 4, L], weight)
return residual_function, guess, bounds
def _martin_ruiz_to_physical(aoi, martin_ruiz_iam, weight, a_r):
# we will optimize for both n and L (recall that K and L always
# appear in the physical model as a product, so it is enough to
# optimize for just L, and to fix K=4)
# set lower bound for n at 1.0 so that x-intercept will be at 90
# order for Powell's method depends on a_r value
bounds = [(1e-6, 0.08), (1.05, 2)] # L, n
guess = [0.002, 1.1] # L, n
# get better results if we reverse order to n, L at high a_r
if a_r > 0.22:
bounds.reverse()
guess.reverse()
# the product of K and L is more important in determining an initial
# guess for the location of the minimum, so we pass L in first
def residual_function(target_params):
# unpack target_params for either search order
if target_params[0] < target_params[1]:
# L will always be less than n
L, n = target_params
else:
n, L = target_params
return _residual(aoi, martin_ruiz_iam, physical, [n, 4, L], weight)
return residual_function, guess, bounds
def _minimize(residual_function, guess, bounds, xtol):
if xtol is not None:
options = {'xtol': xtol}
else:
options = None
with np.errstate(invalid='ignore'):
optimize_result = minimize(residual_function, guess, method="powell",
bounds=bounds, options=options)
if not optimize_result.success:
try:
message = "Optimizer exited unsuccessfully:" \
+ optimize_result.message
except AttributeError:
message = "Optimizer exited unsuccessfully: \
No message explaining the failure was returned. \
If you would like to see this message, please \
update your scipy version (try version 1.8.0 \
or beyond)."
raise RuntimeError(message)
return optimize_result
def _process_return(target_name, optimize_result):
if target_name == "ashrae":
target_params = {'b': optimize_result.x.item()}
elif target_name == "martin_ruiz":
target_params = {'a_r': optimize_result.x.item()}
elif target_name == "physical":
L, n = optimize_result.x
# have to unpack order because search order may be different
if L > n:
L, n = n, L
target_params = {'n': n, 'K': 4, 'L': L}
return target_params
[docs]
def convert(source_name, source_params, target_name, weight=_sin_weight,
fix_n=True, xtol=None):
"""
Convert a source IAM model to a target IAM model.
Parameters
----------
source_name : str
Name of the source model. Must be ``'ashrae'``, ``'martin_ruiz'``, or
``'physical'``.
source_params : dict
A dictionary of parameters for the source model.
If source model is ``'ashrae'``, the dictionary must contain
the key ``'b'``.
If source model is ``'martin_ruiz'``, the dictionary must
contain the key ``'a_r'``.
If source model is ``'physical'``, the dictionary must
contain the keys ``'n'``, ``'K'``, and ``'L'``.
target_name : str
Name of the target model. Must be ``'ashrae'``, ``'martin_ruiz'``, or
``'physical'``.
weight : function, optional
A single-argument function of AOI (degrees) that calculates weights for
the residuals between models. Must return a float or an array-like
object. The default weight function is :math:`f(aoi) = 1 - sin(aoi)`.
fix_n : bool, default True
A flag to determine which method is used when converting from the
ASHRAE model to the physical model.
When ``source_name`` is ``'ashrae'`` and ``target_name`` is
``'physical'``, if `fix_n` is ``True``,
:py:func:`iam.convert` will fix ``n`` so that the returned physical
model has the same x-intercept as the inputted ASHRAE model.
Fixing ``n`` like this improves the fit of the conversion, but often
returns unrealistic values for the parameters of the physical model.
If more physically meaningful parameters are wanted,
set `fix_n` to False.
xtol : float, optional
Passed to scipy.optimize.minimize.
Returns
-------
dict
Parameters for the target model.
If target model is ``'ashrae'``, the dictionary will contain
the key ``'b'``.
If target model is ``'martin_ruiz'``, the dictionary will
contain the key ``'a_r'``.
If target model is ``'physical'``, the dictionary will
contain the keys ``'n'``, ``'K'``, and ``'L'``.
Note
----
Target model parameters are determined by minimizing
.. math::
\\sum_{\\theta=0}^{90} weight \\left(\\theta \\right) \\times
\\| source \\left(\\theta \\right) - target \\left(\\theta \\right) \\|
The sum is over :math:`\\theta = 0, 1, 2, ..., 90`.
References
----------
.. [1] Jones, A. R., Hansen, C. W., Anderson, K. S. Parameter estimation
for incidence angle modifier models for photovoltaic modules. Sandia
report SAND2023-13944 (2023).
See Also
--------
pvlib.iam.fit
pvlib.iam.ashrae
pvlib.iam.martin_ruiz
pvlib.iam.physical
"""
source = _get_model(source_name)
target = _get_model(target_name)
aoi = np.linspace(0, 90, 91)
_check_params(source_name, source_params)
source_iam = source(aoi, **source_params)
if target_name == "physical":
# we can do some special set-up to improve the fit when the
# target model is physical
if source_name == "ashrae":
residual_function, guess, bounds = \
_ashrae_to_physical(aoi, source_iam, weight, fix_n,
source_params['b'])
elif source_name == "martin_ruiz":
residual_function, guess, bounds = \
_martin_ruiz_to_physical(aoi, source_iam, weight,
source_params['a_r'])
else:
# otherwise, target model is ashrae or martin_ruiz, and scipy
# does fine without any special set-up
bounds = [(1e-04, 1)]
guess = [1e-03]
def residual_function(target_param):
return _residual(aoi, source_iam, target, target_param, weight)
optimize_result = _minimize(residual_function, guess, bounds,
xtol=xtol)
return _process_return(target_name, optimize_result)
[docs]
def fit(measured_aoi, measured_iam, model_name, weight=_sin_weight, xtol=None):
"""
Find model parameters that best fit the data.
Parameters
----------
measured_aoi : array-like
Angle of incidence values associated with the
measured IAM values. [degrees]
measured_iam : array-like
IAM values. [unitless]
model_name : str
Name of the model to be fit. Must be ``'ashrae'``, ``'martin_ruiz'``,
or ``'physical'``.
weight : function, optional
A single-argument function of AOI (degrees) that calculates weights for
the residuals between models. Must return a float or an array-like
object. The default weight function is :math:`f(aoi) = 1 - sin(aoi)`.
xtol : float, optional
Passed to scipy.optimize.minimize.
Returns
-------
dict
Parameters for target model.
If target model is ``'ashrae'``, the dictionary will contain
the key ``'b'``.
If target model is ``'martin_ruiz'``, the dictionary will
contain the key ``'a_r'``.
If target model is ``'physical'``, the dictionary will
contain the keys ``'n'``, ``'K'``, and ``'L'``.
References
----------
.. [1] Jones, A. R., Hansen, C. W., Anderson, K. S. Parameter estimation
for incidence angle modifier models for photovoltaic modules. Sandia
report SAND2023-13944 (2023).
Note
----
Model parameters are determined by minimizing
.. math::
\\sum_{AOI} weight \\left( AOI \\right) \\times
\\| IAM \\left( AOI \\right) - model \\left( AOI \\right) \\|
The sum is over ``measured_aoi`` and :math:`IAM \\left( AOI \\right)`
is ``measured_IAM``.
See Also
--------
pvlib.iam.convert
pvlib.iam.ashrae
pvlib.iam.martin_ruiz
pvlib.iam.physical
"""
target = _get_model(model_name)
if model_name == "physical":
bounds = [(0, 0.08), (1, 2)]
guess = [0.002, 1+1e-08]
def residual_function(target_params):
L, n = target_params
return _residual(measured_aoi, measured_iam, target, [n, 4, L],
weight)
# otherwise, target_name is martin_ruiz or ashrae
else:
bounds = [(1e-08, 1)]
guess = [0.05]
def residual_function(target_param):
return _residual(measured_aoi, measured_iam, target,
target_param, weight)
optimize_result = _minimize(residual_function, guess, bounds, xtol)
return _process_return(model_name, optimize_result)
