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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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

From 3561960c61c9728c6e3d0befd3b2be28ba6f192e Mon Sep 17 00:00:00 2001
From: Mark Campanelli <mark.campanelli@gmail.com>
Date: Sat, 3 Oct 2026 11:34:14 -0600
Subject: [PATCH 15/17] Better align perez and perez_driesse

---
 pvlib/irradiance.py | 71 ++++++++++++++++++---------------------------
 1 file changed, 29 insertions(+), 42 deletions(-)

diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index cc30361b07..2e51eba19a 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -1157,7 +1157,7 @@ def king(surface_tilt, dhi, ghi, solar_zenith):
 
 
 def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
-          solar_zenith, solar_azimuth, airmass,
+          solar_zenith, solar_azimuth, airmass=None,
           model='allsitescomposite1990', return_components=False):
     '''
     Determine diffuse irradiance from the sky on a tilted surface using
@@ -1207,11 +1207,13 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     solar_azimuth : numeric
         Solar azimuth angle. See :term:`solar_azimuth`. [°]
 
-    airmass : numeric
-        Relative (not pressure-corrected) airmass values. If AM is a
-        DataFrame it must be of the same size as all other DataFrame
-        inputs. AM must be >=0 (careful using the 1/cos(z) model of AM
-        generation). [unitless]
+    airmass : numeric, optional
+        Relative (not pressure-corrected) airmass values. If ``airmass`` is a
+        DataFrame, then it must be of the same size as all other DataFrame
+        inputs. AM must be >=0 (careful using the 1/cos(z) model of AM 
+        generation). If not provided or None, then the kastenyoung1989 model
+        is used, which corresponds to the calibration of the F1 and F2 
+        coefficients in the original work by Perez. [unitless]
 
     model : string, default 'allsitescomposite1990'
         A string which selects the desired set of Perez coefficients. If
@@ -1281,36 +1283,20 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
        Perez Diffuse Radiation Model". SAND88-7030
     '''
 
-    kappa = 1.041  # for solar_zenith in radians
-    z = np.radians(solar_zenith)  # convert to radians
-
-    # 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.
-    # Assumes:
+    # FIXME Verify _calc_zeta 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):
-        eps = eps.values
+    delta = _calc_delta(dhi, dni_extra, solar_zenith, airmass)
+    zeta, z = _calc_zeta(dhi, dni, solar_zenith)
 
     # Perez et al define clearness bins according to the following
     # rules. 1 = overcast ... 8 = clear (these names really only make
     # sense for small zenith angles, but...) these values will
     # eventually be used as indicies for coeffecient look ups
-    ebin = np.digitize(eps, (0., 1.065, 1.23, 1.5, 1.95, 2.8, 4.5, 6.2))
+    ebin = np.digitize(zeta, (0., 1.065, 1.23, 1.5, 1.95, 2.8, 4.5, 6.2))
     ebin = np.array(ebin)  # GH 642
-    ebin[np.isnan(eps)] = 0
+    ebin[np.isnan(zeta)] = 0
 
     # correct for 0 indexing in coeffecient lookup
     # later, ebin = -1 will yield nan coefficients
@@ -1320,7 +1306,7 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     # in a subfunction to clean up the code.
     F1c, F2c = _get_perez_coefficients(model)
 
-    # results in invalid eps (ebin = -1) being mapped to nans
+    # results in invalid zeta (ebin = -1) being mapped to nans
     nans = np.array([np.nan, np.nan, np.nan])
     F1c = np.vstack((F1c, nans))
     F2c = np.vstack((F2c, nans))
@@ -1382,12 +1368,13 @@ def _calc_delta(dhi, dni_extra, solar_zenith, airmass=None):
     Compute the delta parameter, which represents sky dome "brightness"
     in the Perez and Perez-Driesse models.
 
-    Helper function for perez_driesse transposition.
+    Helper function for perez and perez_driesse transposition.
     '''
     if airmass is None:
         # use the same airmass model as in the original perez work
-        airmass = atmosphere.get_relative_airmass(solar_zenith,
-                                                  'kastenyoung1989')
+        airmass = atmosphere.get_relative_airmass(
+            solar_zenith, 'kastenyoung1989'
+        )
 
     max_airmass = atmosphere.get_relative_airmass(90, 'kastenyoung1989')
     airmass = np.where(solar_zenith >= 90, max_airmass, airmass)
@@ -1398,9 +1385,9 @@ def _calc_delta(dhi, dni_extra, solar_zenith, airmass=None):
 def _calc_zeta(dhi, dni, zenith):
     '''
     Compute the zeta parameter, which represents sky dome "clearness"
-    in the Perez-Driesse model.
+    in the Perez and Perez-Driesse model.
 
-    Helper function for perez_driesse transposition.
+    Helper function for perez and perez_driesse transposition.
     '''
     dhi = np.asarray(dhi)
     dni = np.asarray(dni)
@@ -1414,10 +1401,10 @@ def _calc_zeta(dhi, dni, zenith):
 
     # then apply the kappa correction in a manner analogous to eq. 7
     kappa = 1.041
-    kterm = kappa * np.radians(zenith) ** 3
-    zeta = zeta / (1 - kterm * (zeta - 1))
+    z = np.radians(zenith)
+    zeta = zeta / (1 - kappa * z ** 3 * (zeta - 1))
 
-    return zeta
+    return zeta, z
 
 
 def _f(i, j, zeta):
@@ -1491,9 +1478,11 @@ def perez_driesse(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
 
     airmass : numeric, optional
         Relative (not pressure-corrected) airmass values. If ``airmass`` is a
-        DataFrame it must be of the same size as all other DataFrame
-        inputs. AM must be >=0 (careful using the 1/cos(z) model of AM
-        generation). [unitless]
+        DataFrame, then it must be of the same size as all other DataFrame
+        inputs. AM must be >=0 (careful using the 1/cos(z) model of AM 
+        generation). If not provided or None, then the kastenyoung1989 model
+        is used, which corresponds to the calibration of the F1 and F2 
+        coefficients in the original work by Perez. [unitless]
 
     return_components: bool (optional, default=False)
         Flag used to decide whether to return the calculated diffuse components
@@ -1550,9 +1539,7 @@ def perez_driesse(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     # Contributed by Anton Driesse (@adriesse), PV Performance Labs. Oct., 2023
 
     delta = _calc_delta(dhi, dni_extra, solar_zenith, airmass)
-    zeta = _calc_zeta(dhi, dni, solar_zenith)
-
-    z = np.radians(solar_zenith)
+    zeta, z = _calc_zeta(dhi, dni, solar_zenith)
 
     F1 = _f(1, 1, zeta) + _f(1, 2, zeta) * delta + _f(1, 3, zeta) * z
     F2 = _f(2, 1, zeta) + _f(2, 2, zeta) * delta + _f(2, 3, zeta) * z

From 8d3808d2eb3ba96858208ee33a143575e99df716 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <mark.campanelli@gmail.com>
Date: Sat, 3 Oct 2026 16:16:31 -0600
Subject: [PATCH 16/17] Remove new default value for airmass in perez

---
 pvlib/irradiance.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 2e51eba19a..0a9e67e5f9 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -1156,9 +1156,9 @@ def king(surface_tilt, dhi, ghi, solar_zenith):
     return poa_sky_diffuse
 
 
-def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
-          solar_zenith, solar_azimuth, airmass=None,
-          model='allsitescomposite1990', return_components=False):
+def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra, solar_zenith,
+          solar_azimuth, airmass, model='allsitescomposite1990',
+          return_components=False):
     '''
     Determine diffuse irradiance from the sky on a tilted surface using
     one of the Perez models.
@@ -1207,13 +1207,13 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     solar_azimuth : numeric
         Solar azimuth angle. See :term:`solar_azimuth`. [°]
 
-    airmass : numeric, optional
+    airmass : numeric or None
         Relative (not pressure-corrected) airmass values. If ``airmass`` is a
         DataFrame, then it must be of the same size as all other DataFrame
         inputs. AM must be >=0 (careful using the 1/cos(z) model of AM 
-        generation). If not provided or None, then the kastenyoung1989 model
-        is used, which corresponds to the calibration of the F1 and F2 
-        coefficients in the original work by Perez. [unitless]
+        generation). If None, then the 'kastenyoung1989' model is used, which
+        corresponds to the calibration of the F1 and F2 coefficients in the
+        original work by Perez. [unitless]
 
     model : string, default 'allsitescomposite1990'
         A string which selects the desired set of Perez coefficients. If

From 14090bd2a47f24be1f98c257cd1934a5dec68ed6 Mon Sep 17 00:00:00 2001
From: Mark Campanelli <mark.campanelli@gmail.com>
Date: Sat, 3 Oct 2026 16:24:42 -0600
Subject: [PATCH 17/17] Appease linter

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

diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 0a9e67e5f9..55f9b21bc1 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -1210,7 +1210,7 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra, solar_zenith,
     airmass : numeric or None
         Relative (not pressure-corrected) airmass values. If ``airmass`` is a
         DataFrame, then it must be of the same size as all other DataFrame
-        inputs. AM must be >=0 (careful using the 1/cos(z) model of AM 
+        inputs. AM must be >=0 (careful using the 1/cos(z) model of AM
         generation). If None, then the 'kastenyoung1989' model is used, which
         corresponds to the calibration of the F1 and F2 coefficients in the
         original work by Perez. [unitless]
@@ -1479,9 +1479,9 @@ def perez_driesse(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     airmass : numeric, optional
         Relative (not pressure-corrected) airmass values. If ``airmass`` is a
         DataFrame, then it must be of the same size as all other DataFrame
-        inputs. AM must be >=0 (careful using the 1/cos(z) model of AM 
+        inputs. AM must be >=0 (careful using the 1/cos(z) model of AM
         generation). If not provided or None, then the kastenyoung1989 model
-        is used, which corresponds to the calibration of the F1 and F2 
+        is used, which corresponds to the calibration of the F1 and F2
         coefficients in the original work by Perez. [unitless]
 
     return_components: bool (optional, default=False)