"""
Class for panel calculation
"""

from enum import Enum
import numpy as np

from datetime import datetime, timedelta

import matplotlib.pyplot as plt
import matplotlib.dates as mdates

from .enums import TrackingType
from .utils import gas_parameters_altitude
import numpy as np


class SolarPanel(object):
    """
    Photovoltaic solar panel class

    Input Parameters
    ----------
    latitude : float
        latitude (-90 to 90) in degrees
    longitude : float
        longitude, (-180 to 180) in degrees
    min_altitude : float
        minimal flight altitude above sea level in meters
    date_departure : datetime object
        date of departure
    flight_days : datetime object
        number of expected flight days
    required_power : float
        required daily power capacity
    surface_density : float
        surface density of solar panel [kg/m2]
    safety_coefficient : float
        coefficient of excess power
    efficiency : float
        panel efficiency
    tracking_type : enum
        tracking type for solar panel: "no tracking"; "1 axis tracking", "2 axis tracking"

    Output Parameters
    ----------
    weight : float
        weight of panel in kg
    surface : float
        minimal panel surface in m2

    """

    def __init__(
        self,
        efficiency,
        safety_coefficient,
        surface_density,
        latitude,
        longitude,
        min_altitude,
        tracking_type,
        date_departure,
        flight_days,
    ):
        # Input validations for efficiency, required_power, safety_coefficient, and surface_density
        if not (0 < efficiency < 1):
            raise ValueError("Efficiency must be between 0 and 1, exclusive.")

        if not (0 < safety_coefficient <= 10):
            raise ValueError(
                "Safety coefficient must be between 0 and 10, exclusive of 0 and inclusive of 10."
            )
        if surface_density <= 0:
            raise ValueError("Surface density must be a positive value greater than 0.")

        # Input validations for latitude, longitude, and min_altitude
        if not (-90 <= latitude <= 90):
            raise ValueError("Latitude must be between -90 and 90 degrees.")
        if not (-180 <= longitude <= 180):
            raise ValueError("Longitude must be between -180 and 180 degrees.")
        if min_altitude <= 0:
            raise ValueError(
                "Minimum altitude must be a positive value greater than 0."
            )

        # Input validation for tracking type
        if not isinstance(tracking_type, TrackingType):
            raise ValueError(
                "Invalid tracking type. Must be an instance of TrackingType Enum."
            )
        # Initialize attributes with validated inputs
        self.efficiency = efficiency
        self.safety_coefficient = safety_coefficient
        self.surface_density = surface_density
        self.latitude = latitude
        self.longitude = longitude
        self.altitude = min_altitude
        self.tracking_type = tracking_type
        self.date_departure = date_departure
        self.flight_days = flight_days

        # Finding the day with minimal daylight
        self._calculate_date()
        # Defining ideal solar panel
        self._calculate_energy_power()

    def set_required_power(self, power):
        if power <= 0:
            raise ValueError("Required power must be a positive value greater than 0.")
        self.required_power = power

    def set_required_energy(self, energy):
        if energy <= 0:
            raise ValueError("Required power must be a positive value greater than 0.")
        self.required_energy = energy

    def _calculate_date(self):
        """
        Sets day with lowest amount of daylight

        Parameters
        ----------
        date : datetime object
            date and *solar* time
        """
        dates = [
            self.date_departure + timedelta(days=i) for i in range(self.flight_days)
        ]

        declinations = [declination(i) for i in dates]
        declinations = np.array(declinations)

        if self.latitude > 0:
            self.date = dates[np.argmin(declinations)]
        else:
            self.date = dates[np.argmax(declinations)]

    def _calculate_energy_power(self):
        """
        Returns the average output energy of an ideal 1m2 solar panel on day with least of sunshine
        """
        self.t = [self.date + timedelta(minutes=i) for i in range(0, 24 * 60, 1)]

        self.power_curve_ideal = []
        for i in self.t:
            self.power_curve_ideal.append(self._irradiance_on_plane(i))

        self.power_curve_ideal = np.array(self.power_curve_ideal)
        self.power = max(self.power_curve_ideal)
        # Total energy of ideal solar panel in Wh
        self.energy = sum(self.power_curve_ideal) / (60 / 1)

    # def get_daylight_hours(self):
    #     if all(self.power_curve_ideal == 0):
    #         raise Exception(
    #             "There is no daylight, change location, flight date or duration."
    #         )
    #     power_fraction = 0.1
    #     sunrise = np.argmax(self.power_curve_ideal > 0)
    #     sunset = (
    #         len(self.power_curve_ideal)
    #         - np.argmax(self.power_curve_ideal[::-1] > 0)
    #         - 1
    #     )
    #     daylight_hours = (sunset - sunrise) / 60

    #     if daylight_hours < 2:
    #         raise Exception(
    #             "There is less than 2 hours of daylight, change location, flight date or duration."
    #         )
    #     else:
    #         return daylight_hours
    def get_daylight_hours(self, power_fraction):
        # Check if all power outputs are zero
        if all(self.power_curve_ideal == 0):
            raise Exception(
                "There is no daylight, change location, flight date or duration."
            )

        # Calculate 10% of the maximum power output
        max_power = max(self.power_curve_ideal)
        power_threshold = power_fraction * max_power

        # Find indices where power exceeds 10% of the maximum power
        valid_indices = np.where(self.power_curve_ideal > power_threshold)[0]

        # Check if there are any valid indices
        if len(valid_indices) == 0:
            raise Exception(
                "There is no time when power exceeds 10% of maximum. Adjust conditions."
            )

        # Calculate sunrise and sunset times from valid indices
        self.sunrise = valid_indices[0]
        self.sunset = valid_indices[-1]

        # Convert index difference to hours assuming one index per minute
        daylight_hours = (self.sunset - self.sunrise) / 60

        # Check if there are at least 2 hours of daylight
        if daylight_hours < 2:
            raise Exception(
                "There is less than 2 hours of daylight, change location, flight date or duration."
            )
        else:
            return daylight_hours

    def define_output(self):
        # Defining weight and surface of solar panels based on power and energy requirements

        power_surface = self.required_power * self.safety_coefficient / (self.power)
        energy_surface = (
            self.required_energy
            * self.safety_coefficient
            / (self.energy * self.efficiency)
        )
        if energy_surface > power_surface:
            self.surface = energy_surface
        else:
            self.surface = power_surface
        self.total_power = self.power * self.surface
        self.total_energy = self.energy * self.surface * self.efficiency
        self.weight = self.surface * self.surface_density
        self.power_curve_real = self.power_curve_ideal * self.surface * self.efficiency
        return self.weight, self.surface

    def _irradiance_on_plane(self, date):
        """
        Returns the solar beam irradiance on a plane defined by its unit normal vector at a certain altitude, date, time, and latitude.

        Parameters
        ----------
        tracking_type : TrackingType
            The type of tracking.
        altitude : float
            Altitude above sea level in meters.
        latitude : float
            Latitude (-90 to 90) in degrees.

        Returns
        -------
        G : float
            Beam irradiance in W/m2.
        """
        cos_theta = calculate_angle_of_incidence(
            date, self.tracking_type, self.latitude
        )

        if cos_theta <= 0:
            # No sunlight is hitting the panel (it's night or sun is behind the panel)
            return 0

        # Calculate the beam irradiance based on altitude and atmospheric conditions
        G = beam_irradiance(self.altitude, date, self.latitude) * cos_theta

        return G

    # def _irradiance_on_plane(self, date):
    #     """
    #     Returns the solar beam irradiance on a plane defined by its unit
    #     normal vector in NED frame at a certain altitude, date, *solar*
    #     time and latitude.

    #     Note: it does not take into account the diffuse irradiance

    #     Parameters
    #     ----------
    #     vnorm : array-like
    #         unit vector normal to plane
    #     h : float
    #         altitude above sea level in meters
    #     date : datetime object
    #         date and *solar* time
    #     lat : float
    #         latitude (-90 to 90) in degrees

    #     Returns
    #     -------
    #     G : float
    #         beam irradiance in W/m2
    #     """

    #     def cos_theta_calc(vnorm, vsol):
    #         # TODO: when vnorm = 0

    #         vnorm_abs = np.linalg.norm(vnorm)
    #         vsol_abs = np.linalg.norm(vsol)

    #         return np.dot(vnorm, vsol) / (vnorm_abs * vsol_abs)

    #     vsol = solar_vector_ned(date, self.latitude)

    #     if (vsol == np.array([0, 0, 0])).all():
    #         # in case there is no sun (night or permanent darkness)
    #         G = 0
    #     else:
    #         if self.tracking_type == TrackingType.NO_TRACKING:
    #             # Panel is just looking up - horizontal plane

    #             cos_theta = cos_theta_z(date, self.latitude)

    #         elif self.tracking_type == TrackingType.ONE_AXIS:
    #             # Vertical axis tracking - only azimuth
    #             if self.latitude < 0:
    #                 tilt_angle = np.deg2rad(90 + self.latitude) - declination(self.date)
    #             else:
    #                 tilt_angle = np.deg2rad(90 - self.latitude) + declination(self.date)

    #             # print(f"tilt_angle 1: {tilt_angle}")
    #             # tilt_angle = abs(np.deg2rad(self.latitude) + declination(self.date))
    #             # print(f"tilt_angle 2: {tilt_angle}")
    #             # TODO: Maybe better use solar_azimuth(date, lat) instead of declination, weird shape of graph

    #             # Solar azimuth angle calculation
    #             # azimuth_angle = solar_azimuth(date, self.latitude)
    #             # vnorm = [
    #             #     -np.cos(azimuth_angle) * np.cos(tilt_angle),
    #             #     -np.sin(azimuth_angle) * np.cos(tilt_angle),
    #             #     -np.sin(tilt_angle),
    #             # ]
    #             # print(f"vnorm 1: {vnorm}")
    #             vnorm = [vsol[0], vsol[1], -np.sin(tilt_angle)]
    #             # print(f"vnorm 2: {vnorm}")
    #             cos_theta = cos_theta_calc(vnorm, vsol)
    #         elif self.tracking_type == TrackingType.TWO_AXIS:
    #             # Two axis tracking - tilt and azimuth = maximal solar gain
    #             cos_theta = 1

    #         # If behind the panel -> no irradiance
    #         if cos_theta > 0:
    #             G = beam_irradiance(self.altitude, date, self.latitude) * cos_theta
    #         else:
    #             G = 0

    #     return G

    # def plot(self):
    #     ## Data print
    #     print(f"Required power generation:  {self.required_power} Wh")
    #     print(
    #         f"Power generation:           {self.required_power*self.safety_coefficient} Wh"
    #     )
    #     print(f"Required surface:           {self.surface} m2")
    #     print(f"Solar panel weight:         {self.weight} kg")
    #     ## PLot
    #     # Create a second y-axis for azimuth
    #     fig, ax1 = plt.subplots(figsize=(18, 6))
    #     ax2 = ax1.twinx()
    #     myFmt = mdates.DateFormatter("%H:%M")
    #     # Power curve
    #     power_curve_array = np.array(self.power_curve_ideal)
    #     multiplier = self.efficiency * self.surface * self.safety_coefficient

    #     result_array = power_curve_array * multiplier
    #     # Plot power generation on the first y-axis
    #     ax1.plot(
    #         self.t,
    #         self.power_curve_real,
    #         color="blue",
    #         label="Power",
    #     )
    #     ax1.set_xlabel("Time", fontsize=15)
    #     ax1.set_ylabel("Power (W)", color="blue", fontsize=15)
    #     ax1.xaxis.set_major_formatter(myFmt)
    #     ax1.grid(True)
    #     ax1.legend(loc="upper left")

    #     # # Plot azimuth on the second y-axis
    #     # ax2.plot(t, azimuth, color='red', label='Azimuth')
    #     # ax2.set_ylabel('Azimuth (rad)', color='red', fontsize=15)
    #     # ax2.legend(loc='upper right')

    #     plt.title("Power Generation and Azimuth Over Time", fontsize=18)
    #     plt.show()


# Supporting functions
def beam_irradiance(h, date, lat):
    """
    Returns the solar beam irradiance on a plane normal to the sun vector (not
    taking into account the diffuse component) at a certain altitude, date,
    time and latitude

    Parameters
    ----------
    h : float
        altitude above sea level in meters
    date : datetime object
        date and *solar* time
    lat : float
        latitude (-90 to 90) in degrees

    Returns
    -------
    G : float
        beam irradiance in W/m2

    Notes
    -----
    Aglietti, G.S., Redi, S., Tatnall,A.R., Markvart, T., (2009) "Harnessing
    High-Altitude Solar Power"
    """
    alpha_int = 0.30  # atmospheric extinction, simplified as constant

    # the maximum zenith angle is the one that points to the horizon
    a = 6378137  # [m] Earth equatorial axis

    # Direct calculation of the cosine of the horizon limit angle, the 1.1 is a factor to include also time before sunrise and after sunset
    cos_theta_lim = a / (a + h) - 1.1

    cos_th_z = cos_theta_z(date, lat)  # radians

    pressure_max_alt, _ = gas_parameters_altitude(h)
    pressure_0, _ = gas_parameters_altitude(0)

    if cos_th_z > cos_theta_lim:
        p_relative = pressure_max_alt / pressure_0  # pressure relation
        m = air_mass_young(cos_th_z)
        G = gon(date) * np.exp(-p_relative * m * alpha_int)

    else:
        G = 0

    return G


def calculate_angle_of_incidence(date, tracking_type, latitude):
    """
    Calculates the cosine of the angle of incidence for a solar panel.

    Parameters
    ----------
    date : datetime
        The date and time for which the angle should be calculated.
    tracking_type : TrackingType
        The type of solar tracking implemented.
    latitude : float
        The latitude of the location in degrees.

    Returns
    -------
    cos_theta : float
        The cosine of the angle of incidence of solar radiation on the panel.
    """
    solar_alt = solar_altitude(date, latitude)

    if tracking_type == TrackingType.NO_TRACKING:
        # Panel is just looking up - horizontal plane
        cos_theta = np.sin(solar_alt)
    elif tracking_type == TrackingType.ONE_AXIS:
        # Vertical axis tracking - only azimuth
        if latitude < 0:
            tilt_angle = np.deg2rad(90 + latitude) - declination(date)
        else:
            tilt_angle = np.deg2rad(90 - latitude) + declination(date)
        cos_theta = np.cos(solar_alt - tilt_angle)
    elif tracking_type == TrackingType.TWO_AXIS:
        # Two axis tracking - tilt and azimuth = maximal solar gain
        cos_theta = 1  # perfect alignment with the sun

    return cos_theta


def irradiance_on_plane(date, tracking_type, altitude, latitude):
    """
    Returns the solar beam irradiance on a plane defined by its unit normal vector at a certain altitude, date, time, and latitude.

    Parameters
    ----------
    tracking_type : TrackingType
        The type of tracking.
    altitude : float
        Altitude above sea level in meters.
    latitude : float
        Latitude (-90 to 90) in degrees.

    Returns
    -------
    G : float
        Beam irradiance in W/m2.
    """
    cos_theta = calculate_angle_of_incidence(date, tracking_type, latitude)

    if cos_theta <= 0:
        # No sunlight is hitting the panel (it's night or sun is behind the panel)
        return 0

    # Calculate the beam irradiance based on altitude and atmospheric conditions
    G = beam_irradiance(altitude, date, latitude) * cos_theta

    return G


def solar_azimuth(date, lat):
    """
    Calculate the solar azimuth angle.

    Parameters
    ----------
    date : datetime object
        The date and solar time.
    lat : float
        The latitude in degrees (-90 to 90).

    Returns
    -------
    solar_az : float
        The azimuth angle in radians, positive to the West.
    """
    if abs(lat) == 90:
        lat = np.sign(lat) * 89.999  # adjust extreme latitudes

    w = hour_angle(date)  # hour angle calculation
    dec = declination(date)  # declination of the sun
    lat_rad = np.deg2rad(lat)  # convert latitude to radians

    num = np.sin(w)
    den = np.cos(w) * np.sin(lat_rad) - np.tan(dec) * np.cos(lat_rad)

    solar_az = np.arctan2(num, den)

    return solar_az


def solar_altitude(date, lat):
    """
    * Solar altitude angle *

    Angle between the projection of the sun beam on a horizontal
    surface wrt the beam, for a date, *solar* time and latitude.

    Parameters
    ----------
    date : datetime object
        date and *solar* time
    lat : float
        latitude (-90 to 90) in degrees

    Returns
    -------
    solar_altitude : float
        altitude angle in radians
    """
    cos_th_z = cos_theta_z(date, lat)

    return np.arcsin(cos_th_z)


def cos_theta_z(date, lat):
    """
    * Zenith angle *

    Angle of incidence of the sun beam on a horizontal surface wrt the
    normal to that surface, for a date, *solar* time and latitude.

    Parameters
    ----------
    date : datetime object
        date and *solar* time
    lat : float
        latitude (-90 to 90) in degrees

    Returns
    -------
    cos_th_z : float
        zenith angle of incidence in radians
    """

    dec = declination(date)
    lat = np.deg2rad(lat)
    w = hour_angle(date)

    cos_th_z = np.sin(dec) * np.sin(lat) + np.cos(dec) * np.cos(lat) * np.cos(w)

    return cos_th_z


def declination(date):
    """
    Angular position of the Sun at solar noon on a desired date and time.
    Must comply with -23.45º < declination < 23.45º

    Parameters
    ----------
    date : datetime object
        date of interest

    Returns
    ----------
    declination : float
        declination in radians

    Source
    ----------
    Spencer, J.W., 1971. Fourier series representation of the position of the sun. Search 2 (5), 172.
    """
    B = day_angle(date)
    return (
        0.006918
        - 0.399912 * np.cos(B)
        + 0.070257 * np.sin(B)
        - 0.006758 * np.cos(2 * B)
        + 0.000907 * np.sin(2 * B)
        - 0.002679 * np.cos(3 * B)
        + 0.00148 * np.sin(3 * B)
    )


def air_mass_young(cos_th_z):
    """
    Returns the ratio between air mass crossed by a sun beam to the mass
    it would pass if the sun were in the zenith at any altitude.

    Parameters
    ----------
    cos_th_z : float
        cos of zenith angle of incidence in radians


    Returns
    -------
    m : float
        ratio

    References
    -----
    A. T. Young, "AIR-MASS AND REFRACTION," Applied Optics, vol. 33,
    pp. 1108-1110, Feb 1994.
    """
    numerator = 1.002432 * (cos_th_z) ** 2 + 0.148386 * cos_th_z + 0.0096467
    denominator = (
        cos_th_z**3 + 0.149864 * cos_th_z**2 + 0.0102963 * cos_th_z + 0.000303978
    )
    # to avoid division by zero and unreasonable values
    if denominator <= 0:
        m = float("inf")
    else:
        m = numerator / denominator

    return m


def sunset_hour_angle(date, lat):
    """
    Sunset hour angle for a date and latitude

    Note: theta_z = 90º

    Parameters
    ----------
    date : datetime object
        date (indifferent time)
    lat : float
        latitude (-90 to 90) in degrees

    Returns
    -------
    sunset_hour_angle : float
        hour angle at sunset in radians
    """

    dec = declination(date)
    lat = np.deg2rad(lat)
    cos_ws = (-1) * np.tan(lat) * np.tan(dec)

    if abs(cos_ws) > 1:
        raise Exception("Sun does not set on this date. No sunrise or sunset.")
    else:
        return np.arccos(cos_ws)


def sunrise_hour_angle(date, lat):
    """
    Sunset hour angle for a date and latitude

    Note: theta_z = -90º

    Parameters
    ----------
    date : datetime object
        date (indifferent time)
    lat : float
        latitude (-90 to 90) in degrees

    Returns
    -------
    sunrise_hour_angle : float
        hour angle at sunrise in radians
    """
    return -sunset_hour_angle(date, lat)


def daylight_hours(date, lat):
    """
    Nº of hours of light for a particular day

    Parameters
    ----------
    date : datetime object
        date (indifferent time)
    lat : float
        latitude (-90 to 90) in degrees

    Returns
    -------
    day_hours : float
        number of hours of light within the day

    Note
    ----
        http://mathforum.org/library/drmath/view/56478.html
    """

    dec = declination(date)
    lat = np.deg2rad(lat)

    tmp = -np.tan(lat) * np.tan(dec)

    if tmp < -1.0:
        return 24
    elif tmp > 1.0:
        return 0
    else:
        return 24 * (2 * np.arccos(tmp) / (2 * np.pi))


def gon(date):
    """
    Extraterrestrial radiation on a plane normal to
    the radiation on a desired date and time.

    Parameters
    ----------
    date : datetime object
        date of interest

    Returns
    -------
    gon : float
        extraterrestrial radiation in W/m2
    """
    B = day_angle(date)
    return 1366.1 * (
        1.00011
        + 0.034221 * np.cos(B)
        + 0.00128 * np.sin(B)
        + 0.000719 * np.cos(2 * B)
        + 0.000077 * np.sin(2 * B)
    )


def day_angle(date):
    """
    Day-of-the-year angle on a desired date and time.

    Parameters
    ----------
    date : datetime object
        date of interest

    Returns
    -------
    B : float
        angle of the day of the year in radians
    """
    n = day_of_the_year(date)
    return np.deg2rad((n - 1) * (360 / 365))


def day_of_the_year(date):
    """
    Returns the day of the year

    Parameters
    ----------
    date : datetime object
        date of interest

    Returns
    -------
    day : int
        day of the year (1 to 365)
    """
    if isinstance(date, datetime):
        return date.timetuple().tm_yday
    else:
        msg = "date must be a datetime object or array of datetime objects"
        raise TypeError(msg)


def hour_angle(date):
    """
    Angular displacement of the sun east-west of the local meridian for a
    date, *solar* time and latitude.
    Note: 15 degrees per hour, morning < 0 < afternoon

    Parameters
    ----------
    date : datetime object
        date and *solar* time

    Returns
    -------
    hour angle : float
        local hour angle in radians
    """
    if isinstance(date, datetime):
        w = (date.hour + (date.minute / 60) - 12) * 15
        return np.deg2rad(w)
    else:
        raise TypeError("Date must be a datetime object for calculation of hour angle.")
