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