diff --git a/docs/sphinx/source/whatsnew/v0.9.4.rst b/docs/sphinx/source/whatsnew/v0.9.4.rst
index 1d890e5512..8a67c201f0 100644
--- a/docs/sphinx/source/whatsnew/v0.9.4.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.4.rst
@@ -20,6 +20,8 @@ Enhancements
 Bug fixes
 ~~~~~~~~~
 
+* Fixed bug in :py:func:`pvlib.shading.masking_angle` and :py:func:`pvlib.bifacial.infinite_sheds._ground_angle`
+  where zero ``gcr`` input caused a ZeroDivisionError (:issue:`1576`, :pull:`1589`)
 
 Testing
 ~~~~~~~
@@ -44,5 +46,5 @@ Contributors
 * Christian Orner (:ghuser:`chrisorner`)
 * Saurabh Aneja (:ghuser:`spaneja`)
 * Marcus Boumans (:ghuser:`bowie2211`)
+* Karel De Brabandere (:ghuser:`kdebrab`)
 * Naman Priyadarshi (:ghuser:`Naman-Priyadarshi`)
-
diff --git a/pvlib/bifacial/infinite_sheds.py b/pvlib/bifacial/infinite_sheds.py
index 3dd37e3b05..47cf7c2d4d 100644
--- a/pvlib/bifacial/infinite_sheds.py
+++ b/pvlib/bifacial/infinite_sheds.py
@@ -216,8 +216,8 @@ def _ground_angle(x, surface_tilt, gcr):
     #  :         \  v      *-.\
     #  :          \<-----P---->\
 
-    x1 = x * sind(surface_tilt)
-    x2 = (x * cosd(surface_tilt) + 1 / gcr)
+    x1 = gcr * x * sind(surface_tilt)
+    x2 = gcr * x * cosd(surface_tilt) + 1
     psi = np.arctan2(x1, x2)  # do this first because it handles 0 / 0
     return np.rad2deg(psi)
 
diff --git a/pvlib/shading.py b/pvlib/shading.py
index 9479eb1739..01d725207c 100644
--- a/pvlib/shading.py
+++ b/pvlib/shading.py
@@ -52,8 +52,8 @@ def masking_angle(surface_tilt, gcr, slant_height):
     # The original equation (8 in [1]) requires pitch and collector width,
     # but it's easy to non-dimensionalize it to make it a function of GCR
     # by factoring out B from the argument to arctan.
-    numerator = (1 - slant_height) * sind(surface_tilt)
-    denominator = 1/gcr - (1 - slant_height) * cosd(surface_tilt)
+    numerator = gcr * (1 - slant_height) * sind(surface_tilt)
+    denominator = 1 - gcr * (1 - slant_height) * cosd(surface_tilt)
     phi = np.arctan(numerator / denominator)
     return np.degrees(phi)
 
diff --git a/pvlib/tests/bifacial/test_infinite_sheds.py b/pvlib/tests/bifacial/test_infinite_sheds.py
index 05d1bc58db..696d10568c 100644
--- a/pvlib/tests/bifacial/test_infinite_sheds.py
+++ b/pvlib/tests/bifacial/test_infinite_sheds.py
@@ -106,6 +106,14 @@ def test__ground_angle(test_system):
     assert np.allclose(angles, expected_angles)
 
 
+def test__ground_angle_zero_gcr():
+    surface_tilt = 30.0
+    x = np.array([0.0, 0.5, 1.0])
+    angles = infinite_sheds._ground_angle(x, surface_tilt, 0)
+    expected_angles = np.array([0, 0, 0])
+    assert np.allclose(angles, expected_angles)
+
+
 def test__vf_row_ground(test_system):
     ts, _, _ = test_system
     x = np.array([0., 0.5, 1.0])
diff --git a/pvlib/tests/test_shading.py b/pvlib/tests/test_shading.py
index 8a9fd46a69..0558cb76ad 100644
--- a/pvlib/tests/test_shading.py
+++ b/pvlib/tests/test_shading.py
@@ -45,6 +45,13 @@ def test_masking_angle_scalar(surface_tilt, masking_angle):
         assert np.isclose(masking_angle_actual, angle)
 
 
+def test_masking_angle_zero_gcr(surface_tilt):
+    # scalar inputs and outputs, including zero
+    for tilt in surface_tilt:
+        masking_angle_actual = shading.masking_angle(tilt, 0, 0.25)
+        assert np.isclose(masking_angle_actual, 0)
+
+
 def test_masking_angle_passias_series(surface_tilt, average_masking_angle):
     # pandas series inputs and outputs
     masking_angle_actual = shading.masking_angle_passias(surface_tilt, 0.5)