Update irradiance.aoi to use more reliable formula (#1191)

* clip to [-1, 1]

* add test

* whatsnew

* use new aoi formula

* revert changes to aoi_projection

* update test and whatsnew

* don't set name; break after first Series

* roll back to 2a0425b

* test improvements

Author: kandersolar

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

@@ -133,6 +133,8 @@ Bug fixes
   (:issue:`1065`, :pull:`1140`)
 * Reindl model fixed to generate sky_diffuse=0 when GHI=0.
   (:issue:`1153`, :pull:`1154`)
+* Fix floating point round-off issue in
+  :py:func:`~pvlib.irradiance.aoi_projection` (:issue:`1185`, :pull:`1191`)
 * Update GFS product names for GFS v16. (:issue:`1202`, :pull:`1203`)
 * Corrected methodology error in :py:func:`~pvlib.scaling.wvm`. Tracks with
   fix in PVLib for MATLAB. (:issue:`1206`, :pull:`1213`)

pvlib/irradiance.py

@@ -182,6 +182,9 @@ def aoi_projection(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth):
         tools.sind(surface_tilt) * tools.sind(solar_zenith) *
         tools.cosd(solar_azimuth - surface_azimuth))
 
+    # GH 1185
+    projection = np.clip(projection, -1, 1)
+
     try:
         projection.name = 'aoi_projection'
     except AttributeError:

pvlib/tests/test_irradiance.py

@@ -792,6 +792,27 @@ def test_aoi_and_aoi_projection(surface_tilt, surface_azimuth, solar_zenith,
     assert_allclose(aoi_projection, aoi_proj_expected, atol=1e-6)
 
 
+def test_aoi_projection_precision():
+    # GH 1185 -- test that aoi_projection does not exceed 1.0, and when
+    # given identical inputs, the returned projection is very close to 1.0
+
+    # scalars
+    zenith = 89.26778228223463
+    azimuth = 60.932028605997004
+    projection = irradiance.aoi_projection(zenith, azimuth, zenith, azimuth)
+    assert projection <= 1
+    assert np.isclose(projection, 1)
+
+    # arrays
+    zeniths = np.array([zenith])
+    azimuths = np.array([azimuth])
+    projections = irradiance.aoi_projection(zeniths, azimuths,
+                                            zeniths, azimuths)
+    assert all(projections <= 1)
+    assert all(np.isclose(projections, 1))
+    assert projections.dtype == np.dtype('float64')
+
+
 @pytest.fixture
 def airmass_kt():
     # disc algorithm stopped at am=12. test am > 12 for out of range behavior