223 lines
10 KiB
Python
223 lines
10 KiB
Python
|
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])]
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
|