Incidence Angle Modifier zero instead of np.nan (#338) (#339)

* changed np.nan to 0 in functions ashraeiam, physicaliam, sapm_aoi_loss

* updated docstrings

* updated tests to account for changes in functions

* replace 0-degree aoi in physicaliam with 1e-06, analogous to the calculation of tau0

* updated tests to account for input angles of np.nan

* updated test for ashraeiam to check for possible negative values for input angles close to 90°

* updated docstring with info about nan-inputs resulting in nan.

* updated variables-and-symbols-page with allowable range of aoi

* changed whatsnew for v0.5

* applied @adriesse suggestion to split calculation of tau

* forgot to delete unused line

* typo :-(

* typo in variables_styles_rules

* implemented code suggestion by @adriesse

Author: jkfm

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

@@ -0,0 +1,29 @@
+.. _whatsnew_0500:
+
+v0.5.0 ()
+---------
+
+
+Bug fixes
+~~~~~~~~~
+
+
+Enhancements
+~~~~~~~~~~~~
+
+
+API Changes
+~~~~~~~~~~~
+
+* Changes calculation of the Incidence Angle Modifier to return 0 instead of np.nan for angles >= 90°.
+  This improves the calculation of effective irradiance close to sunrise and sunset. (:issue:`338`)
+
+
+Documentation
+~~~~~~~~~~~~~
+
+
+Contributors
+~~~~~~~~~~~~
+
+* Johannes Kaufmann

pvlib/data/variables_style_rules.csv

@@ -6,7 +6,7 @@ dni;direct normal irradiance
 dni_extra;direct normal irradiance at top of atmosphere (extraterrestrial)
 dhi;diffuse horizontal irradiance
 ghi;global horizontal irradiance
-aoi;angle of incidence
+aoi;angle of incidence between :math:`90\deg` and :math:`90\deg`
 aoi_projection;cos(aoi)
 airmass;airmass
 airmass_relative;relative airmass

pvlib/pvsystem.py

@@ -700,7 +700,7 @@ def ashraeiam(aoi, b=0.05):
     ----------
     aoi : numeric
         The angle of incidence between the module normal vector and the
-        sun-beam vector in degrees.
+        sun-beam vector in degrees. Angles of nan will result in nan.
 
     b : float, default 0.05
         A parameter to adjust the modifier as a function of angle of
@@ -712,7 +712,7 @@ def ashraeiam(aoi, b=0.05):
         The incident angle modifier calculated as 1-b*(sec(aoi)-1) as
         described in [2,3].
 
-        Returns nan for all abs(aoi) >= 90 and for all IAM values that
+        Returns zeros for all abs(aoi) >= 90 and for all IAM values that
         would be less than 0.
 
     References
@@ -735,7 +735,7 @@ def ashraeiam(aoi, b=0.05):
 
     iam = 1 - b*((1/np.cos(np.radians(aoi)) - 1))
 
-    iam = np.where(np.abs(aoi) >= 90, np.nan, iam)
+    iam = np.where(np.abs(aoi) >= 90, 0, iam)
     iam = np.maximum(0, iam)
 
     if isinstance(iam, pd.Series):
@@ -764,7 +764,8 @@ def physicaliam(aoi, n=1.526, K=4., L=0.002):
     ----------
     aoi : numeric
         The angle of incidence between the module normal vector and the
-        sun-beam vector in degrees.
+        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.
 
     n : numeric, default 1.526
         The effective index of refraction (unitless). Reference [1]
@@ -814,27 +815,40 @@ def physicaliam(aoi, n=1.526, K=4., L=0.002):
     spa
     ashraeiam
     '''
+    zeroang = 1e-06
+
+    aoi = np.where(aoi == 0, zeroang, aoi)
+
+    # angle of reflection
     thetar_deg = tools.asind(1.0 / n*(tools.sind(aoi)))
 
-    tau = (np.exp(- 1.0 * (K*L / tools.cosd(thetar_deg))) *
-           ((1 - 0.5*((((tools.sind(thetar_deg - aoi)) ** 2) /
-            ((tools.sind(thetar_deg + aoi)) ** 2) +
-            ((tools.tand(thetar_deg - aoi)) ** 2) /
-            ((tools.tand(thetar_deg + aoi)) ** 2))))))
+    # reflectance and transmittance for normal incidence light
+    rho_zero = ((1-n) / (1+n)) ** 2
+    tau_zero = np.exp(-K*L)
 
-    zeroang = 1e-06
+    # reflectance for parallel and perpendicular polarized light
+    rho_para = (tools.tand(thetar_deg - aoi) /
+                tools.tand(thetar_deg + aoi)) ** 2
+    rho_perp = (tools.sind(thetar_deg - aoi) /
+                tools.sind(thetar_deg + aoi)) ** 2
+
+    # transmittance for non-normal light
+    tau = np.exp(-K*L / tools.cosd(thetar_deg))
 
-    thetar_deg0 = tools.asind(1.0 / n*(tools.sind(zeroang)))
+    # 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)
 
-    tau0 = (np.exp(- 1.0 * (K*L / tools.cosd(thetar_deg0))) *
-            ((1 - 0.5*((((tools.sind(thetar_deg0 - zeroang)) ** 2) /
-             ((tools.sind(thetar_deg0 + zeroang)) ** 2) +
-             ((tools.tand(thetar_deg0 - zeroang)) ** 2) /
-             ((tools.tand(thetar_deg0 + zeroang)) ** 2))))))
+    # 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)
 
-    iam = tau / tau0
+    # 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)
 
-    iam = np.where((np.abs(aoi) >= 90) | (iam < 0), np.nan, iam)
+    # for light coming from behind the plane, none can enter the module
+    iam = np.where(aoi > 90, 0, iam)
 
     if isinstance(aoi, pd.Series):
         iam = pd.Series(iam, index=aoi.index)
@@ -1465,7 +1479,7 @@ def sapm_aoi_loss(aoi, module, upper=None):
     ----------
     aoi : numeric
         Angle of incidence in degrees. Negative input angles will return
-        nan values.
+        zeros.
 
     module : dict-like
         A dict, Series, or DataFrame defining the SAPM performance
@@ -1507,7 +1521,7 @@ def sapm_aoi_loss(aoi, module, upper=None):
 
     aoi_loss = np.polyval(aoi_coeff, aoi)
     aoi_loss = np.clip(aoi_loss, 0, upper)
-    aoi_loss = np.where(aoi < 0, np.nan, aoi_loss)
+    aoi_loss = np.where(aoi < 0, 0, aoi_loss)
 
     if isinstance(aoi, pd.Series):
         aoi_loss = pd.Series(aoi_loss, aoi.index)

pvlib/test/test_pvsystem.py

@@ -94,41 +94,41 @@ def test_systemdef_dict():
 
 @needs_numpy_1_10
 def test_ashraeiam():
-    thetas = np.linspace(-90, 90, 9)
+    thetas = np.array([-90. , -67.5, -45. , -22.5,   0. ,  22.5,  45. ,  67.5, 89.,  90. , np.nan])
     iam = pvsystem.ashraeiam(thetas, .05)
-    expected = np.array([        nan,  0.9193437 ,  0.97928932,  0.99588039,  1.        ,
-        0.99588039,  0.97928932,  0.9193437 ,         nan])
+    expected = np.array([        0,  0.9193437 ,  0.97928932,  0.99588039,  1.        ,
+        0.99588039,  0.97928932,  0.9193437 ,         0, 0,  np.nan])
     assert_allclose(iam, expected, equal_nan=True)
 
 
 @needs_numpy_1_10
 def test_PVSystem_ashraeiam():
     module_parameters = pd.Series({'b': 0.05})
     system = pvsystem.PVSystem(module_parameters=module_parameters)
-    thetas = np.linspace(-90, 90, 9)
+    thetas = np.array([-90. , -67.5, -45. , -22.5,   0. ,  22.5,  45. ,  67.5,  89., 90. , np.nan])
     iam = system.ashraeiam(thetas)
-    expected = np.array([        nan,  0.9193437 ,  0.97928932,  0.99588039,  1.        ,
-        0.99588039,  0.97928932,  0.9193437 ,         nan])
+    expected = np.array([        0,  0.9193437 ,  0.97928932,  0.99588039,  1.        ,
+        0.99588039,  0.97928932,  0.9193437 ,         0, 0,  np.nan])
     assert_allclose(iam, expected, equal_nan=True)
 
 
 @needs_numpy_1_10
 def test_physicaliam():
-    thetas = np.linspace(-90, 90, 9)
+    thetas = np.array([-90. , -67.5, -45. , -22.5,   0. ,  22.5,  45. ,  67.5,  90. , np.nan])
     iam = pvsystem.physicaliam(thetas, 1.526, 0.002, 4)
-    expected = np.array([        nan,  0.8893998 ,  0.98797788,  0.99926198,         nan,
-        0.99926198,  0.98797788,  0.8893998 ,         nan])
+    expected = np.array([        0,  0.8893998,  0.98797788,  0.99926198,         1,
+        0.99926198,  0.98797788,  0.8893998,         0, np.nan])
     assert_allclose(iam, expected, equal_nan=True)
 
 
 @needs_numpy_1_10
 def test_PVSystem_physicaliam():
     module_parameters = pd.Series({'K': 4, 'L': 0.002, 'n': 1.526})
     system = pvsystem.PVSystem(module_parameters=module_parameters)
-    thetas = np.linspace(-90, 90, 9)
+    thetas = np.array([-90. , -67.5, -45. , -22.5,   0. ,  22.5,  45. ,  67.5,  90. , np.nan])
     iam = system.physicaliam(thetas)
-    expected = np.array([        nan,  0.8893998 ,  0.98797788,  0.99926198,         nan,
-        0.99926198,  0.98797788,  0.8893998 ,         nan])
+    expected = np.array([        0,  0.8893998 ,  0.98797788,  0.99926198,         1,
+        0.99926198,  0.98797788,  0.8893998 ,         0, np.nan])
     assert_allclose(iam, expected, equal_nan=True)
 
 
@@ -239,7 +239,7 @@ def test_PVSystem_sapm_spectral_loss(sapm_module_params):
 @pytest.mark.parametrize('aoi,expected', [
     (45, 0.9975036250000002),
     (np.array([[-30, 30, 100, np.nan]]),
-     np.array([[np.nan, 1.007572, 0, np.nan]])),
+     np.array([[0, 1.007572, 0, np.nan]])),
     (pd.Series([80]), pd.Series([0.597472]))
 ])
 def test_sapm_aoi_loss(sapm_module_params, aoi, expected):