shading.py

"""
Contains all functions to do with shading and irradiance
"""

import time
import pandas as pd
import numpy as np
from numpy import linalg
from __utils__ import Array, Module, Cell
from SunAngles import findAngles
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns
from tqdm import tqdm
from shapely import Polygon
import warnings

__all__ = (
    "calculate_self_shading",
    "calculate_crop_shading",
)


def calculate_self_shading(
    pv_array: Array,
    solar_irradiance: np.ndarray,  # I'm assumng solar irradiance from renewables ninja is DNI? Also assuming packaged in a nice way
    Time: tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray],
    spectral_calculation: bool = False,
    plot_graphs: bool = False,
) -> dict[
    dict[Module, np.ndarray],
    dict[Module, np.ndarray],
    dict[Module, np.ndarray],
    dict[Module, np.ndarray] | None,
    dict[Module, np.ndarray] | None,
    dict[Module, np.ndarray] | None,
]:
    """

    Calculates panel on panel shading of a pv array. Returns a dictionary of dictionaries.

    Parameters
    ----------
    - pv_array: Array
        Holds the attributes of the pv array to be modelled
    - solar_irradiance: ndarray
        Solar irridiance for each point in time
    - Time: tuple[ndarray, ndarray, ndarray]
        Tuple that holds all the years, months, days and hours of each individaul time point
        to be calculated. Most importantly, for each hour the respective year, month and day
        must be added to set of arrays. The length of each array in the Time tuple must be the same!
    - spectral calculations: bool, optional
        boolean that says when to calculate the spoectral analysis. Currently means nothing
    - plot_graphs: bool, optional
        Set to false. If set to true will plot: 1 -  the angles of the sun vs the normal,
        2 - plotting the middle of the shadow onto the plane of the panel, 3 - the shadow
        onto the individual modules

    Returns
    ----------
    - illumination: dict
        Dictionary that is keyed to each module in one row and has cell illumination for each time point
    - front row illumination: dict
        dictionary that is keyed to each in one row module and has cell illumination for each time
        point but no panel-panel shading. This is for the front row of the pv array
    - shading fraction: dict
        dictionary keyed to each module in a row and has cell fractional shading
        If the pv_array is bifacial then also emits
    - back illumination: dict
        dictionary that is keyed to each module in one row and has cell
        illumination of the back of the panel. Returned only if the module type is "Bifacial"
    - back row illumination: dict
        dictionary that is keyed to each module in one row and has cell
        illumination of the back of the panel but no shading. Returned only if the module type is "Bifacial"
    - back shading fraction: dict
        dictionary that is keyed to each module in one row and has cell fractional shading on the back.
        Returned only if the module type is "Bifacial"

    """

    # Assuming all modules are at the same angle(!): find angle rays make against the panel.
    #   1. Find Angles
    #   1.1 Define vector width (3, 1) as vector that goes parallel to base of panel
    #   1.2 Define vector length (3, 1) as vector that goes parallel to upright edge of panel
    #   1.3 Define sun ray array (N, 3)as from angles given by the findAngles function
    #   1.4 Take np.arcos of matrix multiply of vector 1 and sun ray to find array x angles (1, N)
    #   1.5 Take np.arcos of matrix multipy of vector 2 and sun ray to find array of y angles (1, N)

    # Defining vectors of the panel that go along the width and length of the module
    vector_width = np.array(
        [
            [np.sin(pv_array.orientation - np.pi / 2)],
            [np.cos(pv_array.orientation - np.pi / 2)],
            [0],
        ]
    )
    vector_length = np.array(
        [
            [np.sin(np.pi + pv_array.orientation) * np.cos(pv_array.tilt)],
            [np.cos(np.pi + pv_array.orientation) * np.cos(pv_array.tilt)],
            [np.sin(pv_array.tilt)],
        ]
    )

    # defing the normal of the module
    # normal used to define the "front"
    # if light is hitting with an angle greater than pi/2 then its striking the back of the module
    # Use and if statement to get effeciency on the back.
    normal = np.cross(vector_width.flatten(), vector_length.flatten())
    sunEl_rad, azimuth_rad = findAngles(
        pv_array.longitude,
        pv_array.latitude,
        Time[0],
        Time[1],
        Time[2],
        Time[3],
        pv_array.TimeZone,
    )
    # sunEl_rad = np.array([0.12517557 , 0.0485066])
    # azimuth_rad = np.array([4.95470605, 4.96098135])
    SunRay = np.array(
        [
            np.cos(sunEl_rad)
            * np.cos(
                np.pi / 2 - azimuth_rad
            ),  # Why np.pi/2 - azimuth  you ask? North is the y axis and azimuth is measured clockwise
            np.cos(sunEl_rad) * np.sin(np.pi / 2 - azimuth_rad),
            np.sin(sunEl_rad),
        ]
    )
    SunRay = SunRay.transpose()

    # finding angle between
    y_flattened_ray = SunRay - np.outer(np.dot(SunRay, vector_width), vector_width)
    y_flattened_ray = y_flattened_ray.transpose() / np.linalg.norm(
        y_flattened_ray, axis=1
    )
    x_flattened_ray = SunRay - np.outer(np.dot(SunRay, vector_length), vector_length)
    x_flattened_ray = x_flattened_ray.transpose() / np.linalg.norm(
        x_flattened_ray, axis=1
    )

    xangles: np.ndarray = np.acos(np.dot(x_flattened_ray.transpose(), vector_width))
    yangles: np.ndarray = np.acos(np.dot(y_flattened_ray.transpose(), vector_length))
    xangles = xangles.flatten()
    yangles = yangles.flatten()

    # Important for when we do inverse 1/tan. Just shifting so the tan values to stop NaN answers
    xangles[np.isclose(np.remainder(xangles, np.pi), 0)] += 1e-16
    yangles[np.isclose(np.remainder(yangles, np.pi), 0)] += 1e-16
    xangles[np.isclose(np.remainder(xangles, np.pi / 2), 0)] -= 1e-16
    yangles[np.isclose(np.remainder(yangles, np.pi / 2), 0)] -= 1e-16

    # calculating the x and z starting values for the front row.
    # The 1000 is to convert between meters and mm
    z_shift = (
        1000 * pv_array.ySpacing * np.sin(pv_array.tilt) / pv_array.Module.cell.xlength
    )  # z wise distance between rows (with the module being at z=0)

    y_shift = (
        1000 * pv_array.ySpacing * np.cos(pv_array.tilt) / pv_array.Module.cell.xlength
    )  # y wise distance between rows (with the module being at y=0)
    # IMPORTANT: I'm measuring things is cell width

    strike_angles = np.arccos(np.dot(SunRay, normal))

    # 2. Understanding where the light is hitting
    #   if strike_angle < pi/2 , light hits front. If < pi/2 light hits back side
    #   do something like strike_bool =  np ones like strike angles and then boolean to set back angles to 0 (aka ignore)
    #   thus strike_bool = 1 when light is hitting and 0 when its not
    direct = np.ones_like(strike_angles)

    direct[strike_angles >= np.pi / 2] = 0  # equals 0 if light hits backside
    direct[sunEl_rad <= 0] = (
        -1
    )  # equals -1 if sun is under the horizon. Presuming no elevation or on a hill.
    # There's probably some spiffy equation that says what altitude you have to be so see the sun at some angle but...
    #
    #
    # 3. Create the shadows
    #   3.1. Define the row ahead of the module
    #   3.2. Project the row down to the z=0 plane aka the plane that the shaded module lies in
    #   3.3. If the array is bifacial then repeat but with a row below and project upwards
    #   3.4. For each point in time: check direct to see where light is shining. 1 = front, 0 = back, -1 = under horizon
    #   3.5. If direct = 1 then find where the shadow strikes if it strikes the row behind
    #   3.6. If direct = 0 and module is bifacial do shading on the back

    # Defining row length. Really just a useful parameter
    row_length = (
        pv_array.columns * pv_array.Module.cells_wide
        + (pv_array.columns - 1)
        * 1000
        * pv_array.xSpacing
        / pv_array.Module.cell.xlength
    )

    # Definig corners of row which's shadow that will be projected down
    bottom_left = np.array([0, -y_shift, z_shift])
    bottom_right = np.array(
        [
            row_length,
            0 - y_shift,
            z_shift,
        ]
    )
    top_left = np.array([0, pv_array.Module.cells_long_in_width - y_shift, z_shift])
    top_right = np.array(
        [
            row_length,
            pv_array.Module.cells_long_in_width - y_shift,
            z_shift,
        ]
    )

    # Projecting each corner onto the module plane (z=0) This creates the points of the shadow at each time point
    proj_bottom_left = z_shift / np.tan([xangles, yangles])
    proj_bottom_left[0, :] = -proj_bottom_left[0, :] + bottom_left[0].astype(float)
    proj_bottom_left[1, :] = -proj_bottom_left[1, :] + bottom_left[1].astype(float)

    proj_bottom_right = z_shift / np.tan([xangles, yangles])
    proj_bottom_right[0, :] = -proj_bottom_right[0, :] + bottom_right[0].astype(float)
    proj_bottom_right[1, :] = -proj_bottom_right[1, :] + bottom_right[1].astype(float)

    proj_top_left = z_shift / np.tan([xangles, yangles])
    proj_top_left[0, :] = -proj_top_left[0, :] + top_left[0].astype(float)
    proj_top_left[1, :] = -proj_top_left[1, :] + top_left[1].astype(float)

    proj_top_right = z_shift / np.tan([xangles, yangles])
    proj_top_right[0, :] = -proj_top_right[0, :] + top_right[0].astype(float)
    proj_top_right[1, :] = -proj_top_right[1, :] + top_right[1].astype(float)

    # This is just for bifacial panels and generates the shadow cast on the panel from the back.
    # We do the same thing but from the back
    if pv_array.Module.type == "Bifacial":
        back_bottom_left = np.array([0, y_shift, -z_shift])
        back_bottom_right = np.array(
            [
                row_length,
                y_shift,
                -z_shift,
            ]
        )
        back_top_left = np.array(
            [0, pv_array.Module.cells_long_in_width + y_shift, -z_shift]
        )
        back_top_right = np.array(
            [
                row_length,
                pv_array.Module.cells_long_in_width + y_shift,
                -z_shift,
            ]
        )

        back_proj_bottom_left = z_shift / np.tan([xangles, yangles])
        back_proj_bottom_left[0, :] = back_proj_bottom_left[0, :] + back_bottom_left[
            0
        ].astype(float)
        back_proj_bottom_left[1, :] = -back_proj_bottom_left[1, :] + back_bottom_left[
            1
        ].astype(float)

        back_proj_bottom_right = z_shift / np.tan([xangles, yangles])
        back_proj_bottom_right[0, :] = back_proj_bottom_right[0, :] + back_bottom_right[
            0
        ].astype(float)
        back_proj_bottom_right[1, :] = -back_proj_bottom_right[
            1, :
        ] + back_bottom_right[1].astype(float)

        back_proj_top_left = z_shift / np.tan([xangles, yangles])
        back_proj_top_left[0, :] = back_proj_top_left[0, :] + back_top_left[0].astype(
            float
        )
        back_proj_top_left[1, :] = -back_proj_top_left[1, :] + back_top_left[1].astype(
            float
        )

        back_proj_top_right = z_shift / np.tan([xangles, yangles])
        back_proj_top_right[0, :] = back_proj_top_right[0, :] + back_top_right[
            0
        ].astype(float)
        back_proj_top_right[1, :] = -back_proj_top_right[1, :] + back_top_right[
            1
        ].astype(float)

    # Initialising the module dictionaries that keep the shading fraction and illumination
    shading_dict = {}
    illumination_dict = {}
    shading_dict2 = {}
    back_shading_dict = {}
    back_illumination_dict = {}
    back_shading_dict2 = {}
    front_row_shading_dict = {}
    front_row_illumination_dict = {}
    front_row_shading_dict2 = {}
    back_row_shading_dict = {}
    back_row_illumination_dict = {}
    back_row_shading_dict2 = {}

    # Making the arrays of zeros that go with each module in the dictionary. Each column represents a cell and each row a point in time
    # dict2 is the fractional shading of each module
    # All the cells are set to zero aka fully sunny
    for i in range(pv_array.columns):
        if (
            pv_array.Module.type == "Bifacial"
        ):  # Generates a dictionary for the front and seperate dictionary for the back
            shading_dict2[f"Module {i+1}"] = np.zeros((len(direct)))
            shading_dict[f"Module {i+1}"] = np.zeros(
                (len(direct), (pv_array.Module.cells_long * pv_array.Module.cells_wide))
            )
            back_shading_dict2[f"Module {i+1}"] = np.zeros((len(direct)))
            back_shading_dict[f"Module {i+1}"] = np.zeros(
                (len(direct), (pv_array.Module.cells_long * pv_array.Module.cells_wide))
            )
            front_row_shading_dict2[f"Module {i+1}"] = np.zeros((len(direct)))
            front_row_shading_dict[f"Module {i+1}"] = np.zeros(
                (len(direct), (pv_array.Module.cells_long * pv_array.Module.cells_wide))
            )
            back_row_shading_dict2[f"Module {i+1}"] = np.zeros((len(direct)))
            back_row_shading_dict[f"Module {i+1}"] = np.zeros(
                (len(direct), (pv_array.Module.cells_long * pv_array.Module.cells_wide))
            )
        else:
            shading_dict2[f"Module {i+1}"] = np.zeros((len(direct)))
            shading_dict[f"Module {i+1}"] = np.zeros(
                (len(direct), pv_array.Module.cells_long * pv_array.Module.cells_wide)
            )
            front_row_shading_dict2[f"Module {i+1}"] = np.zeros((len(direct)))
            front_row_shading_dict[f"Module {i+1}"] = np.zeros(
                (len(direct), (pv_array.Module.cells_long * pv_array.Module.cells_wide))
            )

    # Big loop that runs through each point of time, and where the sun hits at each point and calculates the shading
    for index, val in tqdm(
        enumerate(direct), desc="Calculating Shading Profiles", total=direct.shape[0]
    ):  # index is keeping track of what shadow number we are on
        if val == 1:  # if value is 1 then the sun strikes from the front side
            if (
                pv_array.Module.type == "Bifacial"
            ):  # If sun hits the front then the back side of the bifacial module will be fully shaded
                for i in range(pv_array.columns):
                    back_shading_dict[f"Module {i+1}"][
                        index
                    ] = 1  # aka set all modules to shaded
                    back_shading_dict2[f"Module {i+1}"][
                        index
                    ] = 1  # aka set all modules to shaded
                    back_row_shading_dict[f"Module {i+1}"][index] = 1
                    back_row_shading_dict2[f"Module {i+1}"][index] = 1
            if (
                (
                    0 < proj_top_left[0, index] < row_length
                    and 0
                    < proj_top_left[1, index]
                    < pv_array.Module.cells_long_in_width
                )
                or (
                    0 < proj_top_right[0, index] < row_length
                    and 0
                    < proj_top_right[1, index]
                    < pv_array.Module.cells_long_in_width
                )
                or (
                    0 < proj_bottom_left[0, index] < row_length
                    and 0
                    < proj_bottom_left[1, index]
                    < pv_array.Module.cells_long_in_width
                )
                or (
                    0 < proj_bottom_right[0, index] < row_length
                    and 0
                    < proj_bottom_right[1, index]
                    < pv_array.Module.cells_long_in_width
                )
            ):
                # First if statement finds out if the shadow actually strikes a panel

                # Finds the x and y coordinates of the shadow's left (nearside) and right (farside) corners
                cells_from_nearside = np.array(
                    [proj_top_left[0, index], proj_top_left[1, index]]
                )
                cells_from_farside = np.array(
                    [proj_top_right[0, index], proj_top_right[1, index]]
                )

                # Finds the module that the shadow strikes on
                modules_from_nearside = np.floor(
                    cells_from_nearside[0]
                    / (
                        pv_array.Module.cells_wide
                        + (pv_array.xSpacing * 1000 / pv_array.Module.cell.xlength)
                    )
                )

                modules_from_farside = np.floor(
                    cells_from_farside[0]
                    / (
                        pv_array.Module.cells_wide
                        + (pv_array.xSpacing * 1000 / pv_array.Module.cell.xlength)
                    )
                )

                # nearside_module_intersection and farside_module_intersection
                # give where the shadow strikes on the individual module
                nearside_module_intersection = [
                    cells_from_nearside[0]
                    - modules_from_nearside
                    * (
                        pv_array.Module.cells_wide
                        + (pv_array.xSpacing * 1000 / pv_array.Module.cell.xlength)
                    ),
                    cells_from_nearside[1] / pv_array.Module.cell_ratio,
                ]

                farside_module_intersection = [
                    cells_from_farside[0]
                    - modules_from_farside
                    * (
                        pv_array.Module.cells_wide
                        + (pv_array.xSpacing * 1000 / pv_array.Module.cell.xlength)
                    ),
                    cells_from_farside[1] / pv_array.Module.cell_ratio,
                ]

                # This is essentially tests if the shadow falls between the modules see if the shadow falls between the module.
                # If it does we shift the modules_from variable shift it up by half so i never equals it
                modules_from_nearside = (
                    modules_from_nearside + 0.5
                    if nearside_module_intersection[0] > pv_array.Module.cells_wide
                    else modules_from_nearside
                )
                modules_from_farside = (
                    modules_from_farside + 0.5
                    if farside_module_intersection[0] > pv_array.Module.cells_wide
                    else modules_from_farside
                )
                # Essentially:
                #  module column < modules_from_nearside  -->  identical shading profiles (sunny)
                #  module column = modules_from_nearside  -->  calculate distinct shading profile
                #  modules_from_nearside < module column < modules_from_farside  -->  calculate 1 shading profile then repeat for the rest
                #  module column = modules_from_farside -->  distinct shading profile
                #  module column > modules_from_farside -->  indentical shading profiles (sunny)
                for i in range(pv_array.columns):
                    if i < modules_from_nearside:
                        # set to fully sunny aka do no work and continue the loop. All cells are set to zero shading
                        continue
                    if i == modules_from_nearside:
                        cells = np.zeros(
                            (pv_array.Module.cells_wide, pv_array.Module.cells_long)
                        )  # initialising cells for the ith module. The shading is just added on top

                        # All the cells in the bottom right corner are set to fully shaded
                        cells[
                            int(np.ceil(nearside_module_intersection[0])) :,
                            : int(np.floor(nearside_module_intersection[1])),
                        ] = 1

                        # All the cells in the partially shaded row and column are set to the fraction that is shaded
                        cells[
                            int(np.floor(nearside_module_intersection[0])),
                            : int(np.floor(nearside_module_intersection[1])),
                        ] = nearside_module_intersection[0] - np.floor(
                            nearside_module_intersection[0]
                        )
                        cells[
                            int(np.ceil(nearside_module_intersection[0])) :,
                            int(np.floor(nearside_module_intersection[1])),
                        ] = nearside_module_intersection[1] - np.floor(
                            nearside_module_intersection[1]
                        )
                        # The single cell that is at the corner of the shadow is set to the fraction that is shaded
                        cells[
                            int(np.floor(nearside_module_intersection[0])),
                            int(np.floor(nearside_module_intersection[1])),
                        ] = (
                            nearside_module_intersection[1]
                            - np.floor(nearside_module_intersection[1])
                        ) * (
                            nearside_module_intersection[0]
                            - np.floor(nearside_module_intersection[0])
                        )
                        cells = cells.transpose()
                        shading_dict[f"Module {i+1}"][index] = cells.flatten()
                        shading_dict2[f"Module {i+1}"][index] = np.mean(cells.flatten())
                    if modules_from_nearside < i < modules_from_farside:
                        if np.floor(modules_from_nearside) + 1 == i or i == 0:
                            # This runs for all the modules in the middle that are partially cut off horizontally
                            cells = np.zeros(
                                (pv_array.Module.cells_wide, pv_array.Module.cells_long)
                            )  # initialising cells for the ith module.

                            # Cells below the edge of the shadow are all set to zero
                            cells[
                                :, : int(np.floor(nearside_module_intersection[1]))
                            ] = 1

                            # Cells that are in the single partially shaded row are set to the shaded fraction
                            cells[:, int(np.floor(nearside_module_intersection[1]))] = (
                                nearside_module_intersection[1]
                                - np.floor(nearside_module_intersection[1])
                            )
                            cells = (
                                cells.transpose()
                            )  # We only need to do this calculation once and then keep using the result because the modules are identical

                        shading_dict[f"Module {i+1}"][index] = cells.flatten()
                        shading_dict2[f"Module {i+1}"][index] = np.mean(cells.flatten())
                        # calculate shading for one then add for the rest
                    if i == modules_from_farside:
                        # Initialising cells
                        cells = np.zeros(
                            (pv_array.Module.cells_wide, pv_array.Module.cells_long)
                        )

                        # Sets cells in the bottom right corner of shadow to fully shaded
                        cells[
                            : (int(np.floor(farside_module_intersection[0]))),
                            : int(np.floor(farside_module_intersection[1])),
                        ] = 1

                        # Sets cells in partially shaded row and column to their respective shadings
                        cells[
                            int(np.floor(farside_module_intersection[0])),
                            : int(np.floor(farside_module_intersection[1])),
                        ] = farside_module_intersection[0] - np.floor(
                            farside_module_intersection[0]
                        )

                        cells[
                            : (int(np.floor(farside_module_intersection[0]))),
                            int(np.floor(farside_module_intersection[1])),
                        ] = farside_module_intersection[1] - np.floor(
                            farside_module_intersection[1]
                        )

                        # Sets cell at corner of the shadow to the fractional shading
                        cells[
                            int(np.floor(farside_module_intersection[0])),
                            int(np.floor(farside_module_intersection[1])),
                        ] = (
                            farside_module_intersection[1]
                            - np.floor(farside_module_intersection[1])
                        ) * (
                            farside_module_intersection[0]
                            - np.floor(farside_module_intersection[0])
                        )
                        cells = cells.transpose()
                        shading_dict[f"Module {i+1}"][index] = cells.flatten()
                        shading_dict2[f"Module {i+1}"][index] = np.mean(cells.flatten())
                    if i > modules_from_farside:
                        # All the modules after farside will be fully sunny aka set to zero so we can just break loop
                        break
        if val == 0:
            if pv_array.Module.type == "Bifacial":
                if (
                    (
                        0 < back_proj_top_left[0, index] < row_length
                        and 0
                        < back_proj_top_left[1, index]
                        < pv_array.Module.cells_long_in_width
                    )
                    or (
                        0 < back_proj_top_right[0, index] < row_length
                        and 0
                        < back_proj_top_right[1, index]
                        < pv_array.Module.cells_long_in_width
                    )
                    or (
                        0 < back_proj_bottom_left[0, index] < row_length
                        and 0
                        < back_proj_bottom_left[1, index]
                        < pv_array.Module.cells_long_in_width
                    )
                    or (
                        0 < back_proj_bottom_right[0, index] < row_length
                        and 0
                        < back_proj_bottom_right[1, index]
                        < pv_array.Module.cells_long_in_width
                    )
                ):

                    # The coordinates of the top corners of the shadows (the ones that will shade the modules)
                    cells_from_nearside = np.array(
                        [back_proj_top_left[0, index], back_proj_top_left[1, index]]
                    )
                    cells_from_farside = np.array(
                        [back_proj_top_right[0, index], back_proj_top_right[1, index]]
                    )

                    # Finds the module that the shadow strikes on
                    modules_from_nearside = np.floor(
                        cells_from_nearside[0]
                        / (
                            pv_array.Module.cells_wide
                            + (pv_array.xSpacing * 1000 / pv_array.Module.cell.xlength)
                        )
                    )

                    modules_from_farside = np.floor(
                        cells_from_farside[0]
                        / (
                            pv_array.Module.cells_wide
                            + (pv_array.xSpacing * 1000 / pv_array.Module.cell.xlength)
                        )
                    )

                    # Ensures that the modules_from_ variables are within not higher than the number of columns
                    # or lower than -1. Should prevent the if statement below from breaking
                    # modules_from_nearside = min(modules_from_nearside, pv_array.columns + 1)
                    # modules_from_farside = min(modules_from_farside, pv_array.columns + 1)
                    # modules_from_nearside = max(-1, modules_from_nearside)
                    # modules_from_farside = max(-1, modules_from_farside)

                    # nearside_module_intersection and farside_module_intersection give where the shadow strikes in relative to a single module
                    nearside_module_intersection = [
                        cells_from_nearside[0]
                        - modules_from_nearside
                        * (
                            pv_array.Module.cells_wide
                            + (pv_array.xSpacing * 1000 / pv_array.Module.cell.xlength)
                        ),
                        cells_from_nearside[1] / pv_array.Module.cell_ratio,
                    ]

                    farside_module_intersection = [
                        cells_from_farside[0]
                        - modules_from_farside
                        * (
                            pv_array.Module.cells_wide
                            + (pv_array.xSpacing * 1000 / pv_array.Module.cell.xlength)
                        ),
                        cells_from_farside[1] / pv_array.Module.cell_ratio,
                    ]

                    # Ensures that the intersection variables are within the limits of the module cell (as opposed to in the gap)
                    # nearside_module_intersection[0] = min(nearside_module_intersection[0], pv_array.Module.cells_wide-1)
                    # farside_module_intersection[0] = min(farside_module_intersection[0], pv_array.Module.cells_wide-1)

                    # This is essentially because the shadow falls between the modules. We shift it up by half so i never equals it
                    modules_from_nearside = (
                        modules_from_nearside + 0.5
                        if nearside_module_intersection[0] >= pv_array.Module.cells_wide
                        else modules_from_nearside
                    )
                    modules_from_farside = (
                        modules_from_farside + 0.5
                        if farside_module_intersection[0] >= pv_array.Module.cells_wide
                        else modules_from_farside
                    )

                    for i in range(pv_array.columns):
                        if i < modules_from_nearside:
                            # set to fully sunny aka do no work and continue the loop. All cells are set to zero shading
                            continue
                        if i == modules_from_nearside:
                            cells = np.zeros(
                                (pv_array.Module.cells_wide, pv_array.Module.cells_long)
                            )  # initialising cells for the ith module. The shading is just added on top
                            cells[
                                int(np.ceil(nearside_module_intersection[0])) :,
                                : int(np.floor(nearside_module_intersection[1])),
                            ] = 1
                            cells[
                                int(np.floor(nearside_module_intersection[0])),
                                : int(np.floor(nearside_module_intersection[1])),
                            ] = nearside_module_intersection[0] - np.floor(
                                nearside_module_intersection[0]
                            )
                            cells[
                                int(np.ceil(nearside_module_intersection[0])) :,
                                int(np.floor(nearside_module_intersection[1])),
                            ] = nearside_module_intersection[1] - np.floor(
                                nearside_module_intersection[1]
                            )
                            cells[
                                int(np.floor(nearside_module_intersection[0])),
                                int(np.floor(nearside_module_intersection[1])),
                            ] = (
                                nearside_module_intersection[1]
                                - np.floor(nearside_module_intersection[1])
                            ) * (
                                nearside_module_intersection[0]
                                - np.floor(nearside_module_intersection[0])
                            )
                            cells = cells.transpose()
                            back_shading_dict[f"Module {i+1}"][index] = cells.flatten()
                            back_shading_dict2[f"Module {i+1}"][index] = np.mean(
                                cells.flatten()
                            )
                        if modules_from_nearside < i < modules_from_farside:
                            if np.floor(modules_from_nearside) + 1 == i or i == 0:
                                cells = np.zeros(
                                    (
                                        pv_array.Module.cells_wide,
                                        pv_array.Module.cells_long,
                                    )
                                )  # initialising cells for the ith module.
                                cells[
                                    :, : int(np.floor(nearside_module_intersection[1]))
                                ] = 1
                                cells[
                                    :, int(np.floor(nearside_module_intersection[1]))
                                ] = nearside_module_intersection[1] - np.floor(
                                    nearside_module_intersection[1]
                                )
                                cells = (
                                    cells.transpose()
                                )  # We only need to do this calculation once and then keep using the result
                            back_shading_dict[f"Module {i+1}"][index] = cells.flatten()
                            back_shading_dict2[f"Module {i+1}"][index] = np.mean(
                                cells.flatten()
                            )
                            # calculate shading for one then add for the rest
                        if i == modules_from_farside:
                            cells = np.zeros(
                                (pv_array.Module.cells_wide, pv_array.Module.cells_long)
                            )
                            cells[
                                : (int(np.floor(farside_module_intersection[0]))),
                                : int(np.floor(farside_module_intersection[1])),
                            ] = 1
                            cells[
                                int(np.floor(farside_module_intersection[0])),
                                : int(np.floor(farside_module_intersection[1])),
                            ] = farside_module_intersection[0] - np.floor(
                                farside_module_intersection[0]
                            )
                            cells[
                                : (int(np.floor(farside_module_intersection[0]))),
                                int(np.floor(farside_module_intersection[1])),
                            ] = farside_module_intersection[1] - np.floor(
                                farside_module_intersection[1]
                            )
                            cells[
                                int(np.floor(farside_module_intersection[0])),
                                int(np.floor(farside_module_intersection[1])),
                            ] = (
                                farside_module_intersection[1]
                                - np.floor(farside_module_intersection[1])
                            ) * (
                                farside_module_intersection[0]
                                - np.floor(farside_module_intersection[0])
                            )
                            cells = cells.transpose()
                            back_shading_dict[f"Module {i+1}"][index] = cells.flatten()
                            back_shading_dict2[f"Module {i+1}"][index] = np.mean(
                                cells.flatten()
                            )
                        if i > modules_from_farside:
                            break
            for i in range(
                pv_array.columns
            ):  # Set all the front side modules to be fully shaded aka 1
                shading_dict[f"Module {i+1}"][
                    index
                ] = 1  # aka set all modules to shaded
                shading_dict2[f"Module {i+1}"][
                    index
                ] = 1  # aka set all modules to shaded
                front_row_shading_dict[f"Module {i+1}"][index] = 1
                front_row_shading_dict2[f"Module {i+1}"][index] = 1
        if val == -1:  # Sun below the horizon
            for i in range(pv_array.columns):
                # aka set all modules to shaded
                shading_dict[f"Module {i+1}"][index] = 1
                shading_dict2[f"Module {i+1}"][index] = 1
                front_row_shading_dict[f"Module {i+1}"][index] = 1
                front_row_shading_dict2[f"Module {i+1}"][index] = 1
                if pv_array.Module.type == "Bifacial":
                    back_shading_dict[f"Module {i+1}"][index] = 1
                    back_shading_dict2[f"Module {i+1}"][index] = 1
                    back_row_shading_dict[f"Module {i+1}"][index] = 1
                    back_row_shading_dict2[f"Module {i+1}"][index] = 1

    # Calculating the actual illumination in watts per m^2 after being given solar irradiance
    for i in range(pv_array.columns):
        # Doing this transposed version to make the dimensions and broadcasting work
        # Finding fraction that is illuminated multiplying by irradiance and the cosine
        # of angle of incidence (strike_angles) then adding the diffuse component (taken from some empirical formula. Best to do something smarter...)
        transposed_illumination: np.ndarray = (
            np.transpose((1 - shading_dict[f"Module {i+1}"]))
            * np.cos(strike_angles)
            * solar_irradiance
            + solar_irradiance * 0.1 * (1 - np.cos(pv_array.tilt)) / 2
        )
        illumination_dict[f"Module {i+1}"] = transposed_illumination.transpose()

        transposed_illumination: np.ndarray = (
            np.transpose((1 - front_row_shading_dict[f"Module {i+1}"]))
            * np.cos(strike_angles)
            * solar_irradiance
            + solar_irradiance * 0.1 * (1 - np.cos(pv_array.tilt)) / 2
        )
        front_row_illumination_dict[f"Module {i+1}"] = (
            transposed_illumination.transpose()
        )

        if pv_array.Module.type == "Bifacial":
            transposed_illumination: np.ndarray = (
                np.transpose((1 - back_shading_dict[f"Module {i+1}"]))
                * np.cos(np.pi - strike_angles)
                * solar_irradiance
                + solar_irradiance * 0.1 * (1 - np.cos(np.pi - pv_array.tilt)) / 2
            )
            back_illumination_dict[f"Module {i+1}"] = (
                transposed_illumination.transpose()
            )
            transposed_illumination: np.ndarray = (
                np.transpose((1 - back_row_shading_dict[f"Module {i+1}"]))
                * np.cos(np.pi - strike_angles)
                * solar_irradiance
                + solar_irradiance * 0.1 * (1 - np.cos(np.pi - pv_array.tilt)) / 2
            )
            back_row_illumination_dict[f"Module {i+1}"] = (
                transposed_illumination.transpose()
            )

    # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
    #  Below is plotting for visualisation and verification purposes only #
    # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
    if plot_graphs == True:

        # Graph 1. Sun angles
        SunRay = np.vstack([SunRay, normal.flatten()])
        colors = np.linspace(0, 1, direct.shape[0] + 1)
        fig = plt.figure(figsize=(8, 6))
        ax = fig.add_subplot(111, projection="3d")
        ax.scatter(SunRay[:, 0], SunRay[:, 1], SunRay[:, 2], c=colors, cmap="viridis")
        ax.quiver(
            np.zeros_like(SunRay[:, 0]),
            np.zeros_like(SunRay[:, 0]),
            np.zeros_like(SunRay[:, 0]),
            SunRay[:, 0],
            SunRay[:, 1],
            SunRay[:, 2],
            color="black",
            alpha=0.6,
            arrow_length_ratio=0.1,
        )
        x_surface = np.linspace(-2, 2, 12)
        y_surface = np.linspace(-2, 2, 12)
        x_surface, y_surface = np.meshgrid(x_surface, y_surface)
        ax.plot_surface(x_surface, y_surface, x_surface * 0, alpha=0.2)
        plt.title("Sun Angles with normal of panel")
        plt.show()

        # Graph 2. Plotting the middle of the shadow on the module row plane
        colors = np.linspace(0, 1, direct.shape[0])
        x = proj_bottom_left[0] + row_length / 2
        y = proj_bottom_left[1] + pv_array.Module.cells_long_in_width / 2
        scatter = plt.scatter(x, y, alpha=0.5, c=colors, s=20)
        plt.colorbar(scatter, label="Gradient")
        plt.title("Projection of middle of shadow onto module plane")
        plt.show()

        # Graph 3. Plotting the shadows on the modules
        ones = 0
        for i, v in enumerate(direct):
            if v == 1:
                ones += 1
        front_shadow_palette = sns.color_palette("flare", ones)
        back_shadow_palette = sns.color_palette(
            "crest", n_colors=direct.shape[0] - np.count_nonzero(direct)
        )
        fsp = 0
        bsp = 0
        for index, value in enumerate(direct):
            if value == 1:
                x = np.array(
                    [
                        proj_bottom_left[0, index],
                        proj_bottom_right[0, index],
                        proj_top_right[0, index],
                        proj_top_left[0, index],
                    ]
                )
                y = np.array(
                    [
                        proj_bottom_left[1, index],
                        proj_bottom_right[1, index],
                        proj_top_right[1, index],
                        proj_top_left[1, index],
                    ]
                )
                plt.fill(x, y, alpha=0.5, fc=front_shadow_palette[fsp], ec="black")
                fsp += 1

            if value == 0 and pv_array.Module.type == "Bifacial":
                x = np.array(
                    [
                        back_proj_bottom_left[0, index],
                        back_proj_bottom_right[0, index],
                        back_proj_top_right[0, index],
                        back_proj_top_left[0, index],
                    ]
                )
                y = np.array(
                    [
                        back_proj_bottom_left[1, index],
                        back_proj_bottom_right[1, index],
                        back_proj_top_right[1, index],
                        back_proj_top_left[1, index],
                    ]
                )
                plt.fill(x, y, alpha=0.5, fc=back_shadow_palette[bsp], ec="black")
                bsp += 1

        for i in range(pv_array.columns):
            x = np.array(
                [
                    i
                    * (
                        pv_array.Module.cells_wide
                        + pv_array.xSpacing * 1000 / (pv_array.Module.cell.xlength)
                    ),
                    pv_array.Module.cells_wide
                    + i
                    * (
                        pv_array.Module.cells_wide
                        + pv_array.xSpacing * 1000 / (pv_array.Module.cell.xlength)
                    ),
                    pv_array.Module.cells_wide
                    + i
                    * (
                        pv_array.Module.cells_wide
                        + pv_array.xSpacing * 1000 / (pv_array.Module.cell.xlength)
                    ),
                    i
                    * (
                        pv_array.Module.cells_wide
                        + pv_array.xSpacing * 1000 / (pv_array.Module.cell.xlength)
                    ),
                ]
            )
            y = np.array(
                [
                    0,
                    0,
                    pv_array.Module.cells_long_in_width,
                    pv_array.Module.cells_long_in_width,
                ]
            )
            plt.fill(x, y, alpha=0.5, fc="none", ec="black")
        plt.xlim(-pv_array.Module.cells_wide, row_length + pv_array.Module.cells_wide)
        plt.ylim(-11 * pv_array.Module.cells_wide, pv_array.Module.cells_wide)
        plt.gca().set_aspect("equal")
        plt.title("Modules with projected shadows")

        plt.show()

    # Creating the dictionary of dictionaries
    if pv_array.Module.type == "Bifacial":
        dicts = {
            "illumination": illumination_dict,
            "shading": shading_dict,
            "back illumination": back_illumination_dict,
            "back shading": back_shading_dict,
            "front row illumination": front_row_illumination_dict,
            "back row illumination": back_row_illumination_dict,
        }
    else:
        dicts = {
            "illumination": illumination_dict,
            "shading": shading_dict,
            "front row illumination": front_row_illumination_dict,
        }
    return dicts


def calculate_crop_shading(
    pv_array: Array,
    solar_irradiance: pd.DataFrame,
    Time: tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray],
    resolution: float = 1,
    spectral_calculation: bool = False,
    ground_coordinates: np.ndarray = np.array([]),
    plot_graph: bool = False,
) -> tuple[np.ndarray, np.ndarray]:
    """
    Somewhat untested :(

    Calculates what ground shaded at each time point.

    Parameters
    ----------
    - pv_array: Array
        Holds the attributes of the pv array to be modelled
    - solar_irradiance: np.ndarray
        Solar irridiance for each point in time
    - Time: tuple[np.ndarray, np.ndarray, np.ndarray]
        Tuple that holds all the years, months, days and hours of each individaul time point
        to be calculated. Most importantly, for each hour the respective year, month and day
        must be added to set of arrays. The length of each array in the Time tuple must be the same!
    - resolution: float
        Gives the resolution in meters of the shading. Shading is returned per resolution^2 meter^2 squares
    - spectral calculations: bool, optional
        Boolean that says when to calculate the spoectral analysis. Currently means nothing
    - ground_coordinates: ndarray, optional
        Gives corners of the area of ground you would like to analyse. If not given function
        just calculates area around panels. Should have a shape (4, 2).
        Coordinates are given with respect to the bottom left corner of the bottom left module of the pv array and
        given in anticlockwise direction from there.
     - plot_graph: bool, optional
        Boolean that plots each shadow in the ground

    Returns
    ----------
    - crop_shading: ndarray
        3D array of shaded and non shaded squares of land. 0 axis is time, 1 axis is x axis and 2 axis is y axis
    - crop_illumination: ndarray
        3D array of irradiance of squares of land. 0 axis is time, 1 axis is x axis and 2 axis is y axis
    """
    sunEl_rad, azimuth_rad = findAngles(
        pv_array.longitude,
        pv_array.latitude,
        Time[0],
        Time[1],
        Time[2],
        Time[3],
        pv_array.TimeZone,
    )

    # Defing vectors the point along the edge of the module
    vector_width = np.array(
        [
            [np.sin(pv_array.orientation - np.pi / 2)],
            [np.cos(pv_array.orientation - np.pi / 2)],
            [0],
        ]
    )
    vector_length = np.array(
        [
            [np.sin(np.pi + pv_array.orientation) * np.cos(pv_array.tilt)],
            [np.cos(np.pi + pv_array.orientation) * np.cos(pv_array.tilt)],
            [np.sin(pv_array.tilt)],
        ]
    )

    # Defining the normal (aka front) of the panel
    normal = np.cross(vector_width.flatten(), vector_length.flatten())

    # Defining Sun Ray vector
    SunRay = np.array(
        [
            np.cos(sunEl_rad)
            * np.cos(
                np.pi / 2 - (azimuth_rad)
            ),  # Why np.pi/2 - azimuth  you ask? North is the y axis and azimuth is measured clockwise
            np.cos(sunEl_rad) * np.sin(np.pi / 2 - (azimuth_rad)),
            np.sin(sunEl_rad),
        ]
    )
    SunRay = SunRay.transpose()

    # Defining angles that sun hits the panel
    strike_angles = np.arccos(np.dot(SunRay, normal))
    direct = np.ones_like(strike_angles)
    direct[strike_angles >= np.pi / 2] = 0  # equals 0 if light hits backside
    direct[sunEl_rad <= 0] = -1

    # Redefining SunRay with respect to the panels
    SunRay = np.array(
        [
            np.cos(sunEl_rad)
            * np.cos(np.pi / 2 - (azimuth_rad - (pv_array.orientation - np.pi))),
            np.cos(sunEl_rad)
            * np.sin(np.pi / 2 - (azimuth_rad - (pv_array.orientation - np.pi))),
            np.sin(sunEl_rad),
        ]
    )
    SunRay = SunRay.transpose()

    # Defining length of one row in units of resolution
    row_length = (
        pv_array.columns * pv_array.Module.width
        + (pv_array.columns - 1) * pv_array.xSpacing
    )

    # Defining corners of pv row shading. Assuming column spacing is negligible
    bottom_right_row_corner = (
        np.array([pv_array.Module.length + row_length, 0, pv_array.Module.height])
        / resolution
    )
    bottom_left_row_corner = np.array([0, 0, pv_array.Module.height]) / resolution
    top_right_row_corner = (
        np.array(
            [
                pv_array.Module.length + row_length,
                pv_array.Module.length * np.cos(pv_array.tilt),
                pv_array.Module.height + pv_array.Module.length * np.sin(pv_array.tilt),
            ]
        )
        / resolution
    )
    top_left_row_corner = (
        np.array(
            [
                0,
                pv_array.Module.length * np.cos(pv_array.tilt),
                pv_array.Module.height + pv_array.Module.length * np.sin(pv_array.tilt),
            ]
        )
        / resolution
    )

    # Defing a projection vector
    projector = np.array([SunRay[:, 0], SunRay[:, 1]]) / -SunRay[:, 2]

    # Defining the ground
    # Scaling the ground to be in units of resolution
    if not np.any(
        ground_coordinates,
    ):
        ground_coordinates = np.floor(
            np.array(
                [
                    [-1, -1],
                    [row_length + 1, -1],
                    [row_length + 1, pv_array.ySpacing * pv_array.rows + 1],
                    [-1, pv_array.ySpacing * pv_array.rows + 1],
                ]
            )
            / resolution
        )
    elif ground_coordinates.shape == (2, 4):
        ground_coordinates = ground_coordinates.transpose()
        ground_coordinates: np.ndarray = ground_coordinates / resolution
    elif ground_coordinates.shape == (4, 2):
        ground_coordinates: np.ndarray = ground_coordinates / resolution
    else:
        print("None or invalid ground coordinates input.")
        return None

    ground = Polygon(
        (
            [0, 0],
            [np.max(ground_coordinates[:, 0]) - np.min(ground_coordinates[:, 0]), 0],
            [
                np.max(ground_coordinates[:, 0]) - np.min(ground_coordinates[:, 0]),
                np.max(ground_coordinates[:, 1]) - np.min(ground_coordinates[:, 1]),
            ],
            [
                0,
                np.max(ground_coordinates[:, 1]) - np.min(ground_coordinates[:, 1]),
            ],
        )
    )

    # Defining the crop squares shading.
    # This is an 3D array with time as first dim, x axis as secodn and y axis as third
    # We measure in units of resolution
    crop_squares_shading = np.zeros(
        (
            direct.shape[0],
            int(
                np.floor(
                    np.max(ground_coordinates[:, 0]) - np.min(ground_coordinates[:, 0])
                )
            ),
            int(
                np.floor(
                    np.max(ground_coordinates[:, 1]) - np.min(ground_coordinates[:, 1])
                )
            ),
        )
    )

    # Doing calculations for each shadow cast to find fractional area
    colours = sns.color_palette("flare", int(direct.shape[0]))
    count = 0
    for index, val in tqdm(
        enumerate(direct), desc="Calculating Ground Shading", total=direct.shape[0]
    ):
        if val != -1:
            for row in range(pv_array.rows):
                # Creating shadow of row on the ground
                proj_bottom_right: np.ndarray = (
                    np.array(
                        [
                            bottom_right_row_corner[0] - ground_coordinates[0, 0],
                            bottom_right_row_corner[1]
                            - ground_coordinates[0, 1]
                            + row * pv_array.ySpacing / resolution,
                        ]
                    )
                    + projector[:, index] * bottom_right_row_corner[2]
                )

                proj_bottom_left: np.ndarray = (
                    np.array(
                        [
                            bottom_left_row_corner[0] - ground_coordinates[0, 0],
                            bottom_left_row_corner[1]
                            - ground_coordinates[0, 1]
                            + row * pv_array.ySpacing / resolution,
                        ]
                    )
                    + projector[:, index] * bottom_left_row_corner[2]
                )

                proj_top_left: np.ndarray = (
                    np.array(
                        [
                            top_left_row_corner[0] - ground_coordinates[0, 0],
                            top_left_row_corner[1]
                            + row * pv_array.ySpacing / resolution
                            - ground_coordinates[0, 1],
                        ]
                    )
                    + projector[:, index] * top_left_row_corner[2]
                )

                proj_top_right: np.ndarray = (
                    np.array(
                        [
                            top_right_row_corner[0] - ground_coordinates[0, 0],
                            top_right_row_corner[1]
                            + row * pv_array.ySpacing / resolution
                            - ground_coordinates[0, 1],
                        ]
                    )
                    + projector[:, index] * top_right_row_corner[2]
                )

                # 1. Find whole numbers between top right and bottom right corners y coordinates.
                y_coords = np.array(
                    range(
                        int(
                            np.floor(np.min([proj_bottom_right[1], proj_top_right[1]]))
                        ),
                        int(np.floor(np.max([proj_top_right[1], proj_bottom_right[1]])))
                        + 1,
                    )
                )

                # 2. Use np.interp to interpolate x coords to find the x coords of the whole numbers
                if proj_bottom_left[1] < proj_top_left[1]:
                    x_coords = np.interp(
                        y_coords + 0.5,
                        [proj_bottom_left[1], proj_top_left[1]],
                        [proj_bottom_left[0], proj_top_left[0]],
                    )
                elif proj_bottom_left[1] > proj_top_left[1]:
                    x_coords = np.interp(
                        y_coords + 0.5,
                        [proj_top_left[1], proj_bottom_left[1]],
                        [proj_top_left[0], proj_bottom_left[0]],
                    )
                else:
                    x_coords = proj_bottom_left[0]

                # 3. Rounding down to whole numbers aka matching resolution
                y_coords = np.floor(y_coords)
                x_coords = np.floor(x_coords)

                # Defining the end x coordinates of the shadow
                end_x_coords = x_coords + np.floor(row_length / resolution)

                # 4. Minimising to be in range of ground coordinates

                x_coords[x_coords >= crop_squares_shading.shape[1] - 1] = (
                    crop_squares_shading.shape[1] - 1
                )
                x_coords[x_coords <= 0] = 0
                end_x_coords[end_x_coords >= crop_squares_shading.shape[1] - 1] = (
                    crop_squares_shading.shape[1] - 1
                )
                end_x_coords[end_x_coords <= 0] = 0

                x_coords = x_coords[
                    (0 <= y_coords) & (y_coords <= crop_squares_shading.shape[2] - 1)
                ]
                end_x_coords = end_x_coords[
                    (0 <= y_coords) & (y_coords <= crop_squares_shading.shape[2] - 1)
                ]
                y_coords = y_coords[
                    (0 <= y_coords) & (y_coords <= crop_squares_shading.shape[2] - 1)
                ]
                # 5. Use for each y coord set all indices between x coord and end x coord to shaded
                #   to shaded in crop_squares_shading
                for i, ycoord in enumerate(y_coords.astype("int")):
                    crop_squares_shading[
                        index,
                        x_coords.astype("int")[i] : end_x_coords.astype("int")[i] + 1,
                        ycoord,
                    ] = 1

                # If plotting, making shadow and adding to graph
                if plot_graph == True:
                    shadow = Polygon(
                        (
                            proj_bottom_left,
                            proj_bottom_right,
                            proj_top_right,
                            proj_top_left,
                        )
                    )
                    x, y = shadow.exterior.xy
                    plt.fill(x, y, alpha=0.5, fc=colours[count], ec="black")
            count += 1
    if plot_graph == True:
        x, y = ground.exterior.xy
        plt.fill(x, y, alpha=0.5, fc="None", ec="black")
        plt.gca().set_aspect("equal")
        plt.title("Shadows in the daytime!")

    solar_irradiance = solar_irradiance.reshape(crop_squares_shading.shape[0], 1, 1)
    illumination = (
        1 - crop_squares_shading
    ) * solar_irradiance + 0.1 * solar_irradiance
    # Diffuse light is given by a first order approximation of 10% of direct illumination

    return crop_squares_shading, illumination