336 lines
14 KiB
Python
336 lines
14 KiB
Python
"""
|
||
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])]
|
||
|
||
|
||
|
||
|
||
|