diff --git a/Utilities/Optimisation.py b/Utilities/Optimisation.py
index 3034862..b5fdce3 100644
--- a/Utilities/Optimisation.py
+++ b/Utilities/Optimisation.py
@@ -27,7 +27,11 @@ def optimise_vertical_panel_pitch(c):
     # perform minimization
     initial_pitch = c["array"]["spacing"]
     result = minimize(
-        objective, initial_pitch, bounds=[(0, 20)], tol=1e-8, options={"eps": 1}
+        objective,
+        initial_pitch,
+        bounds=[(c["panel"]["dimensions"]["length"], 20)],
+        tol=1e-8,
+        options={"eps": 1},
     )
     optimal_pitch = result.x[0]
     c["array"]["spacing"] = optimal_pitch
diff --git a/Utilities/Shading.py b/Utilities/Shading.py
index 42ce992..5ac719b 100644
--- a/Utilities/Shading.py
+++ b/Utilities/Shading.py
@@ -2,6 +2,7 @@ import numpy as np
 import pandas as pd
 import logging
 import math
+import matplotlib.pyplot as pl
 
 import pvlib
 from pvlib.bifacial.pvfactors import pvfactors_timeseries
@@ -154,22 +155,15 @@ def calculate_energy_production_vertical(c):
     shading_row_rotation = 90
     surface_azimuth = 90  # east facing
     axis_tilt = 0
-    axis_azimuth = 0
+    axis_azimuth = 180
 
-    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
+    gcr = np.divide(c["panel"]["dimensions"]["length"], pitch)
+    gcr = min(1, gcr)
     logger.info(f"Ground coverage ratio: {gcr}")
 
     # use pvfactors bifacial modelling package
-    for row in range(no_of_rows):
+    POA_data = pd.Series(dtype=object)
+    for row in range(0, 3):
         result = pvfactors_timeseries(
             solar_zenith=solar_positions["apparent_zenith"],
             solar_azimuth=solar_positions["azimuth"],
@@ -183,55 +177,39 @@ def calculate_energy_production_vertical(c):
             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,
+            n_pvrows=3,
             index_observed_pvrow=row,
         )
+        # set negative values to 0
+        poa_front = result[2].clip(lower=0)
+        poa_rear = result[3].clip(lower=0)
+        poa_global = poa_front + poa_rear * c["panel"]["bifaciality"]
+        POA_data.at[row] = poa_global
 
-    # calculate irradiance on plane of array
-    poa_front = pvlib.irradiance.get_total_irradiance(
-        surface_tilt=90,
-        surface_azimuth=90,
-        solar_zenith=morning_projected_solar_zenith,
-        solar_azimuth=solar_positions["azimuth"],
-        dni=clearsky_data["dni"],
-        ghi=clearsky_data["ghi"],
-        dhi=clearsky_data["dhi"],
-        surface_type="urban",
-    )
-    # drop rows with poa_global NaN values
-    poa_front = poa_front.dropna(subset=["poa_global"])
-
-    poa_rear = pvlib.irradiance.get_total_irradiance(
-        surface_tilt=180 - 90,
-        surface_azimuth=90 + 180,
-        solar_zenith=afternoon_projected_solar_zenith,
-        solar_azimuth=solar_positions["azimuth"],
-        dni=clearsky_data["dni"],
-        ghi=clearsky_data["ghi"],
-        dhi=clearsky_data["dhi"],
-        surface_type="urban",
-    )
-    # drop rows with poa_global NaN values
-    poa_rear = poa_rear.dropna(subset=["poa_global"])
-
-    effective_front = poa_front["poa_global"] * (1 - morning_shaded_fraction)
-    effective_rear = (
-        poa_rear["poa_global"]
-        * (1 - afternoon_shaded_fraction)
-        * c["panel"]["bifaciality"]
-    )
-
-    total_hourly_irradiance = effective_front + effective_rear
-    system_size = c["panel"]["peak_power"] * no_of_panels
-    pdc0 = system_size
+    total_hourly_irradiance = POA_data.at[0] + POA_data.at[2] + (POA_data.at[1] * 20)
     gamma_pdc = c["panel"]["temperature_coefficient"]
     temp_cell = c["panel"]["nominal_operating_cell_temperature"]
-    pdc = pvlib.pvsystem.pvwatts_dc(
-        pdc0=pdc0,
+    p_row = no_of_panels_in_row * c["panel"]["peak_power"]
+    p_middle_rows = (no_of_panels - 2 * no_of_panels_in_row) * c["panel"]["peak_power"]
+    pdc_first_row = pvlib.pvsystem.pvwatts_dc(
+        pdc0=p_row,
         gamma_pdc=gamma_pdc,
         temp_cell=temp_cell,
-        g_poa_effective=total_hourly_irradiance,
+        g_poa_effective=POA_data.at[0],
     )
+    pdc_last_row = pvlib.pvsystem.pvwatts_dc(
+        pdc0=p_row,
+        gamma_pdc=gamma_pdc,
+        temp_cell=temp_cell,
+        g_poa_effective=POA_data.at[2],
+    )
+    pdc_middle_rows = pvlib.pvsystem.pvwatts_dc(
+        pdc0=p_middle_rows,
+        gamma_pdc=gamma_pdc,
+        temp_cell=temp_cell,
+        g_poa_effective=POA_data.at[1],
+    )
+    pdc = pdc_first_row + pdc_last_row + pdc_middle_rows
 
     total_hourly_energy = pdc * 15 / 60 / 1e3  # convert to kWh
 
@@ -259,7 +237,7 @@ def calculate_energy_production_horizontal(c):
     shaded_row_rotation = 0
     shading_row_rotation = 0
     axis_tilt = 0
-    axis_azimuth = 90  # south facing
+    axis_azimuth = 270  # south facing surface
 
     projected_solar_zenith = pvlib.shading.projected_solar_zenith_angle(
         solar_zenith=solar_positions["apparent_zenith"],