diff --git a/Utilities/Shading.py b/Utilities/Shading.py
index 62136d7..42ce992 100644
--- a/Utilities/Shading.py
+++ b/Utilities/Shading.py
@@ -4,6 +4,7 @@ import logging
 import math
 
 import pvlib
+from pvlib.bifacial.pvfactors import pvfactors_timeseries
 
 from Utilities.Processes import calculate_no_of_panels, calculate_required_system_size
 
@@ -143,57 +144,48 @@ def calculate_energy_production_vertical(c):
     no_of_rows = np.unique(panel_coordinates["y"]).shape[0]
     no_of_shaded_rows = no_of_rows - 1
 
+    no_of_panels_in_row = np.unique(panel_coordinates["x"]).shape[0]
+
     collector_width = c["panel"]["dimensions"]["length"]
     # calculate delta between unique y coordinates of panels to get pitch
     pitch = np.unique(panel_coordinates["y"])[1] - np.unique(panel_coordinates["y"])[0]
     surface_to_axis_offset = 0
     shaded_row_rotation = 90
     shading_row_rotation = 90
+    surface_azimuth = 90  # east facing
     axis_tilt = 0
-    axis_azimuth = 180
+    axis_azimuth = 0
 
-    morning_projected_solar_zenith = pvlib.shading.projected_solar_zenith_angle(
-        solar_zenith=solar_positions["apparent_zenith"],
-        solar_azimuth=solar_positions["azimuth"],
-        axis_azimuth=axis_azimuth,
-        axis_tilt=axis_tilt,
+    ground_coverage = (
+        no_of_panels
+        * c["panel"]["dimensions"]["width"]
+        * c["panel"]["dimensions"]["thickness"]
     )
+    land_area = (
+        c["environment"]["roof"]["dimensions"]["width"]
+        * c["environment"]["roof"]["dimensions"]["length"]
+    )
+    gcr = ground_coverage / land_area
+    logger.info(f"Ground coverage ratio: {gcr}")
 
-    morning_shaded_fraction = pvlib.shading.shaded_fraction1d(
-        solar_zenith=morning_projected_solar_zenith,
-        solar_azimuth=solar_positions["azimuth"],
-        axis_azimuth=axis_azimuth,
-        shaded_row_rotation=shaded_row_rotation,
-        shading_row_rotation=shading_row_rotation,
-        collector_width=collector_width,
-        pitch=pitch,
-        surface_to_axis_offset=surface_to_axis_offset,
-        axis_tilt=axis_tilt,
-    )
-    morning_shaded_fraction = morning_shaded_fraction * no_of_shaded_rows / no_of_rows
-
-    afternoon_projected_solar_zenith = pvlib.shading.projected_solar_zenith_angle(
-        solar_zenith=solar_positions["apparent_zenith"],
-        solar_azimuth=solar_positions["azimuth"],
-        axis_azimuth=axis_azimuth + 180,
-        axis_tilt=axis_tilt,
-    )
-    afternoon_shaded_fraction = pvlib.shading.shaded_fraction1d(
-        solar_zenith=afternoon_projected_solar_zenith,
-        solar_azimuth=solar_positions["azimuth"],
-        axis_azimuth=axis_azimuth + 180,
-        shaded_row_rotation=shaded_row_rotation,
-        shading_row_rotation=shading_row_rotation,
-        collector_width=collector_width,
-        pitch=pitch,
-        surface_to_axis_offset=surface_to_axis_offset,
-        axis_tilt=axis_tilt,
-    )
-    afternoon_shaded_fraction = (
-        afternoon_shaded_fraction * no_of_shaded_rows / no_of_rows
-    )
-
-    logger.info("Shaded fraction calculated for solar positions")
+    # use pvfactors bifacial modelling package
+    for row in range(no_of_rows):
+        result = pvfactors_timeseries(
+            solar_zenith=solar_positions["apparent_zenith"],
+            solar_azimuth=solar_positions["azimuth"],
+            surface_azimuth=surface_azimuth,
+            surface_tilt=90,
+            axis_azimuth=axis_azimuth,
+            timestamps=solar_positions.index,
+            dni=clearsky_data["dni"],
+            dhi=clearsky_data["dhi"],
+            gcr=gcr,
+            pvrow_height=c["panel"]["dimensions"]["length"],
+            pvrow_width=c["panel"]["dimensions"]["width"] * no_of_panels_in_row,
+            albedo=0.2,
+            n_pvrows=no_of_rows,
+            index_observed_pvrow=row,
+        )
 
     # calculate irradiance on plane of array
     poa_front = pvlib.irradiance.get_total_irradiance(