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: 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): """ A class to calculate the solar angles and solar irradiance on a tilted surface in the City :param city: An object from the City class -> City :param tilt_angle: tilt angle of surface -> float :param surface_azimuth_angle: The orientation of the surface. 0 is North -> float :param standard_meridian: A standard meridian is the meridian whose mean solar time is the basis of the time of day observed in a time zone -> float :param solar_constant: The amount of energy received by a given area one astronomical unit away from the Sun. It is constant and must not be changed :param maximum_clearness_index: This is used to calculate the diffuse fraction of the solar irradiance -> float :param min_cos_zenith: This is needed to avoid unrealistic values in tilted irradiance calculations -> float :param maximum_zenith_angle: This is needed to avoid negative values in tilted irradiance calculations -> float """ 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): 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): 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): 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): 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): 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): 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): 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): 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): 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): 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): 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): 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: for roof in building.roofs: hourly_tilted_irradiance = [] roof_ghi = roof.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]))) roof.global_irradiance_tilted[cte.HOUR] = hourly_tilted_irradiance roof.global_irradiance_tilted[cte.MONTH] = (MonthlyValues.get_total_month( roof.global_irradiance_tilted[cte.HOUR])) roof.global_irradiance_tilted[cte.YEAR] = [sum(roof.global_irradiance_tilted[cte.MONTH])]