[FIX] iam.physical returns nan for aoi > 90° when n = 1 (#1706) (#1707)

* [FIX] iam.physical returns nan for aoi > 90° when n = 1 (#1706)

* address stickler-ci and codecov alerts

* suppress obnoxious numpy warnings

* whatsnew

---------

Co-authored-by: Kevin Anderson <kevin.anderso@gmail.com>

Author: kdebrab

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

@@ -58,6 +58,9 @@ Bug fixes
   :py:meth:`pvlib.modelchain.ModelChain.run_model_from_effective_irradiance`. (:issue:`1713`, :pull:`1720`)
 * ``d2mutau`` and ``NsVbi`` were hardcoded in :py:func:`pvlib.pvsystem.max_power_point` instead of
   passing through the arguments to the function. (:pull:`1733`)
+* :py:func:`pvlib.iam.physical` no longer returns NaN when ``n=1`` and ``aoi>90``.
+  This bug was introduced in v0.9.5.  (:issue:`1706`, :pull:`1707`)
+  
 
 Testing
 ~~~~~~~

pvlib/iam.py

@@ -175,8 +175,12 @@ def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None):
     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
+    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
@@ -208,13 +212,22 @@ def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None):
         tau_0 *= (1 - rho23_0) / (1 - rho23_0 * rho12_0)
 
     # transmittance after absorption in the glass
-    tau_s *= np.exp(-K * L / costheta)
-    tau_p *= np.exp(-K * L / costheta)
+    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
 
 

pvlib/tests/test_iam.py

@@ -51,6 +51,18 @@ def test_physical():
     assert_series_equal(iam, expected)
 
 
+def test_physical_n1_L0():
+    aoi = np.array([0, 22.5, 45, 67.5, 90, 100, np.nan])
+    expected = np.array([1, 1, 1, 1, 0, 0, np.nan])
+    iam = _iam.physical(aoi, n=1, L=0)
+    assert_allclose(iam, expected, equal_nan=True)
+
+    aoi = pd.Series(aoi)
+    expected = pd.Series(expected)
+    iam = _iam.physical(aoi, n=1, L=0)
+    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])