calc_irradiance.py

"""
MIT License

Copyright (c) 2023 RWTH Aachen University

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""

import numpy as np
import pandas as pd
from datetime import timedelta

def calc_total_irradiance(irradiance, timer, beta, psi_f, phi, lambda_st, lambda_1, alpha_albedo=0.2, reflection=1, refraction_index=2):
    """
    check value: alpha_albedo: ground reflectivity
    Reflection = 1?
    Refraction Index = 2?

    beta = plant_parameters[(self.name, self.type)]['beta']
    psi_f = plant_parameters[(self.name, self.type)]['psi_f']
    phi = plant_parameters[(self.name, self.type)]['phi']
    lambda_st = plant_parameters[(self.name, self.type)]['lambda_st']
    lambda_1 = plant_parameters[(self.name, self.type)]['lambda_1']
    """

    def calc_direct_beam_incidence_factor():
        """
        Angle calculation taken from @Erstellung eines Modells zur Simulation der Solarstrahlung auf beliebig
        orientierte Flaechen und deren Trennung in Diffus- und Direktanteil@ Peter Ritzenhoff, March 1992

        :param timer: simulation time in TimeStamp format, it will be transformed to seconds of the year
        :param lambda_1: degree of longitude relative to Greenwich
        :param lambda_st: degree of longitude of defined time zone
        :param phi: degree of latitude
        :param psi_f: azimuth of the pv (see pic.2.4)
        :param beta: surface inclination tilt angle(Neigungswinkel)
        :return:
        """

        daysperyear = 365
        if timer.is_leap_year:
            daysperyear = 366

        real_timer = (timer.to_pydatetime() - pd.Timestamp(
            str(timer.year) + "-01-01 00:00:00").to_pydatetime()).total_seconds()
        # The day number tn, see eq. 2.1
        tn = np.mod(real_timer / 86400, daysperyear)
        # The day angle tw in degree, see eq. 2.3
        tw = np.mod(tn, daysperyear) * 360 / daysperyear
        # Elliptical description for delta see eq.2.2
        delta_rad = np.arcsin(0.3978 * np.sin(tw - 80.2 + 1.92 * np.sin(np.deg2rad(tw - 2.8))))
        # calculate the hour equation zg, see eq.2.6 (hour equation is for Greenwich)
        zg = -0.128 * np.sin(np.deg2rad(tw - 2.8)) - 0.165 * np.sin(np.deg2rad(2 * tw + 19.7))
        # calculate local standard time, time of the day
        lst = real_timer / 3600 - np.floor(tn) * 24
        # calculate the solar time, see eq.2.7
        sz = np.mod(lst + (lambda_1 - lambda_st) / 15 + zg, 24)
        # calculate the hour angle omega, see eq.2.8
        omega_rad = np.deg2rad(15 * (sz - 12))
        phi_rad = np.deg2rad(phi)
        # calculate the sun elevation gamma, see eq.2.9
        gamma_rad = np.arcsin(
            np.sin(phi_rad) * np.sin(delta_rad) + np.cos(omega_rad) * np.cos(phi_rad) * np.cos(delta_rad))
        # calculate the sun azimuth, see eq.2.10
        psi_s_rad = np.arccos(
            (np.sin(phi_rad) * np.sin(gamma_rad) - np.sin(delta_rad)) / (np.cos(phi_rad) * np.cos(gamma_rad)))
        psi_s_rad = -psi_s_rad if 0 <= sz < 12 else psi_s_rad
        # calculate the Surface(PV)-sun-azimuth, see eq.11
        psi_fs_rad = np.deg2rad(np.rad2deg(psi_s_rad) - psi_f)
        beta_rad_local = np.deg2rad(beta)
        # calculate the angle of incidence of the sun relative to the surface, see eq.2.12
        cos_theta = np.cos(gamma_rad) * np.sin(beta_rad_local) * np.cos(psi_fs_rad) + np.sin(gamma_rad) * np.cos(
            beta_rad_local)
        theta_local = np.rad2deg(np.arccos(cos_theta))
        # perhaps a good method, DNI will be considered once the sun height reaches 4 degrees
        f_beam_local = max(cos_theta, 0)
        return theta_local, f_beam_local

    # g_total = g_diffuse + g_reflected + g_beam
    beta_deg = beta
    beta_rad = np.deg2rad(beta_deg)  # tilt angle in [rad]
    R_d = (1 + np.cos(beta_rad)) / 2
    R = alpha_albedo * (1 - np.cos(beta_rad)) / 2  # alpha_albedo: ground reflectivity
    beta_table = [-0.167716, 4.91582, 9.99807, 14.8933, 19.9769, 25.0591, 30.1427, 35.0379, 39.5553, 44.6389,
                  49.7224, 54.6189, 59.892, 64.4095, 69.6826, 74.5791, 79.4744, 84.5605, 89.644]
    diffuse_reflection_loss_table = [0.94538, 0.945765, 0.949564, 0.949935, 0.950321, 0.954119, 0.954505, 0.954876,
                                     0.958631, 0.959017, 0.959403, 0.956361, 0.953348, 0.957103, 0.95409, 0.951049,
                                     0.95142, 0.94498, 0.945365]
    diffuse_reflection_loss = np.interp(beta_deg, beta_table, diffuse_reflection_loss_table)
    theta, f_beam = calc_direct_beam_incidence_factor()
    theta_out = np.rad2deg(np.arcsin(np.sin(np.deg2rad(theta / refraction_index))))
    r_ges = (np.tan(theta - theta_out) ** 2 / np.tan(theta + theta_out) ** 2) * 0.5 * (np.cos(theta) ** 2 + 1) + \
            (np.sin(theta - theta_out) ** 2 / np.sin(theta + theta_out) ** 2) * 0.5 * np.sin(theta) ** 2
    transmission_factor = max(0, 1 - r_ges)
    if reflection == 0:
        diffuse_reflection_loss = 1
        transmission_factor = 1
    g_diffuse = irradiance[0] * R_d * diffuse_reflection_loss
    g_reflected = irradiance[1] * R
    g_beam = irradiance[2] * f_beam * transmission_factor

    g_total = g_diffuse + g_reflected + g_beam
    return g_total

def generate_g_t_series(irradiance, beta, psi_f, phi, lambda_st, lambda_1, t_start, t_horizon, t_step):
    T = list(pd.date_range(pd.Timestamp(t_start), pd.Timestamp(t_start) + timedelta(hours=t_horizon * t_step - t_step), freq=str(t_step) + 'H'))
    g_total_lst = []
    for t in T:
        g_total_lst.append(calc_total_irradiance(irradiance.loc[t], t, beta, psi_f, phi, lambda_st, lambda_1))

    return pd.Series(data=g_total_lst, index=T)