energy_system_modelling_workflow/scripts/solar_angles.py

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