CityBEM-CityLayers-SaeedRay.../scripts/solar_angles.py

147 lines
5.6 KiB
Python

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