pv_workflow/pv_assessment/solar_calculator.py

"""
solar_calculator module
SPDX-License-Identifier: LGPL-3.0-or-later
Copyright © 2022 Concordia CERC group
Project Coder: Saeed Ranjbar saeed.ranjbar@cerc.com
"""
import math
import pandas as pd
from datetime import datetime
import hub.helpers.constants as cte
from hub.helpers.monthly_values import MonthlyValues


class SolarCalculator:
  """
  SolarCalculator class performs solar angle and irradiance calculations for a given city and tilt angles
  """
  def __init__(self, city, tilt_angle, surface_azimuth_angle, standard_meridian=-75,
               solar_constant=1366.1, maximum_clearness_index=1, min_cos_zenith=0.065, maximum_zenith_angle=87):
    """
    Initialize SolarCalculator with city and solar panel configurations.
    :param city: City object containing latitude and longitude
    :param tilt_angle: Tilt angle of the solar panel in degrees
    :param surface_azimuth_angle: Azimuth angle of the solar panel in degrees
    :param standard_meridian: Standard meridian for time zone correction
    :param solar_constant: Extraterrestrial radiation constant in W/m2
    :param maximum_clearness_index: Maximum clearness index for calculation
    :param min_cos_zenith: Minimum cosine zenith value
    :param maximum_zenith_angle: Maximum allowable zenith angle in degrees
    """
    self.city = city
    self.location_latitude = city.latitude
    self.location_longitude = city.longitude
    self.location_latitude_rad = math.radians(self.location_latitude)
    self.surface_azimuth_angle = surface_azimuth_angle
    self.surface_azimuth_rad = math.radians(surface_azimuth_angle)
    self.tilt_angle = tilt_angle
    self.tilt_angle_rad = math.radians(tilt_angle)
    self.standard_meridian = standard_meridian
    self.longitude_correction = (self.location_longitude - standard_meridian) * 4
    self.solar_constant = solar_constant
    self.maximum_clearness_index = maximum_clearness_index
    self.min_cos_zenith = min_cos_zenith
    self.maximum_zenith_angle = maximum_zenith_angle
    timezone_offset = int(-standard_meridian / 15)
    self.timezone = f'Etc/GMT{"+" if timezone_offset < 0 else "-"}{abs(timezone_offset)}'
    self.eot = []
    self.ast = []
    self.hour_angles = []
    self.declinations = []
    self.solar_altitudes = []
    self.solar_azimuths = []
    self.zeniths = []
    self.incidents = []
    self.i_on = []
    self.i_oh = []
    self.times = pd.date_range(start='2023-01-01', end='2023-12-31 23:00', freq='h', tz=self.timezone)
    self.solar_angles = pd.DataFrame(index=self.times)
    self.day_of_year = self.solar_angles.index.dayofyear

  def solar_time(self, datetime_val, day_of_year):
    """
    Calculate apparent solar time for a given datetime and day of the year.
    :param datetime_val: Input datetime value
    :param day_of_year: Day of the year (1-365)
    :return: Apparent solar time (datetime)
    """
    b = (day_of_year - 81) * 2 * math.pi / 364
    eot = 9.87 * math.sin(2 * b) - 7.53 * math.cos(b) - 1.5 * math.sin(b)
    self.eot.append(eot)

    # Calculate Local Solar Time (LST)
    lst_hour = datetime_val.hour
    lst_minute = datetime_val.minute
    lst_second = datetime_val.second
    lst = lst_hour + lst_minute / 60 + lst_second / 3600

    # Calculate Apparent Solar Time (AST) in decimal hours
    ast_decimal = lst + eot / 60 + self.longitude_correction / 60
    ast_hours = int(ast_decimal) % 24  # Adjust hours to fit within 0–23 range
    ast_minutes = round((ast_decimal - ast_hours) * 60)

    # Ensure ast_minutes is within valid range
    if ast_minutes == 60:
      ast_hours += 1
      ast_minutes = 0
    elif ast_minutes < 0:
      ast_minutes = 0
    ast_time = datetime(year=datetime_val.year, month=datetime_val.month, day=datetime_val.day,
                        hour=ast_hours, minute=ast_minutes)
    self.ast.append(ast_time)
    return ast_time

  def declination_angle(self, day_of_year):
    """
    Calculate the solar declination angle for a given day of the year.
    :param day_of_year: Day of the year (1-365)
    :return: Declination angle in radians
    """
    declination = 23.45 * math.sin(math.radians(360 / 365 * (284 + day_of_year)))
    declination_radian = math.radians(declination)
    self.declinations.append(declination)
    return declination_radian

  def hour_angle(self, ast_time):
    """
    Calculate the hour angle for a given apparent solar time (AST).

    The hour angle is a measure of time since solar noon in degrees or radians,
    where solar noon corresponds to 0°. It is negative in the morning and positive in the afternoon.

    :param ast_time: Apparent solar time as a datetime object.
    :return: Hour angle in radians.
    """
    hour_angle = ((ast_time.hour * 60 + ast_time.minute) - 720) / 4
    hour_angle_radian = math.radians(hour_angle)
    self.hour_angles.append(hour_angle)
    return hour_angle_radian

  def solar_altitude(self, declination_radian, hour_angle_radian):
    """
    Calculate the solar altitude angle in radians for a given declination and hour angle.

    The solar altitude angle is the angle between the sun's rays and the horizontal plane
    at a given location and time. It indicates the sun's height in the sky.

    :param declination_radian: Solar declination angle in radians.
    :param hour_angle_radian: Hour angle in radians.
    :return: Solar altitude angle in radians.
    """
    solar_altitude_radians = math.asin(math.cos(self.location_latitude_rad) * math.cos(declination_radian) *
                                       math.cos(hour_angle_radian) + math.sin(self.location_latitude_rad) *
                                       math.sin(declination_radian))
    solar_altitude = math.degrees(solar_altitude_radians)
    self.solar_altitudes.append(solar_altitude)
    return solar_altitude_radians

  def zenith(self, solar_altitude_radians):
    """
    Calculate the solar zenith angle in radians from the solar altitude angle.

    The solar zenith angle is the angle between the vertical direction
    (directly overhead) and the line to the sun. It is complementary to
    the solar altitude angle.

    :param solar_altitude_radians: Solar altitude angle in radians.
    :return: Solar zenith angle in radians.
    """
    solar_altitude = math.degrees(solar_altitude_radians)
    zenith_degree = 90 - solar_altitude
    zenith_radian = math.radians(zenith_degree)
    self.zeniths.append(zenith_degree)
    return zenith_radian

  def solar_azimuth_analytical(self, hourangle, declination, zenith):
    """
    Calculate the solar azimuth angle analytically in radians.

    The solar azimuth angle represents the sun's position relative to true north,
    measured clockwise. The method uses the hour angle, solar declination, and
    solar zenith to compute the azimuth.

    :param hourangle: Hour angle of the sun in radians, indicating its position
                      relative to the solar noon.
    :param declination: Solar declination angle in radians, which is the angle
                        between the sun's rays and the plane of Earth's equator.
    :param zenith: Solar zenith angle in radians, the angle between the vertical
                   direction and the sun's position.
    :return: Solar azimuth angle in radians.
    """
    numer = (math.cos(zenith) * math.sin(self.location_latitude_rad) - math.sin(declination))
    denom = (math.sin(zenith) * math.cos(self.location_latitude_rad))
    if math.isclose(denom, 0.0, abs_tol=1e-8):
      cos_azi = 1.0
    else:
      cos_azi = numer / denom

    cos_azi = max(-1.0, min(1.0, cos_azi))

    sign_ha = math.copysign(1, hourangle)
    solar_azimuth_radians = sign_ha * math.acos(cos_azi) + math.pi
    solar_azimuth_degrees = math.degrees(solar_azimuth_radians)
    self.solar_azimuths.append(solar_azimuth_degrees)
    return solar_azimuth_radians

  def incident_angle(self, solar_altitude_radians, solar_azimuth_radians):
    """
    Calculate the solar incident angle in radians.

    The incident angle represents the angle between the solar rays and the
    normal to a tilted surface. It is a critical parameter for evaluating
    the performance of solar panels.

    :param solar_altitude_radians: Solar altitude angle in radians, indicating
                                    the sun's position above the horizon.
    :param solar_azimuth_radians: Solar azimuth angle in radians, specifying
                                   the sun's position relative to true north.
    :return: Solar incident angle in radians.
    """
    incident_radian = math.acos(math.cos(solar_altitude_radians) *
                                math.cos(abs(solar_azimuth_radians - self.surface_azimuth_rad)) *
                                math.sin(self.tilt_angle_rad) + math.sin(solar_altitude_radians) *
                                math.cos(self.tilt_angle_rad))
    incident_angle_degrees = math.degrees(incident_radian)
    self.incidents.append(incident_angle_degrees)
    return incident_radian

  def dni_extra(self, day_of_year, zenith_radian):
    """
    Calculate extraterrestrial DNI and horizontal irradiance.

    :param day_of_year: Day of the year (1–365/366).
    :param zenith_radian: Solar zenith angle in radians.
    :return: Tuple (i_on, i_oh) where:
             - i_on: Extraterrestrial normal irradiance [W/m²].
             - i_oh: Extraterrestrial horizontal irradiance [W/m²].
    """
    i_on = self.solar_constant * (1 + 0.033 * math.cos(math.radians(360 * day_of_year / 365)))
    i_oh = i_on * max(math.cos(zenith_radian), self.min_cos_zenith)
    self.i_on.append(i_on)
    self.i_oh.append(i_oh)
    return i_on, i_oh

  def clearness_index(self, ghi, i_oh):
    """
    Calculate the clearness index (Kt).

    :param ghi: Global horizontal irradiance [W/m²].
    :param i_oh: Extraterrestrial horizontal irradiance [W/m²].
    :return: Clearness index (Kt).
    """
    k_t = ghi / i_oh
    k_t = max(0, k_t)
    k_t = min(self.maximum_clearness_index, k_t)
    return k_t

  def diffuse_fraction(self, k_t, zenith):
    """
    Estimate the diffuse fraction of irradiance.

    :param k_t: Clearness index (Kt).
    :param zenith: Solar zenith angle in degrees.
    :return: Diffuse fraction.
    """
    if k_t <= 0.22:
      fraction_diffuse = 1 - 0.09 * k_t
    elif k_t <= 0.8:
      fraction_diffuse = (0.9511 - 0.1604 * k_t + 4.388 * k_t ** 2 - 16.638 * k_t ** 3 + 12.336 * k_t ** 4)
    else:
      fraction_diffuse = 0.165
    if zenith > self.maximum_zenith_angle:
      fraction_diffuse = 1
    return fraction_diffuse

  def radiation_components_horizontal(self, ghi, fraction_diffuse, zenith):
    """
    Compute diffuse and beam components of horizontal radiation.

    :param ghi: Global horizontal irradiance [W/m²].
    :param fraction_diffuse: Diffuse fraction.
    :param zenith: Solar zenith angle in degrees.
    :return: Tuple (diffuse_horizontal, dni).
    """
    diffuse_horizontal = ghi * fraction_diffuse
    dni = (ghi - diffuse_horizontal) / math.cos(math.radians(zenith))
    if zenith > self.maximum_zenith_angle or dni < 0:
      dni = 0
    return diffuse_horizontal, dni

  def radiation_components_tilted(self, diffuse_horizontal, dni, incident_angle):
    """
    Compute total radiation on a tilted surface.

    :param diffuse_horizontal: Diffuse horizontal irradiance [W/m²].
    :param dni: Direct normal irradiance [W/m²].
    :param incident_angle: Solar incident angle in degrees.
    :return: Total radiation on tilted surface [W/m²].
    """
    beam_tilted = dni * math.cos(math.radians(incident_angle))
    beam_tilted = max(beam_tilted, 0)
    diffuse_tilted = diffuse_horizontal * ((1 + math.cos(math.radians(self.tilt_angle))) / 2)
    total_radiation_tilted = beam_tilted + diffuse_tilted
    return total_radiation_tilted

  def solar_angles_calculator(self, csv_output=False):
    for i in range(len(self.times)):
      datetime_val = self.times[i]
      day_of_year = self.day_of_year[i]
      declination_radians = self.declination_angle(day_of_year)
      ast_time = self.solar_time(datetime_val, day_of_year)
      hour_angle_radians = self.hour_angle(ast_time)
      solar_altitude_radians = self.solar_altitude(declination_radians, hour_angle_radians)
      zenith_radians = self.zenith(solar_altitude_radians)
      solar_azimuth_radians = self.solar_azimuth_analytical(hour_angle_radians, declination_radians, zenith_radians)
      self.incident_angle(solar_altitude_radians, solar_azimuth_radians)
      self.dni_extra(day_of_year=day_of_year, zenith_radian=zenith_radians)
    self.solar_angles['DateTime'] = self.times
    self.solar_angles['AST'] = self.ast
    self.solar_angles['hour angle'] = self.hour_angles
    self.solar_angles['eot'] = self.eot
    self.solar_angles['declination angle'] = self.declinations
    self.solar_angles['solar altitude'] = self.solar_altitudes
    self.solar_angles['zenith'] = self.zeniths
    self.solar_angles['solar azimuth'] = self.solar_azimuths
    self.solar_angles['incident angle'] = self.incidents
    self.solar_angles['extraterrestrial normal radiation (Wh/m2)'] = self.i_on
    self.solar_angles['extraterrestrial radiation on horizontal (Wh/m2)'] = self.i_oh
    if csv_output:
      self.solar_angles.to_csv('solar_angles_new.csv')

  def tilted_irradiance_calculator(self):
    if self.solar_angles.empty:
      self.solar_angles_calculator()
    for building in self.city.buildings:
      hourly_tilted_irradiance = []
      roof_ghi = building.roofs[0].global_irradiance[cte.HOUR]
      for i in range(len(roof_ghi)):
        k_t = self.clearness_index(ghi=roof_ghi[i], i_oh=self.i_oh[i])
        fraction_diffuse = self.diffuse_fraction(k_t, self.zeniths[i])
        diffuse_horizontal, dni = self.radiation_components_horizontal(ghi=roof_ghi[i],
                                                                       fraction_diffuse=fraction_diffuse,
                                                                       zenith=self.zeniths[i])
        hourly_tilted_irradiance.append(int(self.radiation_components_tilted(diffuse_horizontal=diffuse_horizontal,
                                                                             dni=dni,
                                                                             incident_angle=self.incidents[i])))

      building.roofs[0].global_irradiance_tilted[cte.HOUR] = hourly_tilted_irradiance
      building.roofs[0].global_irradiance_tilted[cte.MONTH] = (MonthlyValues.get_total_month(
        building.roofs[0].global_irradiance_tilted[cte.HOUR]))
      building.roofs[0].global_irradiance_tilted[cte.YEAR] = [sum(building.roofs[0].global_irradiance_tilted[cte.MONTH])]