limit cos(aoi) in diffuse sky functions (#535)

* limit cos(aoi) in diffuse sky functions

* no limit on cos_solar_zenith in hb calculation. change proj limit to -1

* reimplement without projection_minimum kwarg

* fix np.minimum vs. np.maximum

* update tutorial

* add sky diffuse test for zenith close to 90

* fix botched whatsnew merge [skip ci]

Author: wholmgren

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

@@ -111,6 +111,11 @@ Bug fixes
   floats and NumPy arrays. (:issue:`508`)
 * Make GitHub recognize the license, add AUTHORS.md, clarify shared copyright.
   (:issue:`503`)
+* Fix issue with non-zero direct irradiance contribution to Reindl, Klucher,
+  and Hay-Davies diffuse sky algorithms when the sun is behind the array.
+  (:issue:`526`)
+* Fix issue with dividing by near-0 cos(solar_zenith) values in Reindl and
+  Hay-Davies diffuse sky algorithms. (:issue:`432`)
 * Fix argument order of longitude and latitude when querying weather forecasts
   by lonlat bounding box (:issue:`521`)
 

docs/tutorials/irradiance.ipynb

@@ -256,7 +256,7 @@ def poa_horizontal_ratio(surface_tilt, surface_azimuth,
 
     cos_solar_zenith = tools.cosd(solar_zenith)
 
-    # ratio of titled and horizontal beam irradiance
+    # ratio of tilted and horizontal beam irradiance
     ratio = cos_poa_zen / cos_solar_zenith
 
     try:
@@ -746,6 +746,7 @@ def klucher(surface_tilt, surface_azimuth, dhi, ghi, solar_zenith,
     # zenith angle with respect to panel normal.
     cos_tt = aoi_projection(surface_tilt, surface_azimuth,
                             solar_zenith, solar_azimuth)
+    cos_tt = np.maximum(cos_tt, 0)  # GH 526
 
     F = 1 - ((dhi / ghi) ** 2)
     try:
@@ -836,8 +837,9 @@ def haydavies(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     if projection_ratio is None:
         cos_tt = aoi_projection(surface_tilt, surface_azimuth,
                                 solar_zenith, solar_azimuth)
+        cos_tt = np.maximum(cos_tt, 0)  # GH 526
         cos_solar_zenith = tools.cosd(solar_zenith)
-        Rb = cos_tt / cos_solar_zenith
+        Rb = cos_tt / np.maximum(cos_solar_zenith, 0.01745)  # GH 432
     else:
         Rb = projection_ratio
 
@@ -932,11 +934,13 @@ def reindl(surface_tilt, surface_azimuth, dhi, dni, ghi, dni_extra,
 
     cos_tt = aoi_projection(surface_tilt, surface_azimuth,
                             solar_zenith, solar_azimuth)
+    cos_tt = np.maximum(cos_tt, 0)  # GH 526
 
+    # do not apply cos(zen) limit here (needed for HB below)
     cos_solar_zenith = tools.cosd(solar_zenith)
 
     # ratio of titled and horizontal beam irradiance
-    Rb = cos_tt / cos_solar_zenith
+    Rb = cos_tt / np.maximum(cos_solar_zenith, 0.01745)  # GH 432
 
     # Anisotropy Index
     AI = dni / dni_extra

pvlib/irradiance.py

@@ -181,7 +181,9 @@ def test_klucher_series(irrad_data, ephem_data):
     result = irradiance.klucher(40, 180, irrad_data['dhi'], irrad_data['ghi'],
                                 ephem_data['apparent_zenith'],
                                 ephem_data['azimuth'])
-    assert_allclose(result, [0, 37.446276, 109.209347, 56.965916], atol=1e-4)
+    # pvlib matlab 1.4 does not contain the max(cos_tt, 0) correction
+    # so, these values are different
+    assert_allclose(result, [0., 36.789794, 109.209347, 56.965916], atol=1e-4)
     # expect same result for np.array and pd.Series
     expected = irradiance.klucher(
         40, 180, irrad_data['dhi'].values, irrad_data['ghi'].values,
@@ -195,15 +197,17 @@ def test_haydavies(irrad_data, ephem_data, dni_et):
                          dni_et,
                          ephem_data['apparent_zenith'],
                          ephem_data['azimuth'])
-    assert_allclose(result, [0, 14.967008, 102.994862, 33.190865], atol=1e-4)
+    # values from matlab 1.4 code
+    assert_allclose(result, [0, 27.1775, 102.9949, 33.1909], atol=1e-4)
 
 
 def test_reindl(irrad_data, ephem_data, dni_et):
     result = irradiance.reindl(40, 180, irrad_data['dhi'], irrad_data['dni'],
                       irrad_data['ghi'], dni_et,
                       ephem_data['apparent_zenith'],
                       ephem_data['azimuth'])
-    assert_allclose(result, [np.nan, 15.730664, 104.131724, 34.166258], atol=1e-4)
+    # values from matlab 1.4 code
+    assert_allclose(result, [np.nan, 27.9412, 104.1317, 34.1663], atol=1e-4)
 
 
 def test_king(irrad_data, ephem_data):
@@ -265,6 +269,16 @@ def test_perez_arrays(irrad_data, ephem_data, dni_et, relative_airmass):
     assert_allclose(out, expected, atol=1e-2)
 
 
+@pytest.mark.parametrize('model', ['isotropic', 'klucher', 'haydavies',
+                                   'reindl', 'king', 'perez'])
+def test_sky_diffuse_zenith_close_to_90(model):
+    # GH 432
+    sky_diffuse = irradiance.get_sky_diffuse(
+        30, 180, 89.999, 230,
+        dni=10, ghi=51, dhi=50, dni_extra=1360, airmass=12, model=model)
+    assert sky_diffuse < 100
+
+
 def test_liujordan():
     expected = pd.DataFrame(np.
         array([[863.859736967, 653.123094076, 220.65905025]]),