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

from Utilities.Processes import calculate_no_of_panels, calculate_required_system_size

logger = logging.getLogger(__name__)


def get_location(c):
    location = pvlib.location.Location(
        latitude=c["environment"]["location"]["latitude"],
        longitude=c["environment"]["location"]["longitude"],
        tz=c["simulation_date_time"]["tz"],
    )

    return location


def define_grid_layout(c, panel_tilt):
    no_of_panels = calculate_no_of_panels(
        c["array"]["system_size"], c["panel"]["peak_power"]
    )

    # calculate pitch
    pitch = c["array"]["spacing"] + c["panel"]["dimensions"]["thickness"]
    # calculate minimum pitch if we don't want panel overlap at all
    min_pitch = c["panel"]["dimensions"]["length"] * math.cos(
        panel_tilt / 180 * math.pi
    )
    if pitch < min_pitch:
        logger.warning(
            f"Spacing is less than minimum pitch. Setting spacing to {min_pitch}."
        )
        pitch = min_pitch
    logger.info(f"Pitch between panels: {pitch}m")

    # get maximum number of panels based on spacing and dimensions
    max__panels_per_row = np.floor(
        (
            c["environment"]["roof"]["dimensions"]["width"]
            - (2 * c["array"]["edge_setback"] + c["panel"]["dimensions"]["width"])
        )
        / c["panel"]["dimensions"]["width"]
    )
    max_number_of_rows = np.floor(
        (
            c["environment"]["roof"]["dimensions"]["length"]
            - (2 * c["array"]["edge_setback"])
        )
        / pitch
    )
    max_no_of_panels = max__panels_per_row * max_number_of_rows
    logger.info(
        f"Number of panels required: {no_of_panels}, Maximum panels possible: {max_no_of_panels}"
    )
    if no_of_panels > max_no_of_panels:
        no_of_panels = max_no_of_panels
        logger.warning(
            f"Number of panels required exceeds maximum possible. Setting number of panels to {no_of_panels}."
        )
    else:
        logger.info(
            f"Number of panels required is within the maximum possible. Setting number of panels to {no_of_panels}."
        )

    # coordinate of panel determined by bottom left corner
    # x - row wise position, y - column wise position, z - height
    # first panel in row 1 is at (0, 0, 0)
    # nth panel in row 1 is at ((n-1)*panel_width, 0, 0)
    # first panel in nth row is at (0, (n-1)*(panel_thickness + spacing), 0)

    # create matrices for x, y, z coordinates of panels
    x = []
    y = []
    z = []
    counter = 0
    for j in range(int(max_number_of_rows)):
        for i in range(int(max__panels_per_row)):
            if counter < no_of_panels:
                x.append(i * c["panel"]["dimensions"]["width"])
                y.append(j * pitch)
                z.append(0)
                counter += 1
            else:
                break

    coordinates = pd.DataFrame(
        {
            "x": x,
            "y": y,
            "z": z,
        }
    )
    return coordinates, no_of_panels


def get_solar_data(c):
    logger.info(
        f"Getting solar position data for {c['simulation_date_time']['start']} to {c['simulation_date_time']['end']}"
    )
    """
    Function to get solar position from PVLib
    """
    location = get_location(c)

    times = pd.date_range(
        c["simulation_date_time"]["start"],
        c["simulation_date_time"]["end"],
        freq="15min",
        tz=location.tz,
    )

    # Get solar position data using PVLib
    solar_positions = location.get_solarposition(times)
    # filter solar positions to only include times when the sun is above the horizon
    solar_positions = solar_positions[solar_positions["apparent_elevation"] > 0]
    # get datetime range from solar_positions
    day_times = solar_positions.index
    clearsky_data = location.get_clearsky(day_times)

    return solar_positions, clearsky_data


def sanity_check_minimum_pitch(c):
    solar_positions, _ = get_solar_data(c)
    zenith = solar_positions["zenith"].values
    solar_positions["shadow_length"] = (c["panel"]["dimensions"]["length"]) * np.tan(
        zenith / 180 * np.pi
    )
    return solar_positions


def calculate_energy_production_vertical(c):
    c = calculate_required_system_size(c)
    panel_coordinates, no_of_panels = define_grid_layout(c, panel_tilt=90)
    solar_positions, clearsky_data = get_solar_data(c)

    # the first row is always not shaded so exclude
    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

    gcr = np.divide(c["panel"]["dimensions"]["length"], pitch)
    gcr = min(1, gcr)
    logger.info(f"Ground coverage ratio: {gcr}")

    # use pvfactors bifacial modelling package
    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"],
            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=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

    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"]
    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=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

    total_energy = total_hourly_energy.sum()
    logger.info(f"Total energy yield calculated: {total_energy} kWh")

    return total_hourly_energy, no_of_panels


def calculate_energy_production_horizontal(c):
    c["array"]["system_size"] = (
        c["array"]["system_size"] * c["array"]["horizontal_max_capacity"]
    )
    panel_coordinates, no_of_panels = define_grid_layout(c, panel_tilt=0)
    solar_positions, clearsky_data = get_solar_data(c)

    # the first row is always not shaded so exclude
    no_of_rows = np.unique(panel_coordinates["y"]).shape[0]
    no_of_shaded_rows = no_of_rows - 1

    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 = 0
    shading_row_rotation = 0
    axis_tilt = 0
    axis_azimuth = 270  # south facing surface

    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,
    )

    shaded_fraction = pvlib.shading.shaded_fraction1d(
        solar_zenith=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,
    )
    shaded_fraction = shaded_fraction * no_of_shaded_rows / no_of_rows

    poa = pvlib.irradiance.get_total_irradiance(
        surface_tilt=0,
        surface_azimuth=180,
        solar_zenith=projected_solar_zenith,
        solar_azimuth=solar_positions["azimuth"],
        dni=clearsky_data["dni"],
        ghi=clearsky_data["ghi"],
        dhi=clearsky_data["dhi"],
        surface_type="urban",
    )
    poa = poa.dropna(subset=["poa_global"])

    effective_front = poa["poa_global"] * (1 - shaded_fraction)

    system_size = c["panel"]["peak_power"] * no_of_panels

    pdc0 = system_size
    gamma_pdc = c["panel"]["temperature_coefficient"]
    temp_cell = c["panel"]["nominal_operating_cell_temperature"]
    pdc = pvlib.pvsystem.pvwatts_dc(
        pdc0=pdc0,
        gamma_pdc=gamma_pdc,
        temp_cell=temp_cell,
        g_poa_effective=effective_front,
    )

    total_hourly_energy = pdc * 15 / 60 / 1e3  # convert to kWh
    total_energy = total_hourly_energy.sum()
    logger.info(f"Total energy yield calculated: {total_energy} kWh")

    return total_hourly_energy, no_of_panels