diff --git a/.gitignore b/.gitignore index 83b5876a..7f33cea9 100644 --- a/.gitignore +++ b/.gitignore @@ -11,4 +11,4 @@ **/.idea/ cerc_hub.egg-info /out_files -/input_files/output_buildings.geojson +/input_files/output_buildings1.geojson diff --git a/hub/city_model_structure/building.py b/hub/city_model_structure/building.py index 09b5e2b1..4f11b118 100644 --- a/hub/city_model_structure/building.py +++ b/hub/city_model_structure/building.py @@ -92,6 +92,7 @@ class Building(CityObject): logging.error('Building %s [%s] has an unexpected surface type %s.', self.name, self.aliases, surface.type) self._domestic_hot_water_peak_load = None self._fuel_consumption_breakdown = {} + self._pv_generation = {} @property def shell(self) -> Polyhedron: @@ -903,10 +904,13 @@ class Building(CityObject): if cte.COOLING in energy_system.demand_types: for generation_system in energy_system.generation_systems: fuel_breakdown[generation_system.fuel_type][cte.COOLING] = self.cooling_consumption[cte.YEAR][0] / 3600 - - - - - self._fuel_consumption_breakdown = fuel_breakdown return self._fuel_consumption_breakdown + + @property + def pv_generation(self) -> dict: + return self._pv_generation + + @pv_generation.setter + def pv_generation(self, value): + self._pv_generation = value diff --git a/hub/city_model_structure/building_demand/surface.py b/hub/city_model_structure/building_demand/surface.py index 2cd42755..65f90ae3 100644 --- a/hub/city_model_structure/building_demand/surface.py +++ b/hub/city_model_structure/building_demand/surface.py @@ -46,6 +46,8 @@ class Surface: self._vegetation = None self._percentage_shared = None self._solar_collectors_area_reduction_factor = None + self._global_irradiance_tilted = {} + self._installed_solar_collector_area = None @property def name(self): @@ -384,3 +386,35 @@ class Surface: :param value: float """ self._solar_collectors_area_reduction_factor = value + + @property + def global_irradiance_tilted(self) -> dict: + """ + Get global irradiance on a tilted surface in J/m2 + :return: dict + """ + return self._global_irradiance_tilted + + @global_irradiance_tilted.setter + def global_irradiance_tilted(self, value): + """ + Set global irradiance on a tilted surface in J/m2 + :param value: dict + """ + self._global_irradiance_tilted = value + + @property + def installed_solar_collector_area(self): + """ + Get installed solar collector area in m2 + :return: dict + """ + return self._installed_solar_collector_area + + @installed_solar_collector_area.setter + def installed_solar_collector_area(self, value): + """ + Set installed solar collector area in m2 + :return: dict + """ + self._installed_solar_collector_area = value \ No newline at end of file diff --git a/hub/data/costs/montreal_costs_completed.xml b/hub/data/costs/montreal_costs_completed.xml index 81b67768..e3728433 100644 --- a/hub/data/costs/montreal_costs_completed.xml +++ b/hub/data/costs/montreal_costs_completed.xml @@ -25,8 +25,8 @@ - 0 - 0 + 300 + 300 25 diff --git a/main.py b/main.py index 4af9643a..00c7f194 100644 --- a/main.py +++ b/main.py @@ -18,12 +18,14 @@ from scripts.costs.cost import Cost from scripts.costs.constants import SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV, SYSTEM_RETROFIT_AND_PV import hub.helpers.constants as cte from hub.exports.exports_factory import ExportsFactory +from scripts.pv_sizing_and_simulation import PVSizingSimulation +from scripts.solar_angles import CitySolarAngles # Specify the GeoJSON file path -# geojson_file = process_geojson(x=-73.5953602192335, y=45.492414530022515, diff=0.001) -file_path = (Path(__file__).parent / 'input_files' / 'output_buildings.geojson') +# geojson_file = process_geojson(x=-73.5681295982132, y=45.49218262677643, diff=0.0001) +file_path = (Path(__file__).parent / 'input_files' / 'output_buildings1.geojson') # Specify the output path for the PDF file output_path = (Path(__file__).parent / 'out_files').resolve() -# Create city object from GeoJSON file +# # Create city object from GeoJSON file city = GeometryFactory('geojson', path=file_path, height_field='height', @@ -32,10 +34,19 @@ city = GeometryFactory('geojson', function_to_hub=Dictionaries().montreal_function_to_hub_function).city # Enrich city data ConstructionFactory('nrcan', city).enrich() - UsageFactory('nrcan', city).enrich() WeatherFactory('epw', city).enrich() +ExportsFactory('obj', city, output_path).export() +ExportsFactory('sra', city, output_path).export() +sra_path = (output_path / f'{city.name}_sra.xml').resolve() +subprocess.run(['sra', str(sra_path)]) +ResultFactory('sra', city, output_path).enrich() energy_plus_workflow(city) +solar_angles = CitySolarAngles(city.name, + city.latitude, + city.longitude, + tilt_angle=45, + surface_azimuth_angle=180).calculate random_assignation.call_random(city.buildings, random_assignation.residential_systems_percentage) EnergySystemsFactory('montreal_custom', city).enrich() SystemSizing(city.buildings).montreal_custom() @@ -44,9 +55,13 @@ random_assignation.call_random(city.buildings, random_assignation.residential_ne EnergySystemsFactory('montreal_future', city).enrich() for building in city.buildings: EnergySystemsSimulationFactory('archetype1', building=building, output_path=output_path).enrich() - print(building.energy_consumption_breakdown[cte.ELECTRICITY][cte.COOLING] + - building.energy_consumption_breakdown[cte.ELECTRICITY][cte.HEATING] + - building.energy_consumption_breakdown[cte.ELECTRICITY][cte.DOMESTIC_HOT_WATER]) + pv_sizing_simulation = PVSizingSimulation(building, + solar_angles, + tilt_angle=45, + module_height=1, + module_width=2, + ghi=building.roofs[0].global_irradiance[cte.HOUR]) + pv_sizing_simulation.pv_output() new_system = new_system_results(city.buildings) # EnergySystemAnalysisReport(city, output_path).create_report(current_system, new_system) for building in city.buildings: diff --git a/pv_assessment.py b/pv_assessment.py new file mode 100644 index 00000000..ad2480dd --- /dev/null +++ b/pv_assessment.py @@ -0,0 +1,69 @@ +import pandas as pd +from scripts.geojson_creator import process_geojson +from pathlib import Path +import subprocess +from hub.imports.geometry_factory import GeometryFactory +from hub.helpers.dictionaries import Dictionaries +from hub.imports.construction_factory import ConstructionFactory +from hub.imports.usage_factory import UsageFactory +from hub.imports.weather_factory import WeatherFactory +from hub.imports.results_factory import ResultFactory +from scripts.solar_angles import CitySolarAngles +from scripts.ep_run_enrich import energy_plus_workflow +import hub.helpers.constants as cte +from hub.exports.exports_factory import ExportsFactory +from scripts.pv_sizing_and_simulation import PVSizingSimulation +# Specify the GeoJSON file path +# geojson_file = process_geojson(x=-73.5681295982132, y=45.49218262677643, diff=0.0001) +file_path = (Path(__file__).parent / 'input_files' / 'output_buildings1.geojson') +# Specify the output path for the PDF file +output_path = (Path(__file__).parent / 'out_files').resolve() +# Create city object from GeoJSON file +city = GeometryFactory('geojson', + path=file_path, + height_field='height', + year_of_construction_field='year_of_construction', + function_field='function', + function_to_hub=Dictionaries().montreal_function_to_hub_function).city +# Enrich city data +ConstructionFactory('nrcan', city).enrich() + +UsageFactory('nrcan', city).enrich() +WeatherFactory('epw', city).enrich() +ExportsFactory('obj', city, output_path).export() +ExportsFactory('sra', city, output_path).export() +sra_path = (output_path / f'{city.name}_sra.xml').resolve() +subprocess.run(['sra', str(sra_path)]) +ResultFactory('sra', city, output_path).enrich() +energy_plus_workflow(city) +solar_angles = CitySolarAngles(city.name, + city.latitude, + city.longitude, + tilt_angle=45, + surface_azimuth_angle=180).calculate +df = pd.DataFrame() +df.index = ['yearly lighting (kWh)', 'yearly appliance (kWh)', 'yearly heating (kWh)', 'yearly cooling (kWh)', + 'yearly dhw (kWh)', 'roof area (m2)', 'used area for pv (m2)', 'number of panels', 'pv production (kWh)'] +for building in city.buildings: + ghi = [x for x in building.roofs[0].global_irradiance[cte.HOUR]] + pv_sizing_simulation = PVSizingSimulation(building, + solar_angles, + tilt_angle=45, + module_height=1, + module_width=2, + ghi=ghi) + pv_sizing_simulation.pv_output() + yearly_lighting = building.lighting_electrical_demand[cte.YEAR][0] / 1000 + yearly_appliance = building.appliances_electrical_demand[cte.YEAR][0] / 1000 + yearly_heating = building.heating_demand[cte.YEAR][0] / (3.6e6 * 3) + yearly_cooling = building.cooling_demand[cte.YEAR][0] / (3.6e6 * 4.5) + yearly_dhw = building.domestic_hot_water_heat_demand[cte.YEAR][0] / 1000 + roof_area = building.roofs[0].perimeter_area + used_roof = pv_sizing_simulation.available_space() + number_of_pv_panels = pv_sizing_simulation.total_number_of_panels + yearly_pv = building.onsite_electrical_production[cte.YEAR][0] / 1000 + df[f'{building.name}'] = [yearly_lighting, yearly_appliance, yearly_heating, yearly_cooling, yearly_dhw, roof_area, + used_roof, number_of_pv_panels, yearly_pv] + +df.to_csv(output_path / 'pv.csv') + diff --git a/scripts/costs/capital_costs.py b/scripts/costs/capital_costs.py index 68d751d8..6b1e159d 100644 --- a/scripts/costs/capital_costs.py +++ b/scripts/costs/capital_costs.py @@ -68,7 +68,10 @@ class CapitalCosts(CostBase): self._own_capital = 1 - self._configuration.percentage_credit self._surface_pv = 0 for roof in self._building.roofs: - self._surface_pv += roof.solid_polygon.area * roof.solar_collectors_area_reduction_factor + if roof.installed_solar_collector_area is not None: + self._surface_pv += roof.installed_solar_collector_area + else: + self._surface_pv += roof.solid_polygon.area * roof.solar_collectors_area_reduction_factor def calculate(self) -> tuple[pd.DataFrame, pd.DataFrame]: if self._configuration.retrofit_scenario in (SKIN_RETROFIT, SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV): @@ -157,6 +160,7 @@ class CapitalCosts(CostBase): capital_cost_distribution_equipment = 0 capital_cost_lighting = 0 capital_cost_pv = self._surface_pv * chapter.item('D2010_photovoltaic_system').initial_investment[0] + # capital_cost_lighting = self._total_floor_area * \ # chapter.item('D5020_lighting_and_branch_wiring').initial_investment[0] for (i, component) in enumerate(system_components): diff --git a/scripts/costs/total_maintenance_costs.py b/scripts/costs/total_maintenance_costs.py index 81a88de8..b2033d19 100644 --- a/scripts/costs/total_maintenance_costs.py +++ b/scripts/costs/total_maintenance_costs.py @@ -39,8 +39,10 @@ class TotalMaintenanceCosts(CostBase): archetype = self._archetype # todo: change area pv when the variable exists roof_area = 0 + surface_pv = 0 for roof in building.roofs: roof_area += roof.solid_polygon.area + # surface_pv += roof.installed_solar_collector_area surface_pv = roof_area * 0.5 peak_heating = building.heating_peak_load[cte.YEAR][0] / 3.6e6 diff --git a/scripts/costs/total_operational_incomes.py b/scripts/costs/total_operational_incomes.py index 2a110761..50270f18 100644 --- a/scripts/costs/total_operational_incomes.py +++ b/scripts/costs/total_operational_incomes.py @@ -37,7 +37,7 @@ class TotalOperationalIncomes(CostBase): for year in range(1, self._configuration.number_of_years + 1): price_increase_electricity = math.pow(1 + self._configuration.electricity_price_index, year) # todo: check the adequate assignation of price. Pilar - price_export = archetype.income.electricity_export * cte.WATTS_HOUR_TO_JULES * 1000 # to account for unit change + price_export = archetype.income.electricity_export # to account for unit change self._yearly_operational_incomes.loc[year, 'Incomes electricity'] = ( onsite_electricity_production * price_export * price_increase_electricity ) diff --git a/scripts/ep_run_enrich.py b/scripts/ep_run_enrich.py index 24ee4b11..ece4ffa4 100644 --- a/scripts/ep_run_enrich.py +++ b/scripts/ep_run_enrich.py @@ -15,9 +15,6 @@ def energy_plus_workflow(city): out_path = (Path(__file__).parent.parent / 'out_files') files = glob.glob(f'{out_path}/*') - # for file in files: - # if file != '.gitignore': - # os.remove(file) area = 0 volume = 0 for building in city.buildings: diff --git a/scripts/geojson_creator.py b/scripts/geojson_creator.py index f3c28b2e..864320b5 100644 --- a/scripts/geojson_creator.py +++ b/scripts/geojson_creator.py @@ -10,7 +10,7 @@ def process_geojson(x, y, diff): [x - diff, y + diff], [x + diff, y + diff]]) geojson_file = Path('./data/collinear_clean 2.geojson').resolve() - output_file = Path('./input_files/output_buildings.geojson').resolve() + output_file = Path('./input_files/output_buildings1.geojson').resolve() buildings_in_region = [] with open(geojson_file, 'r') as file: diff --git a/scripts/pv_sizing_and_simulation.py b/scripts/pv_sizing_and_simulation.py new file mode 100644 index 00000000..b75f5fbc --- /dev/null +++ b/scripts/pv_sizing_and_simulation.py @@ -0,0 +1,58 @@ +import math + +from scripts.radiation_tilted import RadiationTilted +import hub.helpers.constants as cte +from hub.helpers.monthly_values import MonthlyValues + + +class PVSizingSimulation(RadiationTilted): + def __init__(self, building, solar_angles, tilt_angle, module_height, module_width, ghi): + super().__init__(building, solar_angles, tilt_angle, ghi) + self.module_height = module_height + self.module_width = module_width + self.total_number_of_panels = 0 + self.enrich() + + def available_space(self): + roof_area = self.building.roofs[0].perimeter_area + maintenance_factor = 0.1 + orientation_factor = 0.2 + if self.building.function == cte.RESIDENTIAL: + mechanical_equipment_factor = 0.2 + else: + mechanical_equipment_factor = 0.3 + available_roof = (maintenance_factor + orientation_factor + mechanical_equipment_factor) * roof_area + return available_roof + + def inter_row_spacing(self): + winter_solstice = self.df[(self.df['AST'].dt.month == 12) & + (self.df['AST'].dt.day == 21) & + (self.df['AST'].dt.hour == 12)] + solar_altitude = winter_solstice['solar altitude'].values[0] + solar_azimuth = winter_solstice['solar azimuth'].values[0] + distance = ((self.module_height * abs(math.cos(math.radians(solar_azimuth)))) / + math.tan(math.radians(solar_altitude))) + distance = float(format(distance, '.1f')) + return distance + + def number_of_panels(self, available_roof, inter_row_distance): + space_dimension = math.sqrt(available_roof) + space_dimension = float(format(space_dimension, '.2f')) + panels_per_row = math.ceil(space_dimension / self.module_width) + number_of_rows = math.ceil(space_dimension / inter_row_distance) + self.total_number_of_panels = panels_per_row * number_of_rows + return panels_per_row, number_of_rows + + def pv_output(self): + radiation = self.total_radiation_tilted + pv_module_area = self.module_width * self.module_height + available_roof = self.available_space() + inter_row_spacing = self.inter_row_spacing() + self.number_of_panels(available_roof, inter_row_spacing) + self.building.roofs[0].installed_solar_collector_area = pv_module_area * self.total_number_of_panels + system_efficiency = 0.2 + pv_hourly_production = [x * system_efficiency * self.total_number_of_panels * pv_module_area for x in radiation] + self.building.onsite_electrical_production[cte.HOUR] = pv_hourly_production + self.building.onsite_electrical_production[cte.MONTH] = ( + MonthlyValues.get_total_month(self.building.onsite_electrical_production[cte.HOUR])) + self.building.onsite_electrical_production[cte.YEAR] = [sum(self.building.onsite_electrical_production[cte.MONTH])] \ No newline at end of file diff --git a/scripts/radiation_tilted.py b/scripts/radiation_tilted.py new file mode 100644 index 00000000..d1fabe11 --- /dev/null +++ b/scripts/radiation_tilted.py @@ -0,0 +1,112 @@ +import pandas as pd +import math +import hub.helpers.constants as cte +from hub.helpers.monthly_values import MonthlyValues + + +class RadiationTilted: + def __init__(self, building, solar_angles, tilt_angle, ghi, solar_constant=1366.1, maximum_clearness_index=1, + min_cos_zenith=0.065, maximum_zenith_angle=87): + self.building = building + self.ghi = ghi + self.tilt_angle = tilt_angle + self.zeniths = solar_angles['zenith'].tolist()[:-1] + self.incidents = solar_angles['incident angle'].tolist()[:-1] + self.date_time = solar_angles['DateTime'].tolist()[:-1] + self.ast = solar_angles['AST'].tolist()[:-1] + self.solar_azimuth = solar_angles['solar azimuth'].tolist()[:-1] + self.solar_altitude = solar_angles['solar altitude'].tolist()[:-1] + data = {'DateTime': self.date_time, 'AST': self.ast, 'solar altitude': self.solar_altitude, 'zenith': self.zeniths, + 'solar azimuth': self.solar_azimuth, 'incident angle': self.incidents, 'ghi': self.ghi} + self.df = pd.DataFrame(data) + self.df['DateTime'] = pd.to_datetime(self.df['DateTime']) + self.df['AST'] = pd.to_datetime(self.df['AST']) + self.df.set_index('DateTime', inplace=True) + 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 + self.i_on = [] + self.i_oh = [] + self.k_t = [] + self.fraction_diffuse = [] + self.diffuse_horizontal = [] + self.beam_horizontal = [] + self.dni = [] + self.beam_tilted = [] + self.diffuse_tilted = [] + self.total_radiation_tilted = [] + self.calculate() + + def dni_extra(self): + for i in range(len(self.df)): + self.i_on.append(self.solar_constant * (1 + 0.033 * math.cos(math.radians(360 * self.df.index.dayofyear[i] / 365)))) + + self.df['extraterrestrial normal radiation (Wh/m2)'] = self.i_on + + def clearness_index(self): + for i in range(len(self.df)): + self.i_oh.append(self.i_on[i] * max(math.cos(math.radians(self.zeniths[i])), self.min_cos_zenith)) + self.k_t.append(self.ghi[i] / self.i_oh[i]) + self.k_t[i] = max(0, self.k_t[i]) + self.k_t[i] = min(self.maximum_clearness_index, self.k_t[i]) + self.df['extraterrestrial radiation on horizontal (Wh/m2)'] = self.i_oh + self.df['clearness index'] = self.k_t + + def diffuse_fraction(self): + for i in range(len(self.df)): + if self.k_t[i] <= 0.22: + self.fraction_diffuse.append(1 - 0.09 * self.k_t[i]) + elif self.k_t[i] <= 0.8: + self.fraction_diffuse.append(0.9511 - 0.1604 * self.k_t[i] + 4.388 * self.k_t[i] ** 2 - + 16.638 * self.k_t[i] ** 3 + 12.336 * self.k_t[i] ** 4) + else: + self.fraction_diffuse.append(0.165) + if self.zeniths[i] > self.maximum_zenith_angle: + self.fraction_diffuse[i] = 1 + + self.df['diffuse fraction'] = self.fraction_diffuse + + def radiation_components_horizontal(self): + for i in range(len(self.df)): + self.diffuse_horizontal.append(self.ghi[i] * self.fraction_diffuse[i]) + self.beam_horizontal.append(self.ghi[i] - self.diffuse_horizontal[i]) + self.dni.append((self.ghi[i] - self.diffuse_horizontal[i]) / math.cos(math.radians(self.zeniths[i]))) + if self.zeniths[i] > self.maximum_zenith_angle or self.dni[i] < 0: + self.dni[i] = 0 + + self.df['diffuse horizontal (Wh/m2)'] = self.diffuse_horizontal + self.df['dni (Wh/m2)'] = self.dni + self.df['beam horizontal (Wh/m2)'] = self.beam_horizontal + + def radiation_components_tilted(self): + for i in range(len(self.df)): + self.beam_tilted.append(self.dni[i] * math.cos(math.radians(self.incidents[i]))) + self.beam_tilted[i] = max(self.beam_tilted[i], 0) + self.diffuse_tilted.append(self.diffuse_horizontal[i] * ((1 + math.cos(math.radians(self.tilt_angle))) / 2)) + self.total_radiation_tilted.append(self.beam_tilted[i] + self.diffuse_tilted[i]) + + self.df['beam tilted (Wh/m2)'] = self.beam_tilted + self.df['diffuse tilted (Wh/m2)'] = self.diffuse_tilted + self.df['total radiation tilted (Wh/m2)'] = self.total_radiation_tilted + + def calculate(self) -> pd.DataFrame: + self.dni_extra() + self.clearness_index() + self.diffuse_fraction() + self.radiation_components_horizontal() + self.radiation_components_tilted() + return self.df + + def enrich(self): + tilted_radiation = self.total_radiation_tilted + self.building.roofs[0].global_irradiance_tilted[cte.HOUR] = [x * cte.WATTS_HOUR_TO_JULES for x in + tilted_radiation] + self.building.roofs[0].global_irradiance_tilted[cte.HOUR] = [x * cte.WATTS_HOUR_TO_JULES for x in + tilted_radiation] + self.building.roofs[0].global_irradiance_tilted[cte.MONTH] = ( + MonthlyValues.get_total_month(self.building.roofs[0].global_irradiance_tilted[cte.HOUR])) + self.building.roofs[0].global_irradiance_tilted[cte.YEAR] = \ + [sum(self.building.roofs[0].global_irradiance_tilted[cte.MONTH])] + + diff --git a/scripts/solar_angles.py b/scripts/solar_angles.py new file mode 100644 index 00000000..8361090d --- /dev/null +++ b/scripts/solar_angles.py @@ -0,0 +1,138 @@ +import math +import pandas as pd +from datetime import datetime +from pathlib import Path + + +class CitySolarAngles: + def __init__(self, file_name, location_latitude, location_longitude, tilt_angle, surface_azimuth_angle, + standard_meridian=-75): + self.file_name = file_name + self.location_latitude = location_latitude + self.location_longitude = location_longitude + self.location_latitude_rad = math.radians(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 = (location_longitude - standard_meridian) * 4 + self.timezone = 'Etc/GMT+5' + + self.eot = [] + self.ast = [] + self.hour_angles = [] + self.declinations = [] + self.solar_altitudes = [] + self.solar_azimuths = [] + self.zeniths = [] + self.incidents = [] + self.beam_tilted = [] + self.factor = [] + self.times = pd.date_range(start='2023-01-01', end='2024-01-01', freq='H', tz=self.timezone) + self.df = pd.DataFrame(index=self.times) + self.day_of_year = self.df.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) + 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 + + @property + def calculate(self) -> pd.DataFrame: + 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) + incident_angle_radian = self.incident_angle(solar_altitude_radians, solar_azimuth_radians) + + self.df['DateTime'] = self.times + self.df['AST'] = self.ast + self.df['hour angle'] = self.hour_angles + self.df['eot'] = self.eot + self.df['declination angle'] = self.declinations + self.df['solar altitude'] = self.solar_altitudes + self.df['zenith'] = self.zeniths + self.df['solar azimuth'] = self.solar_azimuths + self.df['incident angle'] = self.incidents + + return self.df \ No newline at end of file