From ae3b8ffb073316cfad1e416d2865b11d84998737 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <mark.campanelli@gmail.com>
Date: Sun, 5 Jul 2026 11:59:08 -0600
Subject: [PATCH 01/14] Handle 0/0 in perez without NaNs

---
 pvlib/irradiance.py      | 17 +++++++++++++----
 tests/test_irradiance.py | 16 ++++++++++++++--
 2 files changed, 27 insertions(+), 6 deletions(-)

diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 9db6da3aa8..4badb5b9fb 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -1147,12 +1147,21 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     kappa = 1.041  # for solar_zenith in radians
     z = np.radians(solar_zenith)  # convert to radians
 
-    # delta is the sky's "brightness"
+    # delta is the sky's "brightness", assumes dni_extra > 0
     delta = dhi * airmass / dni_extra
 
-    # epsilon is the sky's "clearness"
-    with np.errstate(invalid='ignore'):
-        eps = ((dhi + dni) / dhi + kappa * (z ** 3)) / (1 + kappa * (z ** 3))
+    # epsilon is the sky's "clearness". Preserves NaNs for dni or dhi.
+    # Assumes:
+    # - dni >=0 and dhi >= 0.
+    # - dni->0^+ faster than dhi->0^+.
+    irr_ratio = np.zeros(np.broadcast(dni, dhi).shape)
+    with np.errstate(divide="raise"):
+        irr_ratio = np.divide(
+            dni, dhi, out=irr_ratio, where=np.logical_not(
+                np.logical_and(dhi == 0, dni == 0)
+            )
+        )
+    eps = 1 + irr_ratio / (1 + kappa * (z ** 3))
 
     # numpy indexing below will not work with a Series
     if isinstance(eps, pd.Series):
diff --git a/tests/test_irradiance.py b/tests/test_irradiance.py
index acc8495d3b..84b1828e7a 100644
--- a/tests/test_irradiance.py
+++ b/tests/test_irradiance.py
@@ -7,8 +7,9 @@
 import pandas as pd
 
 import pytest
-from numpy.testing import (assert_almost_equal,
-                           assert_allclose)
+from numpy.testing import (
+    assert_almost_equal, assert_allclose, assert_equal, assert_raises
+)
 from pvlib import irradiance, albedo
 
 from .conftest import (
@@ -393,6 +394,17 @@ def test_perez_negative_horizon():
     assert_series_equal(sum_components, expected_for_sum, check_less_precise=2)
 
 
+def test_perez_zero_dhi_and_dni(dni_et):
+    out = irradiance.perez(20, 180, 0.0, 0.0, dni_et, 89.96, 256.28, 37.32)
+    expected = 0.0
+    assert_equal(out, expected)
+
+
+def test_perez_zero_dhi_nonzero_dni(dni_et):
+    with assert_raises(FloatingPointError):
+        irradiance.perez(20, 180, 0.0, 100.0, dni_et, 89.96, 256.28, 37.32)
+
+
 def test_perez_arrays(irrad_data, ephem_data, dni_et, relative_airmass):
     dni = irrad_data['dni'].copy()
     dni.iloc[2] = np.nan

From 8755ecccefc6d52aadedd9f9b6a1bc4e3c13ef15 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <mark.campanelli@gmail.com>
Date: Sun, 5 Jul 2026 13:27:25 -0600
Subject: [PATCH 02/14] Add return_components=True path to tests

---
 tests/test_irradiance.py | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/tests/test_irradiance.py b/tests/test_irradiance.py
index 84b1828e7a..2a9cc9d388 100644
--- a/tests/test_irradiance.py
+++ b/tests/test_irradiance.py
@@ -399,11 +399,29 @@ def test_perez_zero_dhi_and_dni(dni_et):
     expected = 0.0
     assert_equal(out, expected)
 
+    out = irradiance.perez(
+        20, 180, 0.0, 0.0, dni_et, 89.96, 256.28, 37.32, return_components=True
+    )
+    expected = {
+        "poa_sky_diffuse": 0.0,
+        "poa_isotropic": 0.0,
+        "poa_circumsolar": 0.0,
+        "poa_horizon": 0.0,
+    }
+    assert len(out) == len(expected)
+    for key in expected.keys():
+        assert_equal(out[key], expected[key])
+
 
 def test_perez_zero_dhi_nonzero_dni(dni_et):
     with assert_raises(FloatingPointError):
         irradiance.perez(20, 180, 0.0, 100.0, dni_et, 89.96, 256.28, 37.32)
 
+    with assert_raises(FloatingPointError):
+        irradiance.perez(
+            20, 180, 0.0, 100.0, dni_et, 89.96, 256.28, 37.32, return_components=True
+        )
+
 
 def test_perez_arrays(irrad_data, ephem_data, dni_et, relative_airmass):
     dni = irrad_data['dni'].copy()

From 60bf26d40d4f175943dcc826d5855b75fa68fd53 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <mark.campanelli@gmail.com>
Date: Sun, 5 Jul 2026 13:28:46 -0600
Subject: [PATCH 03/14] Appease flake8

---
 tests/test_irradiance.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/tests/test_irradiance.py b/tests/test_irradiance.py
index 2a9cc9d388..8a75e3de26 100644
--- a/tests/test_irradiance.py
+++ b/tests/test_irradiance.py
@@ -400,7 +400,8 @@ def test_perez_zero_dhi_and_dni(dni_et):
     assert_equal(out, expected)
 
     out = irradiance.perez(
-        20, 180, 0.0, 0.0, dni_et, 89.96, 256.28, 37.32, return_components=True
+        20, 180, 0.0, 0.0, dni_et, 89.96, 256.28, 37.32,
+        return_components=True
     )
     expected = {
         "poa_sky_diffuse": 0.0,
@@ -419,7 +420,8 @@ def test_perez_zero_dhi_nonzero_dni(dni_et):
 
     with assert_raises(FloatingPointError):
         irradiance.perez(
-            20, 180, 0.0, 100.0, dni_et, 89.96, 256.28, 37.32, return_components=True
+            20, 180, 0.0, 100.0, dni_et, 89.96, 256.28, 37.32,
+            return_components=True
         )
 
 

From 483c9b8f504efd0bde51881d61aa33c17af3a340 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <mark.campanelli@gmail.com>
Date: Sun, 5 Jul 2026 14:05:10 -0600
Subject: [PATCH 04/14] Add and improve tests

---
 tests/test_irradiance.py | 74 ++++++++++++++++++++++++++++++++--------
 1 file changed, 60 insertions(+), 14 deletions(-)

diff --git a/tests/test_irradiance.py b/tests/test_irradiance.py
index 8a75e3de26..38c2b73a41 100644
--- a/tests/test_irradiance.py
+++ b/tests/test_irradiance.py
@@ -288,7 +288,6 @@ def test_perez_driesse_airmass(irrad_data, ephem_data, dni_et):
     out = irradiance.perez_driesse(40, 180, irrad_data['dhi'], dni,
                                    dni_et, ephem_data['apparent_zenith'],
                                    ephem_data['azimuth'], airmass=None)
-    print(out)
     expected = pd.Series(np.array(
         [0.,   29.991,  np.nan,   47.397]),
         index=irrad_data.index)
@@ -394,15 +393,15 @@ def test_perez_negative_horizon():
     assert_series_equal(sum_components, expected_for_sum, check_less_precise=2)
 
 
-def test_perez_zero_dhi_and_dni(dni_et):
-    out = irradiance.perez(20, 180, 0.0, 0.0, dni_et, 89.96, 256.28, 37.32)
+def test_perez_zero_dhi_and_dni_scalar():
+    # Divides zero by zero.
+    args = (20, 180, 0.0, 0.0, 1366.1, 89.96, 256.28, 37.32)
+
+    out = irradiance.perez(*args)
     expected = 0.0
     assert_equal(out, expected)
 
-    out = irradiance.perez(
-        20, 180, 0.0, 0.0, dni_et, 89.96, 256.28, 37.32,
-        return_components=True
-    )
+    out = irradiance.perez(*args, return_components=True)
     expected = {
         "poa_sky_diffuse": 0.0,
         "poa_isotropic": 0.0,
@@ -411,18 +410,65 @@ def test_perez_zero_dhi_and_dni(dni_et):
     }
     assert len(out) == len(expected)
     for key in expected.keys():
-        assert_equal(out[key], expected[key])
+        assert_equal(out[key], expected[key], err_msg=key)
+
+    assert_equal(
+        out["poa_sky_diffuse"],
+        out["poa_isotropic"] + out["poa_circumsolar"] + out["poa_horizon"],
+    )
+
+
+def test_perez_zero_dhi_and_dni_array():
+    # Divides zero by zero.
+    args = (
+        20, 180, np.array([0.0, 10.0, float("nan")]), 0.0, 1366.1, 89.96,
+        256.28, 37.32
+    )
+
+    out = irradiance.perez(*args)
+    expected = np.array([0.0, 9.424924186619206, float("nan")])
+    assert_allclose(out, expected)
 
+    out = irradiance.perez(*args, return_components=True)
+    expected = {
+        "poa_sky_diffuse": np.array([0.0, 9.424924186619206, float("nan")]),
+        "poa_isotropic": np.array([ 0.0, 9.162258932459126, float("nan")]),
+        "poa_circumsolar": np.array([ 0.0, 0.5187450944545264, float("nan")]),
+        "poa_horizon": np.array([0.0, -0.2560798402944465, float("nan")]),
+    }
+    assert len(out) == len(expected)
+    for key in expected.keys():
+        assert_allclose(out[key], expected[key], err_msg=key)
+
+    assert_almost_equal(
+        out["poa_sky_diffuse"],
+        out["poa_isotropic"] + out["poa_circumsolar"] + out["poa_horizon"],
+    )
+
+
+def test_perez_zero_dhi_nonzero_dni_scalar():
+    # Divides nonzero by zero.
+    args = (20, 180, 0.0, 100.0, 1366.1, 89.96, 256.28, 37.32)
+
+    with assert_raises(FloatingPointError):
+        irradiance.perez(*args)
+
+    with assert_raises(FloatingPointError):
+        irradiance.perez(*args, return_components=True)
+
+
+def test_perez_zero_dhi_nonzero_dni_array():
+    # Divides nonzero by zero.
+    args = (
+        20, 180, np.array([0.0, 10.0, float("nan")]), 100.0, 1366.1, 89.96,
+        256.28, 37.32
+    )
 
-def test_perez_zero_dhi_nonzero_dni(dni_et):
     with assert_raises(FloatingPointError):
-        irradiance.perez(20, 180, 0.0, 100.0, dni_et, 89.96, 256.28, 37.32)
+        irradiance.perez(*args)
 
     with assert_raises(FloatingPointError):
-        irradiance.perez(
-            20, 180, 0.0, 100.0, dni_et, 89.96, 256.28, 37.32,
-            return_components=True
-        )
+        irradiance.perez(*args, return_components=True)
 
 
 def test_perez_arrays(irrad_data, ephem_data, dni_et, relative_airmass):

From e97b83f9d45d3c55bdaac0484349deba1dc19a12 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <mark.campanelli@gmail.com>
Date: Sun, 5 Jul 2026 14:19:28 -0600
Subject: [PATCH 05/14] Make sure NaN irradiances pass through as expected

---
 pvlib/irradiance.py      |  2 +-
 tests/test_irradiance.py | 27 ++++++++++++++++++---------
 2 files changed, 19 insertions(+), 10 deletions(-)

diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 4badb5b9fb..358fa8561f 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -1147,7 +1147,7 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     kappa = 1.041  # for solar_zenith in radians
     z = np.radians(solar_zenith)  # convert to radians
 
-    # delta is the sky's "brightness", assumes dni_extra > 0
+    # delta is the sky's "brightness", NaN airmass ok, assumes dni_extra > 0
     delta = dhi * airmass / dni_extra
 
     # epsilon is the sky's "clearness". Preserves NaNs for dni or dhi.
diff --git a/tests/test_irradiance.py b/tests/test_irradiance.py
index 38c2b73a41..1551fb352a 100644
--- a/tests/test_irradiance.py
+++ b/tests/test_irradiance.py
@@ -418,23 +418,32 @@ def test_perez_zero_dhi_and_dni_scalar():
     )
 
 
-def test_perez_zero_dhi_and_dni_array():
+def test_perez_array_dhi_and_dni_combos():
     # Divides zero by zero.
     args = (
-        20, 180, np.array([0.0, 10.0, float("nan")]), 0.0, 1366.1, 89.96,
-        256.28, 37.32
+        20, 180, np.array([0.0, 10.0, np.nan, 0.0, np.nan]),
+        np.array([0.0, 0.0, 0.0, np.nan, np.nan]), 1366.1, 89.96, 256.28,
+        37.32
     )
 
     out = irradiance.perez(*args)
-    expected = np.array([0.0, 9.424924186619206, float("nan")])
+    expected = np.array([0.0, 9.424924186619206, np.nan, np.nan, np.nan])
     assert_allclose(out, expected)
 
     out = irradiance.perez(*args, return_components=True)
     expected = {
-        "poa_sky_diffuse": np.array([0.0, 9.424924186619206, float("nan")]),
-        "poa_isotropic": np.array([ 0.0, 9.162258932459126, float("nan")]),
-        "poa_circumsolar": np.array([ 0.0, 0.5187450944545264, float("nan")]),
-        "poa_horizon": np.array([0.0, -0.2560798402944465, float("nan")]),
+        "poa_sky_diffuse": np.array(
+            [0.0, 9.424924186619206, np.nan, np.nan, np.nan]
+        ),
+        "poa_isotropic": np.array(
+            [ 0.0, 9.162258932459126, np.nan, np.nan, np.nan]
+        ),
+        "poa_circumsolar": np.array(
+            [ 0.0, 0.5187450944545264, np.nan, np.nan, np.nan]
+        ),
+        "poa_horizon": np.array(
+            [0.0, -0.2560798402944465, np.nan, np.nan, np.nan]
+        ),
     }
     assert len(out) == len(expected)
     for key in expected.keys():
@@ -460,7 +469,7 @@ def test_perez_zero_dhi_nonzero_dni_scalar():
 def test_perez_zero_dhi_nonzero_dni_array():
     # Divides nonzero by zero.
     args = (
-        20, 180, np.array([0.0, 10.0, float("nan")]), 100.0, 1366.1, 89.96,
+        20, 180, np.array([0.0, 10.0, np.nan]), 100.0, 1366.1, 89.96,
         256.28, 37.32
     )
 

From 303ae0101c2e4c6d2ec6a7f33a458f4adae0e0f0 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <mark.campanelli@gmail.com>
Date: Sun, 5 Jul 2026 14:33:09 -0600
Subject: [PATCH 06/14] Improve NaN DNI and DHI division combo testing

---
 tests/test_irradiance.py | 31 ++++++++++++++++---------------
 1 file changed, 16 insertions(+), 15 deletions(-)

diff --git a/tests/test_irradiance.py b/tests/test_irradiance.py
index 1551fb352a..d7fe6567b9 100644
--- a/tests/test_irradiance.py
+++ b/tests/test_irradiance.py
@@ -398,12 +398,12 @@ def test_perez_zero_dhi_and_dni_scalar():
     args = (20, 180, 0.0, 0.0, 1366.1, 89.96, 256.28, 37.32)
 
     out = irradiance.perez(*args)
-    expected = 0.0
-    assert_equal(out, expected)
+    poa_sky_diffuse_expected = 0.0
+    assert_equal(out, poa_sky_diffuse_expected)
 
     out = irradiance.perez(*args, return_components=True)
     expected = {
-        "poa_sky_diffuse": 0.0,
+        "poa_sky_diffuse": poa_sky_diffuse_expected,
         "poa_isotropic": 0.0,
         "poa_circumsolar": 0.0,
         "poa_horizon": 0.0,
@@ -419,30 +419,31 @@ def test_perez_zero_dhi_and_dni_scalar():
 
 
 def test_perez_array_dhi_and_dni_combos():
-    # Divides zero by zero.
+    # Divides zero and non-zero by zero and various NaN division combos.
     args = (
-        20, 180, np.array([0.0, 10.0, np.nan, 0.0, np.nan]),
-        np.array([0.0, 0.0, 0.0, np.nan, np.nan]), 1366.1, 89.96, 256.28,
-        37.32
+        20, 180,
+        np.array([0.0, 10.0, np.nan, np.nan, 0.0, 100.0, np.nan]),
+        np.array([0.0, 0.0, 0.0, 100.0, np.nan, np.nan, np.nan]),
+        1366.1, 89.96, 256.28, 37.32
     )
 
     out = irradiance.perez(*args)
-    expected = np.array([0.0, 9.424924186619206, np.nan, np.nan, np.nan])
-    assert_allclose(out, expected)
+    poa_sky_diffuse_expected = np.array(
+        [0.0, 9.424924186619206, np.nan, np.nan, np.nan, np.nan, np.nan]
+    )
+    assert_allclose(out, poa_sky_diffuse_expected)
 
     out = irradiance.perez(*args, return_components=True)
     expected = {
-        "poa_sky_diffuse": np.array(
-            [0.0, 9.424924186619206, np.nan, np.nan, np.nan]
-        ),
+        "poa_sky_diffuse": poa_sky_diffuse_expected,
         "poa_isotropic": np.array(
-            [ 0.0, 9.162258932459126, np.nan, np.nan, np.nan]
+            [0.0, 9.162258932459126, np.nan, np.nan, np.nan, np.nan, np.nan]
         ),
         "poa_circumsolar": np.array(
-            [ 0.0, 0.5187450944545264, np.nan, np.nan, np.nan]
+            [0.0, 0.5187450944545264, np.nan, np.nan, np.nan, np.nan, np.nan]
         ),
         "poa_horizon": np.array(
-            [0.0, -0.2560798402944465, np.nan, np.nan, np.nan]
+            [0.0, -0.2560798402944465, np.nan, np.nan, np.nan, np.nan, np.nan]
         ),
     }
     assert len(out) == len(expected)

From 46970d219b04766e5666b0e249c064e38ae1e2ed Mon Sep 17 00:00:00 2001
From: Mark Campanelli <mark.campanelli@gmail.com>
Date: Sun, 5 Jul 2026 14:51:05 -0600
Subject: [PATCH 07/14] Push breaking test for discussion purposes

---
 pvlib/irradiance.py      | 17 +++++++++++++----
 tests/test_irradiance.py | 39 +++++++++++++++++++++++++++++++++++++++
 2 files changed, 52 insertions(+), 4 deletions(-)

diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 358fa8561f..b863eb5aa0 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -1042,10 +1042,14 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     Perez models determine the diffuse irradiance from the sky (ground
     reflected irradiance is not included in this algorithm) on a tilted
     surface using the surface tilt angle, surface azimuth angle, diffuse
-    horizontal irradiance, direct normal irradiance, extraterrestrial
-    irradiance, sun zenith angle, sun azimuth angle, and relative (not
-    pressure-corrected) airmass. Optionally a selector may be used to
-    use any of Perez's model coefficient sets.
+    horizontal irradiance (DHI), direct normal irradiance (DNI),
+    extraterrestrial irradiance, sun zenith angle, sun azimuth angle, and
+    relative (not pressure-corrected) airmass. Optionally a selector may be
+    used to use any of Perez's model coefficient sets. It is expected that if
+    DHI is zero, then DNI is also zero, otherwise a FloatingPointError is
+    raised due to a division of a nonzero (and not NaN) value by zero. It is
+    also expected that extraterrestrial irradiance is positive. If airmass
+    is NaN, then the total and all components are zero.
 
     Warning
     -------
@@ -1125,6 +1129,11 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
             * poa_circumsolar
             * poa_horizon
 
+    Raises
+    ------
+    FloatingPointError
+        If dni is zero when dhi is not zero and not NaN.
+
 
     References
     ----------
diff --git a/tests/test_irradiance.py b/tests/test_irradiance.py
index d7fe6567b9..7d15e5cd3d 100644
--- a/tests/test_irradiance.py
+++ b/tests/test_irradiance.py
@@ -456,6 +456,45 @@ def test_perez_array_dhi_and_dni_combos():
     )
 
 
+def test_perez_array_dhi_and_dni_combos_nan_airmass():
+    # Divides zero and non-zero by zero and various NaN division combos, when
+    # airmass is NaN (e.g., for sun below horizon).
+    args = (
+        20, 180,
+        np.array([0.0, 10.0, np.nan, np.nan, 0.0, 100.0, np.nan]),
+        np.array([0.0, 0.0, 0.0, 100.0, np.nan, np.nan, np.nan]),
+        1366.1, 89.96, 256.28, np.nan
+    )
+
+    out = irradiance.perez(*args)
+    poa_sky_diffuse_expected = np.array(
+        [0.0, 0.0, np.nan, np.nan, np.nan, np.nan, np.nan]
+    )
+    assert_allclose(out, poa_sky_diffuse_expected)
+
+    out = irradiance.perez(*args, return_components=True)
+    expected = {
+        "poa_sky_diffuse": poa_sky_diffuse_expected,
+        "poa_isotropic": np.array(
+            [0.0, 9.162258932459126, np.nan, np.nan, np.nan, np.nan, np.nan]
+        ),
+        "poa_circumsolar": np.array(
+            [0.0, 0.5187450944545264, np.nan, np.nan, np.nan, np.nan, np.nan]
+        ),
+        "poa_horizon": np.array(
+            [0.0, -0.2560798402944465, np.nan, np.nan, np.nan, np.nan, np.nan]
+        ),
+    }
+    assert len(out) == len(expected)
+    for key in expected.keys():
+        assert_allclose(out[key], expected[key], err_msg=key)
+
+    assert_almost_equal(
+        out["poa_sky_diffuse"],
+        out["poa_isotropic"] + out["poa_circumsolar"] + out["poa_horizon"],
+    )
+
+
 def test_perez_zero_dhi_nonzero_dni_scalar():
     # Divides nonzero by zero.
     args = (20, 180, 0.0, 100.0, 1366.1, 89.96, 256.28, 37.32)

From 6f6d44caaad42db2d48495451652018ad1314a86 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <mark.campanelli@gmail.com>
Date: Sun, 5 Jul 2026 15:11:08 -0600
Subject: [PATCH 08/14] Fix remaining expected values in broken test

---
 tests/test_irradiance.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/test_irradiance.py b/tests/test_irradiance.py
index 7d15e5cd3d..ecee03a805 100644
--- a/tests/test_irradiance.py
+++ b/tests/test_irradiance.py
@@ -476,13 +476,13 @@ def test_perez_array_dhi_and_dni_combos_nan_airmass():
     expected = {
         "poa_sky_diffuse": poa_sky_diffuse_expected,
         "poa_isotropic": np.array(
-            [0.0, 9.162258932459126, np.nan, np.nan, np.nan, np.nan, np.nan]
+            [0.0, 0.0, np.nan, np.nan, np.nan, np.nan, np.nan]
         ),
         "poa_circumsolar": np.array(
-            [0.0, 0.5187450944545264, np.nan, np.nan, np.nan, np.nan, np.nan]
+            [0.0, 0.0, np.nan, np.nan, np.nan, np.nan, np.nan]
         ),
         "poa_horizon": np.array(
-            [0.0, -0.2560798402944465, np.nan, np.nan, np.nan, np.nan, np.nan]
+            [0.0, 0.0, np.nan, np.nan, np.nan, np.nan, np.nan]
         ),
     }
     assert len(out) == len(expected)

From 4f1b6641b6e0d8112d049e17dd276911dd12f3fe Mon Sep 17 00:00:00 2001
From: Mark Campanelli <mark.campanelli@gmail.com>
Date: Sun, 5 Jul 2026 15:15:03 -0600
Subject: [PATCH 09/14] Update docstring for expected NaN airmass behavior

---
 pvlib/irradiance.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index b863eb5aa0..1c403cb2f9 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -1049,7 +1049,8 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     DHI is zero, then DNI is also zero, otherwise a FloatingPointError is
     raised due to a division of a nonzero (and not NaN) value by zero. It is
     also expected that extraterrestrial irradiance is positive. If airmass
-    is NaN, then the total and all components are zero.
+    is NaN, then the total and all components are zero if they should not 
+    otherwise be NaN.
 
     Warning
     -------

From f6c80709c55c24c6e8673cbc6024cad4e7049fcb Mon Sep 17 00:00:00 2001
From: Mark Campanelli <mark.campanelli@gmail.com>
Date: Sun, 5 Jul 2026 16:39:50 -0600
Subject: [PATCH 10/14] Do not overwrite approapriate NaNs with zero when
 airmass is NaN

---
 pvlib/irradiance.py      | 13 +++++++++----
 tests/test_irradiance.py |  2 +-
 2 files changed, 10 insertions(+), 5 deletions(-)

diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 1c403cb2f9..20df3966ba 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -1049,8 +1049,8 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     DHI is zero, then DNI is also zero, otherwise a FloatingPointError is
     raised due to a division of a nonzero (and not NaN) value by zero. It is
     also expected that extraterrestrial irradiance is positive. If airmass
-    is NaN, then the total and all components are zero if they should not 
-    otherwise be NaN.
+    is NaN, then the total and all components are zero if they should not
+    otherwise be NaN because DHI or DNI was NaN.
 
     Warning
     -------
@@ -1218,10 +1218,15 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     sky_diffuse = np.maximum(dhi * (term1 + term2 + term3), 0)
 
     # we've preserved the input type until now, so don't ruin it!
+    airmass_nan_idx = np.logical_and(
+        np.isnan(airmass), np.logical_not(
+            np.logical_or(np.isnan(dhi), np.isnan(dni))
+        )
+    )
     if isinstance(sky_diffuse, pd.Series):
-        sky_diffuse[np.isnan(airmass)] = 0
+        sky_diffuse[airmass_nan_idx] = 0
     else:
-        sky_diffuse = np.where(np.isnan(airmass), 0, sky_diffuse)
+        sky_diffuse = np.where(airmass_nan_idx, 0, sky_diffuse)
 
     if return_components:
         diffuse_components = OrderedDict()
diff --git a/tests/test_irradiance.py b/tests/test_irradiance.py
index ecee03a805..94057b0ed7 100644
--- a/tests/test_irradiance.py
+++ b/tests/test_irradiance.py
@@ -463,7 +463,7 @@ def test_perez_array_dhi_and_dni_combos_nan_airmass():
         20, 180,
         np.array([0.0, 10.0, np.nan, np.nan, 0.0, 100.0, np.nan]),
         np.array([0.0, 0.0, 0.0, 100.0, np.nan, np.nan, np.nan]),
-        1366.1, 89.96, 256.28, np.nan
+        1366.1, 91, 256.28, np.nan
     )
 
     out = irradiance.perez(*args)

From 2a2185f27ad253b5cf460db34a962285b9c4150e Mon Sep 17 00:00:00 2001
From: Mark Campanelli <mark.campanelli@gmail.com>
Date: Sun, 5 Jul 2026 17:23:06 -0600
Subject: [PATCH 11/14] Update some comments and whatsnew

---
 docs/sphinx/source/whatsnew/v0.15.3.rst | 6 +++++-
 pvlib/irradiance.py                     | 9 +++++----
 2 files changed, 10 insertions(+), 5 deletions(-)

diff --git a/docs/sphinx/source/whatsnew/v0.15.3.rst b/docs/sphinx/source/whatsnew/v0.15.3.rst
index eb36b659ca..0923bf1cd3 100644
--- a/docs/sphinx/source/whatsnew/v0.15.3.rst
+++ b/docs/sphinx/source/whatsnew/v0.15.3.rst
@@ -14,7 +14,10 @@ Deprecations
 
 Bug fixes
 ~~~~~~~~~
-
+* Fix :py:func:`pvlib.irradiance.perez` to not return NaN when both DHI and
+  DNI are zero. (:issue:`2801`, :pull:`2808`)
+* Fix :py:func:`pvlib.irradiance.perez` to not return zero (instead NaN) when
+  airmass is NaN and either DHI or DNI are NaN. (:issue:`2801`, :pull:`2808`)
 
 Enhancements
 ~~~~~~~~~~~~
@@ -48,3 +51,4 @@ Contributors
 * Eesh Saxena (:ghuser:`eeshsaxena`)
 * Karl Hill (:ghuser:`karlhillx`)
 * Yonry Zhu (:ghuser:`yonryzhu`)
+* Mark Campanelli (:ghuser:`markcampanelli`)
diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 20df3966ba..41d45a4a42 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -1211,18 +1211,19 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     B = np.maximum(B, tools.cosd(85))
 
     # Calculate Diffuse POA from sky dome
-    term1 = 0.5 * (1 - F1) * (1 + tools.cosd(surface_tilt))
-    term2 = F1 * A / B
-    term3 = F2 * tools.sind(surface_tilt)
+    term1 = 0.5 * (1 - F1) * (1 + tools.cosd(surface_tilt)) # isotropic
+    term2 = F1 * A / B  # circumsolar
+    term3 = F2 * tools.sind(surface_tilt)  # horizon
 
     sky_diffuse = np.maximum(dhi * (term1 + term2 + term3), 0)
 
-    # we've preserved the input type until now, so don't ruin it!
+    # Use NaN airmass to coerce to zero values that are not otherwise NaN.
     airmass_nan_idx = np.logical_and(
         np.isnan(airmass), np.logical_not(
             np.logical_or(np.isnan(dhi), np.isnan(dni))
         )
     )
+    # we've preserved the input type until now, so don't ruin it!
     if isinstance(sky_diffuse, pd.Series):
         sky_diffuse[airmass_nan_idx] = 0
     else:

From acaa15f9e6846e7f2a8834331e0b2a39c1e17b13 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <mark.campanelli@gmail.com>
Date: Sun, 5 Jul 2026 17:28:38 -0600
Subject: [PATCH 12/14] Appease flake8 yet again

---
 pvlib/irradiance.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 41d45a4a42..23307f9328 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -1211,7 +1211,7 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     B = np.maximum(B, tools.cosd(85))
 
     # Calculate Diffuse POA from sky dome
-    term1 = 0.5 * (1 - F1) * (1 + tools.cosd(surface_tilt)) # isotropic
+    term1 = 0.5 * (1 - F1) * (1 + tools.cosd(surface_tilt))  # isotropic
     term2 = F1 * A / B  # circumsolar
     term3 = F2 * tools.sind(surface_tilt)  # horizon
 

From b3b5d75a7db72d7c6d24da2ad7d323c04b560367 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <mark.campanelli@gmail.com>
Date: Tue, 7 Jul 2026 08:23:16 -0600
Subject: [PATCH 13/14] Update docs/sphinx/source/whatsnew/v0.15.3.rst

Co-authored-by: Cliff Hansen <cwhanse@sandia.gov>
---
 docs/sphinx/source/whatsnew/v0.15.3.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/sphinx/source/whatsnew/v0.15.3.rst b/docs/sphinx/source/whatsnew/v0.15.3.rst
index 9c3618d5d5..459743f3ce 100644
--- a/docs/sphinx/source/whatsnew/v0.15.3.rst
+++ b/docs/sphinx/source/whatsnew/v0.15.3.rst
@@ -14,7 +14,7 @@ Deprecations
 
 Bug fixes
 ~~~~~~~~~
-* Fix :py:func:`pvlib.irradiance.perez` to not return NaN when both DHI and
+* Fix :py:func:`pvlib.irradiance.perez` to return 0 when both DHI and
   DNI are zero. (:issue:`2801`, :pull:`2808`)
 * Fix :py:func:`pvlib.irradiance.perez` to not return zero (instead NaN) when
   airmass is NaN and either DHI or DNI are NaN. (:issue:`2801`, :pull:`2808`)

From cb265433b2d5588f0b06d42e14360c0f10318f4e Mon Sep 17 00:00:00 2001
From: Mark Campanelli <mark.campanelli@gmail.com>
Date: Tue, 7 Jul 2026 08:24:00 -0600
Subject: [PATCH 14/14] Update docs/sphinx/source/whatsnew/v0.15.3.rst

Co-authored-by: Cliff Hansen <cwhanse@sandia.gov>
---
 docs/sphinx/source/whatsnew/v0.15.3.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/sphinx/source/whatsnew/v0.15.3.rst b/docs/sphinx/source/whatsnew/v0.15.3.rst
index 459743f3ce..7b1d2dbb15 100644
--- a/docs/sphinx/source/whatsnew/v0.15.3.rst
+++ b/docs/sphinx/source/whatsnew/v0.15.3.rst
@@ -16,7 +16,7 @@ Bug fixes
 ~~~~~~~~~
 * Fix :py:func:`pvlib.irradiance.perez` to return 0 when both DHI and
   DNI are zero. (:issue:`2801`, :pull:`2808`)
-* Fix :py:func:`pvlib.irradiance.perez` to not return zero (instead NaN) when
+* Fix :py:func:`pvlib.irradiance.perez` to return NaN when
   airmass is NaN and either DHI or DNI are NaN. (:issue:`2801`, :pull:`2808`)
 
 Enhancements