From 6916a5c034d6e4250e78f319e5a9ed1d7501738e Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Thu, 20 Mar 2025 11:46:23 +0100
Subject: [PATCH 01/30] [feature] First script implementation of "CIE general
 sky standard defining luminance distributions"

Functions to compute gradation, scattering indicatrix, relative and absolute luminance.
Results validated against the example from the paper by Darula and Kittler "CIE general sky standard defining luminance distributions"
---
 TEMPORARY_FILES/perez_standard_skies.py | 156 ++++++++++++++++++++++++
 1 file changed, 156 insertions(+)
 create mode 100644 TEMPORARY_FILES/perez_standard_skies.py

diff --git a/TEMPORARY_FILES/perez_standard_skies.py b/TEMPORARY_FILES/perez_standard_skies.py
new file mode 100644
index 0000000..8a22375
--- /dev/null
+++ b/TEMPORARY_FILES/perez_standard_skies.py
@@ -0,0 +1,156 @@
+import logging
+import math
+import numpy as np
+import os
+import pandas as pd
+from matplotlib import pyplot as plt
+import matplotlib
+matplotlib.use('TkAgg')
+
+from MODULES.conversion_functions import sph_to_cart
+
+logger = logging.getLogger(__name__)
+
+def compute_sph_angular_dist(Z, Z_s, az, az_s):
+    """
+
+    Args:
+        Z: Sky patch zenith angle [deg]
+        Z_s: Solar zenith angle [deg]
+        az: Sky patch azimuth angle from north [deg]
+        az_s: Solar azimuth angle from north [deg]
+
+    Returns:
+        chi: Spherical angular distance [rad]
+    """
+
+    # Convert all angles from degrees to radians
+    Z = np.deg2rad(Z)
+    Z_s = np.deg2rad(Z_s)
+    az = np.deg2rad(az)
+    az_s = np.deg2rad(az_s)
+
+    # Absolute difference in azimuth angles
+    Az = np.abs(az-az_s)
+
+    # Angular distance chi
+    chi = np.arccos(np.cos(Z_s) * np.cos(Z) + np.sin(Z_s)*np.sin(Z)*np.cos(Az))
+
+    return chi
+
+def gradation_fun(a, b, Z):
+    """
+    Gradation function defined to be vectorized in compute_gradation().
+    Args:
+        a: gradation parameter [-]
+        b: gradation parameter [-]
+        Z: zenith angle [rad]
+
+    Returns:
+        Gradation value.
+    """
+    if math.isclose(Z, np.pi/2):
+        return 1
+    else:
+        return 1 + a*np.exp(b/np.cos(Z))
+
+def compute_gradation(a, b, Z):
+    """
+    Computes the gradation function based on vectorized gradation_fun function.
+    Args:
+        a: gradation parameter [-]
+        b: gradation parameter [-]
+        Z: Zenith angle [deg]. Type : Nx1 array
+
+    Returns:
+        gradation function result
+    """
+    v_gradation = np.vectorize(gradation_fun)
+
+    Z = np.deg2rad(Z)
+
+    return v_gradation(a, b, Z)
+
+def scattering_indicatrix_fun(c, d, e, x):
+    """
+
+    Args:
+        c: scattering parameter [-]
+        d: scattering parameter [-]
+        e: scattering parameter [-]
+        x: angle (sphericla angular distance or zenith, i.e. 0) [rad]
+    Returns:
+        Scattering indicatrix value
+    """
+    return 1 + c * (np.exp(d*x) - np.exp(d*np.pi/2)) + e * (np.cos(x))**2
+
+def get_sky_params(standard_sky):
+    """
+
+    Args:
+        standard_sky: Dataframe of length 1 with standard sky information
+
+    Returns:
+        a: gradation parameter [-]
+        b: gradation parameter [-]
+        c: scattering parameter [-]
+        d: scattering parameter [-]
+        e: scattering parameter [-]
+    """
+
+    a = standard_sky['a'].values[0]
+    b = standard_sky['b'].values[0]
+    c = standard_sky['c'].values[0]
+    d = standard_sky['d'].values[0]
+    e = standard_sky['e'].values[0]
+
+    return a, b, c, d, e
+
+def luminance_fun(sky_type, patch_az, patch_el, sun_az, sun_el, zenith_lum):
+    # Standard sky parameters and description
+    standard_sky = standard_skies.loc[standard_skies['Type'] == sky_type]
+
+    # Sky model parameters
+    a, b, c, d, e = get_sky_params(standard_sky)
+
+    # Elevation to zenith (sky patch)
+    Z = 90 - patch_el  # [deg]
+
+    # Elevation to zenith (Sun)
+    Z_s = 90 - sun_el  # [deg]
+
+    # Spherical angular dist
+    chi = compute_sph_angular_dist(Z, Z_s, patch_az, sun_az)
+
+    # Gradation
+    phi = compute_gradation(a, b, Z)
+    phi_zenith = compute_gradation(a, b, Z=0)
+
+    # Indicatrix
+    indicatrix = scattering_indicatrix_fun(c, d, e, chi)
+    indicatrix_zenith = scattering_indicatrix_fun(c, d, e, np.deg2rad(Z_s))
+
+    # Relative luminance
+    rel_lum = phi * indicatrix / (phi_zenith * indicatrix_zenith)
+
+    # Absolute luminance
+    lum = rel_lum * zenith_luminance
+
+    return lum
+
+standard_skies = pd.read_csv(os.path.join('..', 'INPUTS', 'CIE_standard_skies.csv'))
+
+print(compute_gradation(1, 1, 90) == 1)
+
+# Example from the paper
+sky_type = 11
+
+az = np.array([130, 120])  # [°]
+el = np.array([10, 20])  # [°]
+az, el = np.meshgrid(az, el)
+
+el_s = 38.02  # [°]
+az_s = 147.67  # [°]
+zenith_luminance = 4404  # cd/m²
+
+lum = luminance_fun(sky_type, az, el, az_s, el_s, zenith_luminance)
-- 
GitLab


From 0beb7e4a34dfb163022a75bb1cf85a7e317b5a49 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Thu, 20 Mar 2025 13:37:24 +0100
Subject: [PATCH 02/30] [feature] Vectorize the luminance function

To accept vectors of azimuth, elevations as inputs.
---
 TEMPORARY_FILES/perez_standard_skies.py | 13 +++++++++----
 1 file changed, 9 insertions(+), 4 deletions(-)

diff --git a/TEMPORARY_FILES/perez_standard_skies.py b/TEMPORARY_FILES/perez_standard_skies.py
index 8a22375..d30bc39 100644
--- a/TEMPORARY_FILES/perez_standard_skies.py
+++ b/TEMPORARY_FILES/perez_standard_skies.py
@@ -138,6 +138,9 @@ def luminance_fun(sky_type, patch_az, patch_el, sun_az, sun_el, zenith_lum):
 
     return lum
 
+def compute_luminance(sky_type, patch_az, patch_el, sun_az, sun_el, zenith_lum):
+    v_luminance_fun = np.vectorize(luminance_fun)
+    return v_luminance_fun(sky_type, patch_az, patch_el, sun_az, sun_el, zenith_lum)
 standard_skies = pd.read_csv(os.path.join('..', 'INPUTS', 'CIE_standard_skies.csv'))
 
 print(compute_gradation(1, 1, 90) == 1)
@@ -145,12 +148,14 @@ print(compute_gradation(1, 1, 90) == 1)
 # Example from the paper
 sky_type = 11
 
-az = np.array([130, 120])  # [°]
-el = np.array([10, 20])  # [°]
-az, el = np.meshgrid(az, el)
+az = np.arange(0, 180, step=10)  # [°]
+el = np.arange(0, 90, step=5)  # [°]
+az_meshgrid, el_meshgrid = np.meshgrid(az, el)
 
 el_s = 38.02  # [°]
-az_s = 147.67  # [°]
+az_s = 147.67  # [°] 145° is south-east
 zenith_luminance = 4404  # cd/m²
 
 lum = luminance_fun(sky_type, az, el, az_s, el_s, zenith_luminance)
+lum = compute_luminance(sky_type, az_meshgrid.flatten(), el_meshgrid.flatten(), az_s, el_s, zenith_luminance)
+
-- 
GitLab


From 44fac173f4ed010fe1aba58b67002c15e2270de2 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Thu, 20 Mar 2025 13:38:09 +0100
Subject: [PATCH 03/30] [feature] Plot the result with pcolor in a rectangular
 grid

Plot the luminance distribution and sun position with pcolor
---
 TEMPORARY_FILES/perez_standard_skies.py | 30 +++++++++++++++++++++++--
 1 file changed, 28 insertions(+), 2 deletions(-)

diff --git a/TEMPORARY_FILES/perez_standard_skies.py b/TEMPORARY_FILES/perez_standard_skies.py
index d30bc39..6d51d58 100644
--- a/TEMPORARY_FILES/perez_standard_skies.py
+++ b/TEMPORARY_FILES/perez_standard_skies.py
@@ -134,13 +134,38 @@ def luminance_fun(sky_type, patch_az, patch_el, sun_az, sun_el, zenith_lum):
     rel_lum = phi * indicatrix / (phi_zenith * indicatrix_zenith)
 
     # Absolute luminance
-    lum = rel_lum * zenith_luminance
+    lum = rel_lum * zenith_lum
 
     return lum
 
 def compute_luminance(sky_type, patch_az, patch_el, sun_az, sun_el, zenith_lum):
     v_luminance_fun = np.vectorize(luminance_fun)
     return v_luminance_fun(sky_type, patch_az, patch_el, sun_az, sun_el, zenith_lum)
+
+def plot_luminance_map(az, el, lum, az_s=None, el_s=None, luminance_or_radiance='radiance'):
+    lum_min = np.min(lum)
+    lum_max = np.max(lum)
+
+    fig = plt.figure()
+    ax = fig.gca()
+
+    c = ax.pcolor(az, el, lum, vmin=lum_min, vmax=lum_max)
+    ax.set_xlabel('Azimuth [deg]')
+    ax.set_ylabel('Elevation [deg]')
+
+    cbar = fig.colorbar(c, ax=ax)
+    if luminance_or_radiance == 'luminance':
+        cbar_label = 'Luminance [cd/m²]'
+    elif luminance_or_radiance == 'radiance':
+        cbar_label = 'Radiance [W/(m²⋅sr)]'
+    cbar.set_label(cbar_label)
+
+    if (az_s is not None) and (el_s is not None):
+        ax.plot(az_s, el_s, 'ro', label='Sun position')
+        ax.legend()
+
+    plt.show()
+
 standard_skies = pd.read_csv(os.path.join('..', 'INPUTS', 'CIE_standard_skies.csv'))
 
 print(compute_gradation(1, 1, 90) == 1)
@@ -156,6 +181,7 @@ el_s = 38.02  # [°]
 az_s = 147.67  # [°] 145° is south-east
 zenith_luminance = 4404  # cd/m²
 
-lum = luminance_fun(sky_type, az, el, az_s, el_s, zenith_luminance)
 lum = compute_luminance(sky_type, az_meshgrid.flatten(), el_meshgrid.flatten(), az_s, el_s, zenith_luminance)
 
+# Plot with pcolor
+plot_luminance_map(az, el, lum.reshape(az_meshgrid.shape), az_s, el_s, luminance_or_radiance='luminance')
-- 
GitLab


From 70ada4689873b9ab1ed78843749fb1b877673cd2 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Thu, 20 Mar 2025 15:02:51 +0100
Subject: [PATCH 04/30] [docs] More accurate description for parameters a, b,
 c, d, e

From A. Marsh website
---
 TEMPORARY_FILES/perez_standard_skies.py | 24 ++++++++++++------------
 1 file changed, 12 insertions(+), 12 deletions(-)

diff --git a/TEMPORARY_FILES/perez_standard_skies.py b/TEMPORARY_FILES/perez_standard_skies.py
index 6d51d58..2bd7b69 100644
--- a/TEMPORARY_FILES/perez_standard_skies.py
+++ b/TEMPORARY_FILES/perez_standard_skies.py
@@ -42,8 +42,8 @@ def gradation_fun(a, b, Z):
     """
     Gradation function defined to be vectorized in compute_gradation().
     Args:
-        a: gradation parameter [-]
-        b: gradation parameter [-]
+        a: Horizon-zenith gradient (gradation parameter) [-]
+        b: Gradient intensity (gradation parameter) [-]
         Z: zenith angle [rad]
 
     Returns:
@@ -58,8 +58,8 @@ def compute_gradation(a, b, Z):
     """
     Computes the gradation function based on vectorized gradation_fun function.
     Args:
-        a: gradation parameter [-]
-        b: gradation parameter [-]
+        a: Horizon-zenith gradient (gradation parameter) [-]
+        b: Gradient intensity (gradation parameter) [-]
         Z: Zenith angle [deg]. Type : Nx1 array
 
     Returns:
@@ -75,9 +75,9 @@ def scattering_indicatrix_fun(c, d, e, x):
     """
 
     Args:
-        c: scattering parameter [-]
-        d: scattering parameter [-]
-        e: scattering parameter [-]
+        c: Circumsolar intensity (scattering parameter) [-]
+        d: Circumsolar radius (scattering parameter) [-]
+        e: Backscattering effect (scattering parameter) [-]
         x: angle (sphericla angular distance or zenith, i.e. 0) [rad]
     Returns:
         Scattering indicatrix value
@@ -91,11 +91,11 @@ def get_sky_params(standard_sky):
         standard_sky: Dataframe of length 1 with standard sky information
 
     Returns:
-        a: gradation parameter [-]
-        b: gradation parameter [-]
-        c: scattering parameter [-]
-        d: scattering parameter [-]
-        e: scattering parameter [-]
+        a: Horizon-zenith gradient (gradation parameter) [-]
+        b: Gradient intensity (gradation parameter) [-]
+        c: Circumsolar intensity (scattering parameter) [-]
+        d: Circumsolar radius (scattering parameter) [-]
+        e: Backscattering effect (scattering parameter) [-]
     """
 
     a = standard_sky['a'].values[0]
-- 
GitLab


From 965c9e2f8a4a26562c14faa6ead589577c910ce4 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Thu, 20 Mar 2025 15:04:19 +0100
Subject: [PATCH 05/30] [feature] Plot results in polar coordinates

To compare against external references ; they are all in polar coordinates.
---
 TEMPORARY_FILES/perez_standard_skies.py | 26 ++++++++++++++++++-------
 1 file changed, 19 insertions(+), 7 deletions(-)

diff --git a/TEMPORARY_FILES/perez_standard_skies.py b/TEMPORARY_FILES/perez_standard_skies.py
index 2bd7b69..512fad7 100644
--- a/TEMPORARY_FILES/perez_standard_skies.py
+++ b/TEMPORARY_FILES/perez_standard_skies.py
@@ -142,16 +142,21 @@ def compute_luminance(sky_type, patch_az, patch_el, sun_az, sun_el, zenith_lum):
     v_luminance_fun = np.vectorize(luminance_fun)
     return v_luminance_fun(sky_type, patch_az, patch_el, sun_az, sun_el, zenith_lum)
 
-def plot_luminance_map(az, el, lum, az_s=None, el_s=None, luminance_or_radiance='radiance'):
+def plot_luminance_map(az, el, lum, az_s=None, el_s=None, luminance_or_radiance='radiance', polar=True):
     lum_min = np.min(lum)
     lum_max = np.max(lum)
 
     fig = plt.figure()
-    ax = fig.gca()
+    if polar:
+        ax = fig.add_subplot(111, polar=True)
+        c = ax.pcolor(np.deg2rad(az), 90-el, lum, vmin=lum_min, vmax=lum_max)
+    else:
+        ax = fig.gca()
+
+        c = ax.pcolor(az, el, lum, vmin=lum_min, vmax=lum_max)
 
-    c = ax.pcolor(az, el, lum, vmin=lum_min, vmax=lum_max)
-    ax.set_xlabel('Azimuth [deg]')
-    ax.set_ylabel('Elevation [deg]')
+        ax.set_xlabel('Azimuth [deg]')
+        ax.set_ylabel('Zenith angle [deg]')
 
     cbar = fig.colorbar(c, ax=ax)
     if luminance_or_radiance == 'luminance':
@@ -161,9 +166,16 @@ def plot_luminance_map(az, el, lum, az_s=None, el_s=None, luminance_or_radiance=
     cbar.set_label(cbar_label)
 
     if (az_s is not None) and (el_s is not None):
-        ax.plot(az_s, el_s, 'ro', label='Sun position')
+        if polar:
+            ax.plot(np.deg2rad(az_s), 90-el_s, 'ro', label='Sun position')
+        else:
+            ax.plot(az_s, el_s, 'ro', label='Sun position')
         ax.legend()
 
+    # Set N at the top, E to the right
+    ax.set_theta_zero_location('N')
+    ax.set_theta_direction(-1)
+
     plt.show()
 
 standard_skies = pd.read_csv(os.path.join('..', 'INPUTS', 'CIE_standard_skies.csv'))
@@ -173,7 +185,7 @@ print(compute_gradation(1, 1, 90) == 1)
 # Example from the paper
 sky_type = 11
 
-az = np.arange(0, 180, step=10)  # [°]
+az = np.arange(0, 360, step=5)  # [°]
 el = np.arange(0, 90, step=5)  # [°]
 az_meshgrid, el_meshgrid = np.meshgrid(az, el)
 
-- 
GitLab


From f0cd59d548cdc208c79f6b3381937fd3a2aa771d Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Thu, 20 Mar 2025 15:05:09 +0100
Subject: [PATCH 06/30] [feature] Plot all 15 sky types

In a (5, 3) subplot to have an overview. The results seem consistent when compared to the cie-sky webapp
---
 TEMPORARY_FILES/perez_standard_skies.py | 62 ++++++++++++++++++++++++-
 1 file changed, 61 insertions(+), 1 deletion(-)

diff --git a/TEMPORARY_FILES/perez_standard_skies.py b/TEMPORARY_FILES/perez_standard_skies.py
index 512fad7..59015b5 100644
--- a/TEMPORARY_FILES/perez_standard_skies.py
+++ b/TEMPORARY_FILES/perez_standard_skies.py
@@ -8,8 +8,10 @@ import matplotlib
 matplotlib.use('TkAgg')
 
 from MODULES.conversion_functions import sph_to_cart
+from MODULES.ENVIRONMENT.sky_model import ReinhartSky
 
 logger = logging.getLogger(__name__)
+logger.setLevel('DEBUG')
 
 def compute_sph_angular_dist(Z, Z_s, az, az_s):
     """
@@ -178,6 +180,62 @@ def plot_luminance_map(az, el, lum, az_s=None, el_s=None, luminance_or_radiance=
 
     plt.show()
 
+def plot_all_skies(polar=True, luminance_or_radiance='luminance', colormap='jet'):
+    az = np.arange(0, 360, step=5)  # [°]
+    el = np.arange(0, 90, step=5)  # [°]
+    az_meshgrid, el_meshgrid = np.meshgrid(az, el)
+
+    el_s = 38.02  # [°]
+    az_s = 147.67  # [°] 145° is south-east
+    zenith_luminance = 4404  # cd/m²
+
+    fig, axs = plt.subplots(3, 5,
+                            figsize=(14, 8),
+                            subplot_kw=dict(polar=polar))
+
+    sky_types = np.arange(1, 16)
+    for sky_type in sky_types:
+        lum = compute_luminance(sky_type, az_meshgrid.flatten(), el_meshgrid.flatten(), az_s, el_s, zenith_luminance)
+
+        lum_min = np.min(lum)
+        lum_max = np.max(lum)
+
+        ax_row = int(np.floor((sky_type-1)/5))
+        ax_col = (sky_type-1) % 5
+        print(f'{sky_type=}', f'{ax_row=}', f'{ax_col=}')
+
+        ax = axs[ax_row, ax_col]
+        if polar:
+            c = ax.pcolor(np.deg2rad(az), 90-el, lum.reshape(az_meshgrid.shape), vmin=lum_min, vmax=lum_max, cmap=colormap)
+        else:
+            c = ax.pcolor(az, el, lum.reshape(az_meshgrid.shape), vmin=lum_min, vmax=lum_max, cmap=colormap)
+
+            ax.set_xlabel('Azimuth [deg]')
+            ax.set_ylabel('Zenith [deg]')
+
+        ax.title.set_text(f'Sky type {sky_type}')
+
+        # Set N at the top, E to the right
+        ax.set_theta_zero_location('N')
+        ax.set_theta_direction(-1)
+
+        fig.colorbar(c, ax=ax)
+
+        if (az_s is not None) and (el_s is not None):
+            if polar:
+                ax.plot(np.deg2rad(az_s), 90-el_s, 'ro', markeredgecolor='yellow', label='Sun position')
+            else:
+                ax.plot(az_s, 90-el_s, 'ro', markeredgecolor='yellow', label='Sun position')
+
+    if luminance_or_radiance == 'luminance':
+        fig_title = 'Luminance [cd/m²]'
+    elif luminance_or_radiance == 'radiance':
+        fig_title = 'Radiance [W/(m²⋅sr)]'
+    fig.suptitle(fig_title + ' (azimuth-zenith)')
+
+    plt.tight_layout()
+    plt.show()
+
 standard_skies = pd.read_csv(os.path.join('..', 'INPUTS', 'CIE_standard_skies.csv'))
 
 print(compute_gradation(1, 1, 90) == 1)
@@ -196,4 +254,6 @@ zenith_luminance = 4404  # cd/m²
 lum = compute_luminance(sky_type, az_meshgrid.flatten(), el_meshgrid.flatten(), az_s, el_s, zenith_luminance)
 
 # Plot with pcolor
-plot_luminance_map(az, el, lum.reshape(az_meshgrid.shape), az_s, el_s, luminance_or_radiance='luminance')
+plot_luminance_map(az, el, lum.reshape(az_meshgrid.shape), az_s, el_s, luminance_or_radiance='luminance', polar=True)
+
+plot_all_skies()
-- 
GitLab


From 85e467ef2b317ad9fa8a4c14953a66d154a0306a Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Thu, 20 Mar 2025 16:01:46 +0100
Subject: [PATCH 07/30] [feature] Switch to pcolormesh for efficiency

Matplotlib documentation states pcolormesh is more efficient.
---
 TEMPORARY_FILES/perez_standard_skies.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/TEMPORARY_FILES/perez_standard_skies.py b/TEMPORARY_FILES/perez_standard_skies.py
index 59015b5..89e8fbd 100644
--- a/TEMPORARY_FILES/perez_standard_skies.py
+++ b/TEMPORARY_FILES/perez_standard_skies.py
@@ -151,11 +151,11 @@ def plot_luminance_map(az, el, lum, az_s=None, el_s=None, luminance_or_radiance=
     fig = plt.figure()
     if polar:
         ax = fig.add_subplot(111, polar=True)
-        c = ax.pcolor(np.deg2rad(az), 90-el, lum, vmin=lum_min, vmax=lum_max)
+        c = ax.pcolormesh(np.deg2rad(az), 90-el, lum, vmin=lum_min, vmax=lum_max)
     else:
         ax = fig.gca()
 
-        c = ax.pcolor(az, el, lum, vmin=lum_min, vmax=lum_max)
+        c = ax.pcolormesh(az, el, lum, vmin=lum_min, vmax=lum_max)
 
         ax.set_xlabel('Azimuth [deg]')
         ax.set_ylabel('Zenith angle [deg]')
@@ -206,9 +206,9 @@ def plot_all_skies(polar=True, luminance_or_radiance='luminance', colormap='jet'
 
         ax = axs[ax_row, ax_col]
         if polar:
-            c = ax.pcolor(np.deg2rad(az), 90-el, lum.reshape(az_meshgrid.shape), vmin=lum_min, vmax=lum_max, cmap=colormap)
+            c = ax.pcolormesh(np.deg2rad(az), 90-el, lum.reshape(az_meshgrid.shape), vmin=lum_min, vmax=lum_max, cmap=colormap)
         else:
-            c = ax.pcolor(az, el, lum.reshape(az_meshgrid.shape), vmin=lum_min, vmax=lum_max, cmap=colormap)
+            c = ax.pcolormesh(az, el, lum.reshape(az_meshgrid.shape), vmin=lum_min, vmax=lum_max, cmap=colormap)
 
             ax.set_xlabel('Azimuth [deg]')
             ax.set_ylabel('Zenith [deg]')
@@ -253,7 +253,7 @@ zenith_luminance = 4404  # cd/m²
 
 lum = compute_luminance(sky_type, az_meshgrid.flatten(), el_meshgrid.flatten(), az_s, el_s, zenith_luminance)
 
-# Plot with pcolor
+# Plot with pcolormesh
 plot_luminance_map(az, el, lum.reshape(az_meshgrid.shape), az_s, el_s, luminance_or_radiance='luminance', polar=True)
 
 plot_all_skies()
-- 
GitLab


From 2266f9e91e665cb93ee009c69c9a19245244b7d2 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Tue, 25 Mar 2025 12:54:56 +0100
Subject: [PATCH 08/30] [refactor] Refactor Reinhart sky discretization

There is a pd.concat() statement in a for loop, which is inefficient.

Solution: in the for loop, append dictionaries to a list, then after the loop, convert the list of dictionaries into a pandas dataframe. It reduces computing time by a factor of 10 on my computer.
---
 MODULES/ENVIRONMENT/sky_model.py | 23 +++++++++++++----------
 1 file changed, 13 insertions(+), 10 deletions(-)

diff --git a/MODULES/ENVIRONMENT/sky_model.py b/MODULES/ENVIRONMENT/sky_model.py
index 04da300..64639ac 100644
--- a/MODULES/ENVIRONMENT/sky_model.py
+++ b/MODULES/ENVIRONMENT/sky_model.py
@@ -112,12 +112,14 @@ class ReinhartSky:
         return reinhart_sky, reinhart_num_total
 
     def build_patches(self):
-        self.reinhart_patches = pd.DataFrame()
         self.patch_id = 0
         self.cols = ['row', 'patch_id', 'patch_id_in_row', 'd_el', 'd_az', 'el',
                      'az', 'solid_angle_sr']
+
+        # temporary list that will be appended to (dictionaries) in the loop,
+        # then converted to Dataframe after the loop
+        temp_list = []
         for row in self.reinhart_sky['row']:
-            # print(f'{row=}')
             d_az = 360 / \
                    self.reinhart_sky[self.reinhart_sky['row'] == row][
                        'num_segments'].values[
@@ -125,21 +127,22 @@ class ReinhartSky:
             el = self.alpha * (row + 1 / 2)  # [deg]
 
             for segment in range(self.reinhart_sky['num_segments'].iloc[row]):
-                # print(f'{segment=}')
                 az = d_az * (segment + 1 / 2)  # [deg]
                 solid_angle = self.compute_solid_angle_from_elev_azim(el,
-                                                                      d_az)
+                                                                      d_az)  # [sr]
 
-                df = pd.DataFrame(
-                    [[row, self.patch_id, segment, self.alpha, d_az, el, az,
-                      solid_angle]],
-                    columns=self.cols)
+                # Dictionary that will be appended to temp_list
+                dict = {'row': row, 'patch_id': self.patch_id, 'patch_id_in_row': segment,
+                        'd_el': self.alpha, 'd_az': d_az, 'el': el,
+                     'az': az, 'solid_angle_sr': solid_angle}
 
-                self.reinhart_patches = pd.concat([self.reinhart_patches, df],
-                                             ignore_index=True)
+                temp_list.append(dict)
 
                 self.patch_id += 1
 
+            # Convert the list of dictionaries to DataFrame
+            self.reinhart_patches = pd.DataFrame.from_records(temp_list)
+
     def add_top_patch(self):
         top_row = max(self.reinhart_sky['row']) + 1
         top_patch_solid_angle = self.compute_solid_angle_cone(self.alpha / 2)
-- 
GitLab


From 5db58ca4369e95c20f014a6509636b3d8787e011 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Tue, 25 Mar 2025 12:57:32 +0100
Subject: [PATCH 09/30] [feature] Adjust tolerance when checking sum of solid
 angles in Reinhart sky discretization

The tolerance was very small (default value, rel_tol = 1e-9) which almost always raised an error.

Set rel_tol=5e-3 to give some leeway.
---
 MODULES/ENVIRONMENT/sky_model.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/MODULES/ENVIRONMENT/sky_model.py b/MODULES/ENVIRONMENT/sky_model.py
index 64639ac..18e1d1f 100644
--- a/MODULES/ENVIRONMENT/sky_model.py
+++ b/MODULES/ENVIRONMENT/sky_model.py
@@ -205,9 +205,9 @@ class ReinhartSky:
         sum_solid_angles = self.reinhart_patches['solid_angle_sr'].sum()
         sum_solid_angles_div_pi = sum_solid_angles / np.pi
         logger.debug(f'Sum of solid angles = {sum_solid_angles_div_pi:.4g} * pi')
-        if not math.isclose(sum_solid_angles, 2*np.pi):
+        if not math.isclose(sum_solid_angles, 2*np.pi, rel_tol=5e-3):
             logger.warning('Sum of solid angles seems off, check the sky discretization scheme')
-            logger.warning(f'Sum of solid angles = {sum_solid_angles_div_pi:.4g} * pi')
+            logger.warning(f'Sum of solid angles = {sum_solid_angles_div_pi:.6g} * pi')
 
     def get_patches_xyz(self):
         # Compute xyz coords on a unit sphere (actually half sphere since it's the sky dome)
-- 
GitLab


From 764ce4c0c3b07fe0fb5521d7d7545072f8542448 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Tue, 25 Mar 2025 14:06:32 +0100
Subject: [PATCH 10/30] [perf] Implement profiling with cProfile

cProfile is a module of the standard Python library and provides stats on function calls within a module or script.
---
 .gitignore                             |  8 ++++++--
 MODULES/ENVIRONMENT/sky_model.py       | 24 +++++++++++++++++++-----
 TEMPORARY_FILES/parse_cprofile_data.py | 18 ++++++++++++++++++
 3 files changed, 43 insertions(+), 7 deletions(-)
 create mode 100644 TEMPORARY_FILES/parse_cprofile_data.py

diff --git a/.gitignore b/.gitignore
index f93affb..37f53ca 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,13 +4,17 @@
 
 *.png
 
+# PyCharm project folder
 .idea/
 
 *.csv
 
-
+# STICS related files
 *.sti
 *param_files/*
 *default/*
 *logs/*
-*STICS/*
\ No newline at end of file
+*STICS/*
+
+# Profiling outputs
+*.prof
\ No newline at end of file
diff --git a/MODULES/ENVIRONMENT/sky_model.py b/MODULES/ENVIRONMENT/sky_model.py
index 18e1d1f..28f0e16 100644
--- a/MODULES/ENVIRONMENT/sky_model.py
+++ b/MODULES/ENVIRONMENT/sky_model.py
@@ -5,9 +5,11 @@ Paper openly available on Researchgate: https://www.google.com/url?sa=t&source=w
 
 """
 
+import cProfile
 import logging
 import math
 import numpy as np
+import os
 import pandas as pd
 from matplotlib import pyplot as plt
 import matplotlib
@@ -266,8 +268,20 @@ class ReinhartSky:
 
 
 if __name__ == '__main__':
-    sky = ReinhartSky(MF=1)
-    sky.scatter_2D()
-    sky.scatter_2D(text_id=True)
-    sky.scatter_3D()
-    sky.scatter_3D(text_id=True)
+    MF = 32
+    profile_flag = True
+    MFs = [1, 2]
+    for MF in MFs:
+        if profile_flag:
+            profile_output = os.path.join('..', '..', 'OUTPUTS', 'profiling', f'cprofile_output_sky_model_MF{MF}.prof')
+            cProfile.run(f"ReinhartSky(MF={MF})", profile_output, sort='tottime')
+        else:
+            sky = ReinhartSky(MF=MF)
+            if MF == 1:
+                text_id = True
+            else:
+                text_id = False
+
+            if MF < 8:
+                sky.scatter_2D(text_id=text_id)
+                sky.scatter_3D(text_id=text_id)
diff --git a/TEMPORARY_FILES/parse_cprofile_data.py b/TEMPORARY_FILES/parse_cprofile_data.py
new file mode 100644
index 0000000..6ccc77a
--- /dev/null
+++ b/TEMPORARY_FILES/parse_cprofile_data.py
@@ -0,0 +1,18 @@
+import os
+import pstats
+from pstats import SortKey
+
+from MODULES.ENVIRONMENT.sky_model import ReinhartSky
+
+# Profiles of sky_model Reinhart module
+MF = 8
+MFs = [16, 32]
+for MF in MFs:
+    profile_output = os.path.join('..', 'OUTPUTS', 'profiling', f'cprofile_output_sky_model_MF{MF}.prof')
+    p = pstats.Stats(profile_output)
+    p.strip_dirs().sort_stats('tottime').print_stats(10)
+
+# Profile of example.py
+profile_output = os.path.join('..', 'OUTPUTS', 'profiling', f'example.prof')
+p = pstats.Stats(profile_output)
+p.strip_dirs().sort_stats('tottime').print_stats(100)
-- 
GitLab


From 62624f43661b2141ed71ed288f89652c3b34ef3e Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Tue, 25 Mar 2025 14:39:21 +0100
Subject: [PATCH 11/30] [feature] Return relative luminance instead of absolute
 in luminance_fun()

It is more useful considering the integration in the PASE modules
---
 TEMPORARY_FILES/perez_standard_skies.py | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/TEMPORARY_FILES/perez_standard_skies.py b/TEMPORARY_FILES/perez_standard_skies.py
index 89e8fbd..670aa6a 100644
--- a/TEMPORARY_FILES/perez_standard_skies.py
+++ b/TEMPORARY_FILES/perez_standard_skies.py
@@ -135,10 +135,7 @@ def luminance_fun(sky_type, patch_az, patch_el, sun_az, sun_el, zenith_lum):
     # Relative luminance
     rel_lum = phi * indicatrix / (phi_zenith * indicatrix_zenith)
 
-    # Absolute luminance
-    lum = rel_lum * zenith_lum
-
-    return lum
+    return rel_lum
 
 def compute_luminance(sky_type, patch_az, patch_el, sun_az, sun_el, zenith_lum):
     v_luminance_fun = np.vectorize(luminance_fun)
-- 
GitLab


From 6a61ca9449f18791d7f4faf5ce1654385d03130a Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Tue, 25 Mar 2025 14:41:37 +0100
Subject: [PATCH 12/30] [refactor] Place script part in a __name__ guard

Improves the layout of the file
---
 TEMPORARY_FILES/perez_standard_skies.py | 26 ++++++++++++++-----------
 1 file changed, 15 insertions(+), 11 deletions(-)

diff --git a/TEMPORARY_FILES/perez_standard_skies.py b/TEMPORARY_FILES/perez_standard_skies.py
index 670aa6a..0026d1b 100644
--- a/TEMPORARY_FILES/perez_standard_skies.py
+++ b/TEMPORARY_FILES/perez_standard_skies.py
@@ -233,22 +233,26 @@ def plot_all_skies(polar=True, luminance_or_radiance='luminance', colormap='jet'
     plt.tight_layout()
     plt.show()
 
-standard_skies = pd.read_csv(os.path.join('..', 'INPUTS', 'CIE_standard_skies.csv'))
+if __name__ == '__main__':
 
-print(compute_gradation(1, 1, 90) == 1)
+    standard_skies = pd.read_csv(os.path.join('..', 'INPUTS', 'CIE_standard_skies.csv'))
 
-# Example from the paper
-sky_type = 11
+    # Example from the paper
+    sky_type = 2
 
-az = np.arange(0, 360, step=5)  # [°]
-el = np.arange(0, 90, step=5)  # [°]
-az_meshgrid, el_meshgrid = np.meshgrid(az, el)
+    # az = np.arange(0, 360, step=5)  # [°]
+    # el = np.arange(0, 90, step=5)  # [°]
+    sky = ReinhartSky(MF=2).reinhart_patches
+    az = sky['az']
+    el = sky['el']
+    az_meshgrid, el_meshgrid = np.meshgrid(az, el)
 
-el_s = 38.02  # [°]
-az_s = 147.67  # [°] 145° is south-east
-zenith_luminance = 4404  # cd/m²
+    el_s = 38.02  # [°]
+    az_s = 147.67  # [°] 145° is south-east
+    zenith_luminance = 4404  # cd/m²
 
-lum = compute_luminance(sky_type, az_meshgrid.flatten(), el_meshgrid.flatten(), az_s, el_s, zenith_luminance)
+    rel_lum = compute_luminance(sky_type, az_meshgrid.flatten(), el_meshgrid.flatten(), az_s, el_s, zenith_luminance)
+    lum = rel_lum * zenith_luminance
 
 # Plot with pcolormesh
 plot_luminance_map(az, el, lum.reshape(az_meshgrid.shape), az_s, el_s, luminance_or_radiance='luminance', polar=True)
-- 
GitLab


From 3f785edaa4153a91bc228a017f0f33779e3e4611 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Tue, 25 Mar 2025 14:47:47 +0100
Subject: [PATCH 13/30] [feature] Plot 3D with PyVista

Plot the sky dome in 3D and the relative luminance/radiance as a colorscale on the half-sphere. The sun direction is represented as a small yellow sphere.
---
 TEMPORARY_FILES/perez_standard_skies.py | 61 +++++++++++++++++++++++--
 1 file changed, 56 insertions(+), 5 deletions(-)

diff --git a/TEMPORARY_FILES/perez_standard_skies.py b/TEMPORARY_FILES/perez_standard_skies.py
index 0026d1b..bbab7f8 100644
--- a/TEMPORARY_FILES/perez_standard_skies.py
+++ b/TEMPORARY_FILES/perez_standard_skies.py
@@ -3,6 +3,7 @@ import math
 import numpy as np
 import os
 import pandas as pd
+import pyvista as pv
 from matplotlib import pyplot as plt
 import matplotlib
 matplotlib.use('TkAgg')
@@ -108,7 +109,7 @@ def get_sky_params(standard_sky):
 
     return a, b, c, d, e
 
-def luminance_fun(sky_type, patch_az, patch_el, sun_az, sun_el, zenith_lum):
+def luminance_fun(sky_type, patch_az, patch_el, sun_az, sun_el):
     # Standard sky parameters and description
     standard_sky = standard_skies.loc[standard_skies['Type'] == sky_type]
 
@@ -178,8 +179,12 @@ def plot_luminance_map(az, el, lum, az_s=None, el_s=None, luminance_or_radiance=
     plt.show()
 
 def plot_all_skies(polar=True, luminance_or_radiance='luminance', colormap='jet'):
+    # sky = ReinhartSky().reinhart_patches
+    # az = sky['az']
+    # el = sky['el']
     az = np.arange(0, 360, step=5)  # [°]
     el = np.arange(0, 90, step=5)  # [°]
+
     az_meshgrid, el_meshgrid = np.meshgrid(az, el)
 
     el_s = 38.02  # [°]
@@ -233,6 +238,47 @@ def plot_all_skies(polar=True, luminance_or_radiance='luminance', colormap='jet'
     plt.tight_layout()
     plt.show()
 
+def plot_3D(sky, rel_lum):
+
+    sky['patch_bound_el_low'] = sky['el'] - sky['d_el'] / 2
+    sky['patch_bound_el_up'] = sky['el'] + sky['d_el'] / 2
+    sky['patch_bound_az_low'] = sky['az'] - sky['d_az'] / 2
+    sky['patch_bound_az_up'] = sky['az'] + sky['d_az'] / 2
+    print(sky['patch_bound_el_up'].iloc[-1])
+    # Fix elevation for top patch
+    sky['patch_bound_el_up'].iloc[-1] = 90
+    print(sky['patch_bound_el_up'].iloc[len(sky) - 1])
+
+    # Vertical levels
+    # in this case a single level slightly above the surface of a sphere
+    RADIUS = 1
+    levels = [RADIUS * 1.01]
+
+    xx_bounds = np.concatenate([sky['patch_bound_az_low'].values,
+                                [sky['patch_bound_az_up'].values[
+                                     -1]]])  # az
+    xx_bounds[-2:-1] += 180
+    yy_bounds = np.concatenate([sky['patch_bound_el_low'].values,
+                                [sky['patch_bound_el_up'].values[
+                                     -1]]])  # el
+
+    grid_scalar = pv.grid_from_sph_coords(xx_bounds, 90 - yy_bounds, levels)
+
+    # And fill its cell arrays with the scalar data
+    grid_scalar.cell_data["Luminance relative to zenith"] = np.array(
+        rel_lum.reshape(az_meshgrid.shape)).swapaxes(-2, -1).ravel("C")
+
+    # Prepare the Sun
+    xyz_sun = sph_to_cart('deg', azimut=az_s, elev=el_s, zenith_angle=None,
+                          dist=RADIUS)
+    sun_mesh = pv.Sphere(radius=RADIUS / 20, center=xyz_sun)
+
+    # Make a plot
+    p = pv.Plotter()
+    p.add_mesh(sun_mesh, color='yellow')
+    p.add_mesh(grid_scalar, opacity=0.8, cmap="jet")
+    p.show()
+
 if __name__ == '__main__':
 
     standard_skies = pd.read_csv(os.path.join('..', 'INPUTS', 'CIE_standard_skies.csv'))
@@ -242,7 +288,7 @@ if __name__ == '__main__':
 
     # az = np.arange(0, 360, step=5)  # [°]
     # el = np.arange(0, 90, step=5)  # [°]
-    sky = ReinhartSky(MF=2).reinhart_patches
+    sky = ReinhartSky(MF=1).reinhart_patches
     az = sky['az']
     el = sky['el']
     az_meshgrid, el_meshgrid = np.meshgrid(az, el)
@@ -254,7 +300,12 @@ if __name__ == '__main__':
     rel_lum = compute_luminance(sky_type, az_meshgrid.flatten(), el_meshgrid.flatten(), az_s, el_s, zenith_luminance)
     lum = rel_lum * zenith_luminance
 
-# Plot with pcolormesh
-plot_luminance_map(az, el, lum.reshape(az_meshgrid.shape), az_s, el_s, luminance_or_radiance='luminance', polar=True)
+    # Plot with pcolormesh
+    plot_luminance_map(az, el, lum.reshape(az_meshgrid.shape), az_s, el_s, luminance_or_radiance='luminance', polar=True)
+
+    # Plot all 15 standard skies in a (5, 3) subplots figure (takes a couple of minutes)
+    # plot_all_skies()
+
+    # Plot with PyVista
+    plot_3D(sky, rel_lum)
 
-plot_all_skies()
-- 
GitLab


From 731a0aea558dff4e26a83e5f4a7bebc18f3e5af4 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Fri, 28 Mar 2025 09:01:01 +0100
Subject: [PATCH 14/30] [refactor] Move fibonacci_half_sphere function from
 light.py to sky_model.py

Regroup the function in the same module as ReinhartSky discretization class
---
 MODULES/ENVIRONMENT/light.py     | 21 +--------------------
 MODULES/ENVIRONMENT/sky_model.py | 21 +++++++++++++++++++++
 2 files changed, 22 insertions(+), 20 deletions(-)

diff --git a/MODULES/ENVIRONMENT/light.py b/MODULES/ENVIRONMENT/light.py
index c09c7a2..0c6ee9e 100644
--- a/MODULES/ENVIRONMENT/light.py
+++ b/MODULES/ENVIRONMENT/light.py
@@ -14,29 +14,10 @@ import matplotlib.pyplot as plt
 import os
 
 import MODULES.conversion_functions as cf
-from MODULES.ENVIRONMENT.sky_model import ReinhartSky
+from MODULES.ENVIRONMENT.sky_model import ReinhartSky, fibonacci_half_sphere
 
 logger = logging.getLogger(__name__)
 
-def fibonacci_half_sphere(samples=18):
-    """
-    Function computing a number of direction sampling a virtual upper half sphere with a center at 0,0,0
-    
-    Parameters:
-        samples (int): Number of directions
-
-    Returns:
-       Directions (np.array): Matrix (samples x 3) providing the direction
-    """
-    phi = np.pi * (3. - np.sqrt(5.))
-    i = np.linspace(0,samples-1,num=samples)
-    yp = (1 - i/float(samples-1))
-    radius = np.sqrt(1-yp**2) 
-    theta = phi * i 
-    xp = np.cos(theta) * radius
-    zp = np.sin(theta) * radius
-    return np.column_stack([xp,zp,yp])
-
 
 class Sun_positions:
     
diff --git a/MODULES/ENVIRONMENT/sky_model.py b/MODULES/ENVIRONMENT/sky_model.py
index 28f0e16..2bb95fa 100644
--- a/MODULES/ENVIRONMENT/sky_model.py
+++ b/MODULES/ENVIRONMENT/sky_model.py
@@ -19,6 +19,27 @@ from MODULES.conversion_functions import sph_to_cart
 
 logger = logging.getLogger(__name__)
 
+
+def fibonacci_half_sphere(samples=18):
+    """
+    Function computing a number of direction sampling a virtual upper half sphere with a center at 0,0,0
+
+    Parameters:
+        samples (int): Number of directions
+
+    Returns:
+       Directions (np.array): Matrix (samples x 3) providing the direction
+    """
+    phi = np.pi * (3. - np.sqrt(5.))
+    i = np.linspace(0, samples - 1, num=samples)
+    yp = (1 - i / float(samples - 1))
+    radius = np.sqrt(1 - yp ** 2)
+    theta = phi * i
+    xp = np.cos(theta) * radius
+    zp = np.sin(theta) * radius
+    return np.column_stack([xp, zp, yp])
+
+
 class ReinhartSky:
 
     def __init__(self, MF=1):
-- 
GitLab


From 4c573f54dcf0cafc116ac448f7dfd23a386ff0dc Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Mon, 7 Apr 2025 15:13:50 +0200
Subject: [PATCH 15/30] [fix] Complete the switch from absolute to relative
 luminance

I forgot some pieces of code in commit 62624f43661b2141ed71ed288f89652c3b34ef3e
---
 TEMPORARY_FILES/perez_standard_skies.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/TEMPORARY_FILES/perez_standard_skies.py b/TEMPORARY_FILES/perez_standard_skies.py
index bbab7f8..716df8e 100644
--- a/TEMPORARY_FILES/perez_standard_skies.py
+++ b/TEMPORARY_FILES/perez_standard_skies.py
@@ -138,9 +138,9 @@ def luminance_fun(sky_type, patch_az, patch_el, sun_az, sun_el):
 
     return rel_lum
 
-def compute_luminance(sky_type, patch_az, patch_el, sun_az, sun_el, zenith_lum):
+def compute_luminance(sky_type, patch_az, patch_el, sun_az, sun_el):
     v_luminance_fun = np.vectorize(luminance_fun)
-    return v_luminance_fun(sky_type, patch_az, patch_el, sun_az, sun_el, zenith_lum)
+    return v_luminance_fun(sky_type, patch_az, patch_el, sun_az, sun_el)
 
 def plot_luminance_map(az, el, lum, az_s=None, el_s=None, luminance_or_radiance='radiance', polar=True):
     lum_min = np.min(lum)
@@ -197,7 +197,7 @@ def plot_all_skies(polar=True, luminance_or_radiance='luminance', colormap='jet'
 
     sky_types = np.arange(1, 16)
     for sky_type in sky_types:
-        lum = compute_luminance(sky_type, az_meshgrid.flatten(), el_meshgrid.flatten(), az_s, el_s, zenith_luminance)
+        lum = compute_luminance(sky_type, az_meshgrid.flatten(), el_meshgrid.flatten(), az_s, el_s)
 
         lum_min = np.min(lum)
         lum_max = np.max(lum)
-- 
GitLab


From b2f133b9b7e98aac8561f29407e05dc4f842d486 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Mon, 7 Apr 2025 16:12:48 +0200
Subject: [PATCH 16/30] [feature] CIEStandardSky class

Takes as input the discretized sky, sun position angles (azimuth and elevation) and sky type to compute the radiance distribution relative to the value at the zenith.

There are currently no plotting capabilities because the pcolor function that I would use is not plug-and-play with the sky discretization scheme.
---
 MODULES/ENVIRONMENT/sky_model.py | 178 +++++++++++++++++++++++++++++++
 1 file changed, 178 insertions(+)

diff --git a/MODULES/ENVIRONMENT/sky_model.py b/MODULES/ENVIRONMENT/sky_model.py
index 2bb95fa..dd6848e 100644
--- a/MODULES/ENVIRONMENT/sky_model.py
+++ b/MODULES/ENVIRONMENT/sky_model.py
@@ -288,6 +288,184 @@ class ReinhartSky:
         plt.show()
 
 
+class CIEStandardSky:
+    """
+    Generate a standard radiance distribution among the 15 CIE General Skies.
+    """
+
+    def __init__(self, discrete_sky,
+                 sun_az, sun_el,
+                 sky_type=5):
+
+        self.sky = discrete_sky
+
+        self.az = discrete_sky['az'].values
+        self.el = discrete_sky['el'].values
+
+        self.az_s = sun_az  # [°] Sun azimuth (145° is south-east)
+        self.el_s = sun_el  # [°]  Sun elevation
+
+        self.sky_type = sky_type
+
+        self.rel_radiance_distribution = self.compute_rel_radiance()
+
+
+    def compute_sph_angular_dist(self, Z, Z_s, az, az_s):
+        """
+
+        :param Z: Sky patch zenith angle [deg]
+        :param Z: float
+        :param Z_s: Solar zenith angle [deg]
+        :param Z_s: float
+        :param az: Sky patch azimuth angle from north [deg]
+        :param az: float
+        :param az_s: Solar azimuth angle from north [deg]
+        :param az_s: float
+        :return: chi = Spherical angular distance [rad]
+        """
+
+        # Convert all angles from degrees to radians
+        Z = np.deg2rad(Z)
+        Z_s = np.deg2rad(Z_s)
+        az = np.deg2rad(az)
+        az_s = np.deg2rad(az_s)
+
+        # Absolute difference in azimuth angles
+        Az = np.abs(az - az_s)
+
+        # Angular distance chi
+        chi = np.arccos(
+            np.cos(Z_s) * np.cos(Z) + np.sin(Z_s) * np.sin(Z) * np.cos(Az))
+
+        return chi
+
+    def gradation_fun(self, a, b, Z):
+        """
+        Gradation function. It is not used as is but vectorized in compute_gradation().
+
+        :param a: Horizon-zenith gradient (gradation parameter) [-]
+        :type a: float
+        :param b: Gradient intensity (gradation parameter) [-]
+        :type b: float
+        :param Z: Zenith angle [deg]
+        :type Z: float
+        :return: gradation function value [-]
+        """
+        if math.isclose(Z, np.pi / 2):
+            return 1
+        else:
+            return 1 + a * np.exp(b / np.cos(Z))
+
+    def compute_gradation(self, a, b, Z):
+        """
+        Vectorizes and computes the gradation function
+
+        :param a: Horizon-zenith gradient (gradation parameter) [-]
+        :type a: float
+        :param b: Gradient intensity (gradation parameter) [-]
+        :type b: float
+        :param Z: Zenith angle [deg]
+        :type Z: Nx1 array of float
+        :return: gradation function [-], same shape as Z
+        """
+
+        v_gradation = np.vectorize(self.gradation_fun)
+
+        Z = np.deg2rad(Z)
+
+        return v_gradation(a, b, Z)
+
+    def scattering_indicatrix_fun(self, c, d, e, x):
+        """
+        Compute the scattering indicatrix value.
+
+        :param c: Circumsolar intensity (scattering parameter) [-]
+        :type c: float
+        :param d: Circumsolar radius (scattering parameter) [-]
+        :type d: float
+        :param e: Backscattering effect (scattering parameter) [-]
+        :type e: float
+        :param x: angle (spherical angular distance or zenith, i.e. 0) [rad]
+        :type x: array of float or float
+        :return: Scattering indicatrix value [-], same shape as x
+        """
+        return 1 + c * (np.exp(d * x) - np.exp(d * np.pi / 2)) + e * (
+            np.cos(x)) ** 2
+
+    def get_sky_params(self, standard_sky):
+        """
+
+        :param standard_sky: Dataframe of length 1 with standard sky information
+        :return a: Horizon-zenith gradient (gradation parameter) [-]
+        :return b: Gradient intensity (gradation parameter) [-]
+        :return c: Circumsolar intensity (scattering parameter) [-]
+        :return d: Circumsolar radius (scattering parameter) [-]
+        :return e: Backscattering effect (scattering parameter) [-]
+        """
+
+        a = standard_sky['a'].values[0]
+        b = standard_sky['b'].values[0]
+        c = standard_sky['c'].values[0]
+        d = standard_sky['d'].values[0]
+        e = standard_sky['e'].values[0]
+
+        return a, b, c, d, e
+
+    def relative_radiance_fun(self, sky_type=None, az=None, el=None):
+        """
+
+        :return: relative quantity (radiance or luminance) with respect to
+                 value at the zenith
+        """
+
+        # Read the csv containing the descriptions of the 15 standard skies
+        standard_skies = pd.read_csv(os.path.join('..', '..',
+                                                  'INPUTS',
+                                                  'CIE_standard_skies.csv'))
+
+        # Extract the parameters of the desired sky type
+        if sky_type is None:
+            standard_sky = standard_skies.loc[standard_skies['Type'] == self.sky_type]
+        else:
+            standard_sky = standard_skies.loc[standard_skies['Type'] == sky_type]
+
+        # Sky model parameters
+        a, b, c, d, e = self.get_sky_params(standard_sky)
+
+        # Elevation to zenith (Sun)
+        Z_s = 90 - self.el_s  # [deg]
+
+        # Elevation to zenith (sky patch)
+        if el is None:
+            Z = 90 - self.el  # [deg]
+        else:
+            Z = 90 - el  # [deg]
+
+
+        # Spherical angular dist
+        if az is None:
+            chi = self.compute_sph_angular_dist(Z, Z_s, self.az, self.az_s)
+        else:
+            chi = self.compute_sph_angular_dist(Z, Z_s, az, self.az_s)
+
+        # Gradation
+        phi = self.compute_gradation(a, b, Z)
+        phi_zenith = self.compute_gradation(a, b, Z=0)  # TODO do not recompute each time !
+
+        # Scattering indicatrix
+        indicatrix = self.scattering_indicatrix_fun(c, d, e, chi)
+        indicatrix_zenith = self.scattering_indicatrix_fun(c, d, e, np.deg2rad(Z_s))  # TODO do not recompute each time !
+
+        # Relative luminance/radiance (called "rel_quantity" for genericity)
+        rel_quantity = phi * indicatrix / (phi_zenith * indicatrix_zenith)
+
+        return rel_quantity
+
+    def compute_rel_radiance(self, sky_type=None, az=None, el=None):
+        v_radiance_fun = np.vectorize(self.relative_radiance_fun)
+        return v_radiance_fun(sky_type, az, el)
+
+
 if __name__ == '__main__':
     MF = 32
     profile_flag = True
-- 
GitLab


From 9fe25cfcb1f91c0b86d88e8dfe1a5df5e8f5f526 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Mon, 7 Apr 2025 16:14:36 +0200
Subject: [PATCH 17/30] [feature] Add the call to CIEStandardSky class

and also some minor refactor of the script part of the file
---
 MODULES/ENVIRONMENT/sky_model.py | 34 +++++++++++++++++++-------------
 1 file changed, 20 insertions(+), 14 deletions(-)

diff --git a/MODULES/ENVIRONMENT/sky_model.py b/MODULES/ENVIRONMENT/sky_model.py
index dd6848e..05f601c 100644
--- a/MODULES/ENVIRONMENT/sky_model.py
+++ b/MODULES/ENVIRONMENT/sky_model.py
@@ -467,20 +467,26 @@ class CIEStandardSky:
 
 
 if __name__ == '__main__':
-    MF = 32
-    profile_flag = True
-    MFs = [1, 2]
-    for MF in MFs:
-        if profile_flag:
+
+    profile_flag = False
+
+    if profile_flag:
+        MFs = [1, 2]
+        for MF in MFs:
             profile_output = os.path.join('..', '..', 'OUTPUTS', 'profiling', f'cprofile_output_sky_model_MF{MF}.prof')
             cProfile.run(f"ReinhartSky(MF={MF})", profile_output, sort='tottime')
+    else:
+        MF = 1
+
+        sky = ReinhartSky(MF=MF)
+        if MF == 1:
+            text_id = True
         else:
-            sky = ReinhartSky(MF=MF)
-            if MF == 1:
-                text_id = True
-            else:
-                text_id = False
-
-            if MF < 8:
-                sky.scatter_2D(text_id=text_id)
-                sky.scatter_3D(text_id=text_id)
+            text_id = False
+
+        if MF < 8:
+            sky.scatter_2D(text_id=text_id)
+            sky.scatter_3D(text_id=text_id)
+
+    sky.reinhart_patches['Relative radiance distribution'] = CIEStandardSky(sky.reinhart_patches,
+                               38.02, 147.67, 15).rel_radiance_distribution
-- 
GitLab


From 10e347bc93ad075525f7f129565000b66ca6c75f Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Mon, 7 Apr 2025 16:25:40 +0200
Subject: [PATCH 18/30] [feature] Script used to generate figures of the
 Standard skies

A bit messy, but it's not really in the code base. It is only intended to be used to generate figures if necessary.
---
 TEMPORARY_FILES/perez_standard_skies.py | 71 +++++++++++++++++++------
 1 file changed, 54 insertions(+), 17 deletions(-)

diff --git a/TEMPORARY_FILES/perez_standard_skies.py b/TEMPORARY_FILES/perez_standard_skies.py
index 716df8e..7937019 100644
--- a/TEMPORARY_FILES/perez_standard_skies.py
+++ b/TEMPORARY_FILES/perez_standard_skies.py
@@ -4,9 +4,9 @@ import numpy as np
 import os
 import pandas as pd
 import pyvista as pv
-from matplotlib import pyplot as plt
 import matplotlib
 matplotlib.use('TkAgg')
+from matplotlib import pyplot as plt
 
 from MODULES.conversion_functions import sph_to_cart
 from MODULES.ENVIRONMENT.sky_model import ReinhartSky
@@ -142,18 +142,25 @@ def compute_luminance(sky_type, patch_az, patch_el, sun_az, sun_el):
     v_luminance_fun = np.vectorize(luminance_fun)
     return v_luminance_fun(sky_type, patch_az, patch_el, sun_az, sun_el)
 
-def plot_luminance_map(az, el, lum, az_s=None, el_s=None, luminance_or_radiance='radiance', polar=True):
+def print_n_max_values(vect, n):
+    a = vect.copy()
+    b = np.unique(a)
+    c = np.sort(b)
+    c = c[::-1]
+    print(c[:n])
+
+def plot_luminance_map(az, el, lum, az_s=None, el_s=None, luminance_or_radiance='radiance', polar=True, colormap=None):
     lum_min = np.min(lum)
     lum_max = np.max(lum)
 
     fig = plt.figure()
     if polar:
         ax = fig.add_subplot(111, polar=True)
-        c = ax.pcolormesh(np.deg2rad(az), 90-el, lum, vmin=lum_min, vmax=lum_max)
+        c = ax.pcolormesh(np.deg2rad(az), 90-el, lum, vmin=lum_min, vmax=lum_max, cmap=colormap)
     else:
         ax = fig.gca()
 
-        c = ax.pcolormesh(az, el, lum, vmin=lum_min, vmax=lum_max)
+        c = ax.pcolormesh(az, el, lum, vmin=lum_min, vmax=lum_max, cmap=colormap)
 
         ax.set_xlabel('Azimuth [deg]')
         ax.set_ylabel('Zenith angle [deg]')
@@ -172,9 +179,10 @@ def plot_luminance_map(az, el, lum, az_s=None, el_s=None, luminance_or_radiance=
             ax.plot(az_s, el_s, 'ro', label='Sun position')
         ax.legend()
 
-    # Set N at the top, E to the right
-    ax.set_theta_zero_location('N')
-    ax.set_theta_direction(-1)
+    if polar:
+        # Set N at the top, E to the right
+        ax.set_theta_zero_location('N')
+        ax.set_theta_direction(-1)
 
     plt.show()
 
@@ -284,28 +292,57 @@ if __name__ == '__main__':
     standard_skies = pd.read_csv(os.path.join('..', 'INPUTS', 'CIE_standard_skies.csv'))
 
     # Example from the paper
-    sky_type = 2
+    sky_type = 12
 
-    # az = np.arange(0, 360, step=5)  # [°]
-    # el = np.arange(0, 90, step=5)  # [°]
-    sky = ReinhartSky(MF=1).reinhart_patches
-    az = sky['az']
-    el = sky['el']
-    az_meshgrid, el_meshgrid = np.meshgrid(az, el)
+    sky_mesh = 'arange'  # {'arange', 'linspace', 'Reinhart'}
+
+    if sky_mesh.lower() == 'reinhart':
+        sky = ReinhartSky(MF=1).reinhart_patches
+        az = sky['az'].values
+        el = sky['el'].values
+        # meshgrid_flag = False
+        meshgrid_flag = True
+        az_meshgrid, el_meshgrid = np.meshgrid(az, el)
+    else:
+        if sky_mesh.lower() == 'arange':
+            step = 5
+            az = np.arange(0, 360+step, step=step)  # [°]
+            el = np.arange(0, 90+step, step=step)  # [°]
+        elif sky_mesh.lower() == 'linspace':
+            az = np.linspace(0, 360, num=145)  # [°]
+            el = np.linspace(0, 90, num=145)  # [°]
+
+        else:
+            raise ValueError('Unrecognized sky discretization')
+
+        meshgrid_flag = True
+        az_meshgrid, el_meshgrid = np.meshgrid(az, el)
 
     el_s = 38.02  # [°]
     az_s = 147.67  # [°] 145° is south-east
     zenith_luminance = 4404  # cd/m²
 
-    rel_lum = compute_luminance(sky_type, az_meshgrid.flatten(), el_meshgrid.flatten(), az_s, el_s, zenith_luminance)
+    if meshgrid_flag:
+        rel_lum = compute_luminance(sky_type, az_meshgrid.flatten(), el_meshgrid.flatten(), az_s, el_s).reshape(az_meshgrid.shape)
+    else:
+        # rel_lum = compute_luminance(sky_type, az, el, az_s, el_s)
+        rel_lum = compute_luminance(sky_type, np.append(az, [360]), np.append(el, [90]), az_s, el_s)
+
     lum = rel_lum * zenith_luminance
 
+    print_n_max_values(rel_lum, 3)
+
     # Plot with pcolormesh
-    plot_luminance_map(az, el, lum.reshape(az_meshgrid.shape), az_s, el_s, luminance_or_radiance='luminance', polar=True)
+    plot_luminance_map(az, el, rel_lum,
+                       az_s, el_s,
+                       luminance_or_radiance='radiance',
+                       polar=True,
+                       colormap='jet')
 
     # Plot all 15 standard skies in a (5, 3) subplots figure (takes a couple of minutes)
     # plot_all_skies()
 
-    # Plot with PyVista
+if sky_mesh.lower() == 'reinhart':
+    # Plot with PyVista (only for Reinhart scheme)
     plot_3D(sky, rel_lum)
 
-- 
GitLab


From 00483a7141ad5dec4dc2bbb70d89773931eeeab9 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Mon, 7 Apr 2025 17:00:09 +0200
Subject: [PATCH 19/30] [feature] Compute Clear Sky Index and Cloudless Index
 and derive Sky type

Clear Sky Index (Kc) and Cloudless Index (Cle) are two meteorological indices from Igawa (2014). They are computed from weather data in the Light.data dataframe and the CIE sky type is derived from their values, based on Figure 11 of the paper (and I chose a standard sky type for each class : overcast, intermediate overcast, intermediate, intermediate clear and clear).
---
 MODULES/ENVIRONMENT/light.py | 94 +++++++++++++++++++++++++++++++++---
 1 file changed, 88 insertions(+), 6 deletions(-)

diff --git a/MODULES/ENVIRONMENT/light.py b/MODULES/ENVIRONMENT/light.py
index 0c6ee9e..c5e87e2 100644
--- a/MODULES/ENVIRONMENT/light.py
+++ b/MODULES/ENVIRONMENT/light.py
@@ -5,11 +5,14 @@
 #Authors : Roxane Bruhwyler (roxane.bruhwyler@uliege.be or roxane.bruhwyler@hotmail.com) and Nicolas De Cock (nicolas.decock1@gmail.com)
 #This file is part of the PASE software, and is distributed under the MIT license.
 
+from calendar import isleap
 import logging
 import numpy as np
 import pandas as pd
 import pyvista as pyV
 import pvlib.solarposition as pvlibSP
+from pvlib.atmosphere import get_relative_airmass
+from pvlib.irradiance import get_extra_radiation
 import matplotlib.pyplot as plt
 import os
 
@@ -84,7 +87,12 @@ class Sun_positions:
         return solar_vector
         
     def get_top_of_atm_radiation(self, index, n):
-        
+        """
+
+        :param index: datetime index
+        :param n: seems unused ?
+        :return: irradiance at top of atmosphere on a horizontal surface [W/m²]
+        """
         day_of_year = np.array(index.dayofyear)
         n_days_in_year = day_of_year[len(index)-1]
         solar_declination = self.get_solar_declination(day_of_year)
@@ -118,21 +126,37 @@ class Light:
     def __init__(self, WD, SP, ghi_multiplier=1):
         
         self.data = {}
+        sky_type_lut_path = os.path.join('INPUTS', 'Igawa-5_sky_types_lut.csv')
+        self.sky_type_lut = pd.read_csv(sky_type_lut_path, sep=';')
        
         for year in WD.keys():
             
             GHI = WD[year]['G(h)'].to_numpy() * ghi_multiplier
             
-            if int(year)%4 == 0:                
+            if isleap(int(year)):
                 rad_top_atm = SP.sp_leapY['Top_atm_radiation'].to_numpy()
+                apparent_sun_zenith = SP.sp_leapY['apparent_zenith'].to_numpy()
+                sun_elevation = SP.sp_leapY['elevation'].to_numpy()
             else:
                 rad_top_atm = SP.sp_nonleapY['Top_atm_radiation'].to_numpy()
-                
+                apparent_sun_zenith = SP.sp_nonleapY['apparent_zenith'].to_numpy()
+                sun_elevation = SP.sp_nonleapY['elevation'].to_numpy()
+
             kt = self.get_clearness_sky_index(rad_top_atm, GHI)    
             DHI = self.get_diffuse_horizontal_radiation(kt, GHI)
             BHI = self.get_beam_horizontal_radiation(GHI, DHI)
             Ai = self.get_anisotropy_index(rad_top_atm, BHI)
             f = self.get_modulating_factor(GHI, BHI)
+
+            # Extraterrestrial Normal Irradiance (used for the Kc and Cle below)
+            ENI = get_extra_radiation(WD[year]['DateTime'].dt.dayofyear.values)
+
+            # Meteorological indices Clear Sky Index (Kc) and Cloudless index
+            # (Cle) from Igawa (2014)
+            Kc = self.get_clear_sky_index(ENI, GHI, apparent_sun_zenith)
+            Cle = self.get_cloudless_index(DHI, GHI, sun_elevation)
+
+            cie_sky_type = self.get_sky_type(Kc, Cle)
             
             df = pd.DataFrame({'GHI': GHI.tolist(),
                                'BHI': BHI.tolist(),
@@ -140,7 +164,10 @@ class Light:
                                'rad_top_atm': rad_top_atm.tolist(),
                                'kt': kt.tolist(),
                                'Ai': Ai.tolist(),
-                               'f': f.tolist()}, index=WD[year].index)
+                               'f': f.tolist(),
+                               'Kc': Kc.tolist(),
+                               'Cle': Cle.tolist(),
+                               'CIE Sky Type': cie_sky_type}, index=WD[year].index)
             
             self.data[year] = df            
                     
@@ -195,8 +222,63 @@ class Light:
         f[GHI>0] = np.sqrt(BHI[GHI>0]/GHI[GHI>0])
         
         return f
-    
-        
+
+    def get_clear_sky_index(self, ENI, GHI, apparent_sun_zenith):
+
+        # Get relative airmass from pvlib method
+        # (takes apparent zenith angle in [deg] for Kasten and Young 1989 model)
+        m = get_relative_airmass(apparent_sun_zenith, model='kastenyoung1989')
+
+        GHI_non_zero = GHI != 0
+
+        Kc = np.empty(len(GHI))
+        Kc[:] = np.nan
+
+        Kc = np.divide(GHI, 0.84 * ENI / m * np.exp(-0.054 * m), out=Kc, where=GHI_non_zero)
+
+        return Kc
+
+    def get_cloudless_index(self, DHI, GHI, sun_elevation):
+
+        sun_elevation_rad = np.deg2rad(sun_elevation)
+        GHI_non_zero = GHI != 0
+
+        # Cloud ratio [-]
+        Ce = self.get_cloud_ratio(DHI, GHI, GHI_non_zero)
+
+        # Standard cloud ratio [-]
+        Ces = (0.08302
+               + 0.5358 * np.exp(-17.3 * sun_elevation_rad)
+               + 0.3818 * np.exp(-3.2899 * sun_elevation_rad))
+
+        Cle = np.zeros(len(GHI))
+        Cle = np.divide((1 - Ce), (1 - Ces), out=Cle, where=GHI_non_zero)
+
+        return Cle
+
+    def get_cloud_ratio(self, DHI, GHI, GHI_non_zero):
+
+        Ce = np.zeros((len(GHI),))
+
+        Ce = np.divide(DHI, GHI, out=Ce, where=GHI_non_zero)
+
+        return Ce
+
+    def get_sky_type(self, Kc_target, Cle_target):
+
+        mydf = pd.DataFrame({'Kc': Kc_target, 'Cle': Cle_target})
+        temp_list = []
+        for k, v in mydf.iterrows():
+            if np.isnan(v['Kc']):
+                temp_list.append(np.nan)
+            else:
+                i = ((self.sky_type_lut['Kc']-v['Kc']) *
+                     (self.sky_type_lut['Cle']-v['Cle'])).abs().idxmin()
+                temp_list.append(int(self.sky_type_lut['CIE Sky Type'].iloc[i]))
+
+        return temp_list
+
+
 class Sun_positions_sampled:
     
     def __init__(self, lat, long, precision_lvl, loc_name, freq_deter, TZ):
-- 
GitLab


From 0e1b75dd87e2b94cca9aa51a5d123ee5a49bd068 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Tue, 22 Apr 2025 16:27:00 +0200
Subject: [PATCH 20/30] [style] Some whitespaces

Add/remove some whitespaces
---
 MODULES/ENVIRONMENT/light.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/MODULES/ENVIRONMENT/light.py b/MODULES/ENVIRONMENT/light.py
index c5e87e2..29de7c6 100644
--- a/MODULES/ENVIRONMENT/light.py
+++ b/MODULES/ENVIRONMENT/light.py
@@ -169,7 +169,7 @@ class Light:
                                'Cle': Cle.tolist(),
                                'CIE Sky Type': cie_sky_type}, index=WD[year].index)
             
-            self.data[year] = df            
+            self.data[year] = df
                     
     def get_clearness_sky_index(self, rad_top_atm, GHI):
             
@@ -437,6 +437,7 @@ class Ray_casting_scene:
     def __init__(self,mesh,geometry):
         self.mesh = mesh
         self.sourcepoints = self.mesh.get_sourcepoints()
+
         self.geometry = geometry
         self.n_sourcepoints = self.sourcepoints.shape[0]
         self.sources_flag_dict = mesh.get_sources_flag_dict()
@@ -704,9 +705,9 @@ class Ray_casting_scene:
         dfShade = SP_sampled.tz_localize(None).reset_index()
         dfShade['doy'] = dfShade['index'].dt.dayofyear
         dfShade['RefDate'] = pd.to_datetime(dfShade['index'].dt.date)
- #       dfShade['hour'] = dfShade['index'].dt.hour
+        # dfShade['hour'] = dfShade['index'].dt.hour
         dfShade['second'] = pd.to_timedelta(dfShade['index'].dt.time.astype(str)).dt.total_seconds()
- #       dfShade['RefHour'] = dfShade['hour'] 
+        # dfShade['RefHour'] = dfShade['hour']
         dfShade['RefSecond'] = dfShade['second'] 
 
         for year in light_data: 
-- 
GitLab


From 5630ea52f6a3ceb963deca4e83a37f43fa144504 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Tue, 22 Apr 2025 16:28:28 +0200
Subject: [PATCH 21/30] [fix] Path of CIE standard skies csv file more robust

The definition of the file CIE_standard_skies.csv did not work with the module's implementation. pathlib Path object used to fix that.
---
 MODULES/ENVIRONMENT/sky_model.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/MODULES/ENVIRONMENT/sky_model.py b/MODULES/ENVIRONMENT/sky_model.py
index 05f601c..1cf6efc 100644
--- a/MODULES/ENVIRONMENT/sky_model.py
+++ b/MODULES/ENVIRONMENT/sky_model.py
@@ -11,6 +11,7 @@ import math
 import numpy as np
 import os
 import pandas as pd
+from pathlib import Path
 from matplotlib import pyplot as plt
 import matplotlib
 matplotlib.use('TkAgg')
@@ -419,7 +420,9 @@ class CIEStandardSky:
         """
 
         # Read the csv containing the descriptions of the 15 standard skies
-        standard_skies = pd.read_csv(os.path.join('..', '..',
+
+        pase_dir_path = Path(__file__).parents[2]
+        standard_skies = pd.read_csv(os.path.join(pase_dir_path,
                                                   'INPUTS',
                                                   'CIE_standard_skies.csv'))
 
-- 
GitLab


From 26e2913f038c1018d56980e06d9b92a9ee9f377c Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Tue, 22 Apr 2025 16:29:57 +0200
Subject: [PATCH 22/30] [feature] Add cos(zenith angle) column to the CIE sky
 dataframe

It is used many times in the irradiance integration (at each time step), which makes it more efficient to integrate it in the dataframe so we don't need to recompute that many times.
---
 MODULES/ENVIRONMENT/sky_model.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/MODULES/ENVIRONMENT/sky_model.py b/MODULES/ENVIRONMENT/sky_model.py
index 1cf6efc..cc49898 100644
--- a/MODULES/ENVIRONMENT/sky_model.py
+++ b/MODULES/ENVIRONMENT/sky_model.py
@@ -62,6 +62,9 @@ class ReinhartSky:
         # Compute x, y, z coordinates of unit vectors pointing towards the sky patches
         self.get_patches_xyz()
 
+        # Compute cos(zenith angle), used in irradiance computations
+        self.compute_cos_zenith()
+
     def get_original_tregenza(self):
         original_tregenza_segments = pd.DataFrame()
 
@@ -241,6 +244,9 @@ class ReinhartSky:
                                                  azimut=self.reinhart_patches['az'],
                                                  elev=self.reinhart_patches['el'])
 
+    def compute_cos_zenith(self):
+        self.reinhart_patches['cos(z)'] = np.cos(np.deg2rad(90-self.reinhart_patches['el']))
+
     def test_MFs(self):
         MFs = np.arange(1, 33)
         for MF in MFs:
-- 
GitLab


From 867ad5a2e19150a114ac76d5197a439eadb93929 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Tue, 22 Apr 2025 16:31:42 +0200
Subject: [PATCH 23/30] [fix] Precise column name "Normalized surf area" and
 fix type of variable sky_type

* The previous name "surf area" was misleading because it is actually a normalized area (i.e. the sum of all patches' areas is 1, not 2*pi as it should be for a half-sphere)
* the sky_type variable should be integer since it corresponds to an ID number
---
 MODULES/ENVIRONMENT/sky_model.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/MODULES/ENVIRONMENT/sky_model.py b/MODULES/ENVIRONMENT/sky_model.py
index cc49898..d0cd70d 100644
--- a/MODULES/ENVIRONMENT/sky_model.py
+++ b/MODULES/ENVIRONMENT/sky_model.py
@@ -221,9 +221,9 @@ class ReinhartSky:
 
     def compute_surface_areas(self):
         # Compute surf area of each patch
-        self.reinhart_patches['surf area'] = self.reinhart_patches['solid_angle_sr'] / (
+        self.reinhart_patches['Normalized surf area'] = self.reinhart_patches['solid_angle_sr'] / (
                 2 * np.pi)
-        sum_surf_areas = self.reinhart_patches['surf area'].sum()
+        sum_surf_areas = self.reinhart_patches['Normalized surf area'].sum()
 
         logger.debug(f'Sum of surface areas = {sum_surf_areas:.4g}')
 
@@ -312,7 +312,7 @@ class CIEStandardSky:
         self.az_s = sun_az  # [°] Sun azimuth (145° is south-east)
         self.el_s = sun_el  # [°]  Sun elevation
 
-        self.sky_type = sky_type
+        self.sky_type = int(sky_type)
 
         self.rel_radiance_distribution = self.compute_rel_radiance()
 
-- 
GitLab


From 817802b069e105a67aa6473b30a79e9664b65b72 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Tue, 22 Apr 2025 16:55:15 +0200
Subject: [PATCH 24/30] [feature] nclude sky radiance anisotropic distribution
 in diffuse irradiance computation

Before, the diffuse irradiance was computed based on :
- isotropic sky radiance
- computing which parts of the sky (sky patches) contribute to the diffuse irradiance by ray casting
- weighting the contribution of each sky patch with cos(zenith angle)

Now, the sky radiance is defined by the CIE sky type, based on Meteorological Indices "Clear Sky Index" and "Cloudless Index" (see Igawa, 2014 - link in issue #117).
- Each sky patch has its own radiance value, computed at each time step
- the weighting is now based on sky patch area and cos(zenith angle)

BREAKING CHANGE

Sun tracking is likely (untested) temporarily broken by this commit
---
 MODULES/ENVIRONMENT/light.py | 134 +++++++++++++++++------------------
 1 file changed, 67 insertions(+), 67 deletions(-)

diff --git a/MODULES/ENVIRONMENT/light.py b/MODULES/ENVIRONMENT/light.py
index 29de7c6..7f942cb 100644
--- a/MODULES/ENVIRONMENT/light.py
+++ b/MODULES/ENVIRONMENT/light.py
@@ -18,6 +18,7 @@ import os
 
 import MODULES.conversion_functions as cf
 from MODULES.ENVIRONMENT.sky_model import ReinhartSky, fibonacci_half_sphere
+from MODULES.ENVIRONMENT.sky_model import CIEStandardSky
 
 logger = logging.getLogger(__name__)
 
@@ -434,18 +435,20 @@ class Ray_casting_scene:
     The class is initiate with a mesh containing the source points and a geometry
     '''
     #The class light shade scene init with a geometry (pyvista.polydata) and a mesh instance
-    def __init__(self,mesh,geometry):
+    def __init__(self, mesh, geometry, discrete_sky, sun_pos_prec_lvl):
         self.mesh = mesh
         self.sourcepoints = self.mesh.get_sourcepoints()
 
         self.geometry = geometry
         self.n_sourcepoints = self.sourcepoints.shape[0]
         self.sources_flag_dict = mesh.get_sources_flag_dict()
+
+        self.discrete_sky = discrete_sky
+        self.sun_pos_prec_lvl = sun_pos_prec_lvl
         
-        
-    def get_light_maps(self, sun_P, scheme='Reinhart', MF=1, n_small_suns=180):
+    def get_light_maps(self, sun_P):
     
-        if type(self.geometry) == list:
+        if type(self.geometry) == list:  # sun tracking
             sv = np.zeros((1,3))
             self.dir_map = np.zeros((len(self.sourcepoints),len(sun_P[:,0])))
             self.diff_map = np.zeros((len(self.sourcepoints),len(sun_P[:,0])))
@@ -453,19 +456,13 @@ class Ray_casting_scene:
             for time in range(len(sun_P[:,0])):
                 print(time)
                 geometry = self.geometry[time]
-                difff_map = self.diffuse_map(geometry,
-                                             scheme=scheme,
-                                             MF=MF,
-                                             n_small_suns=n_small_suns)
+                difff_map = self.diffuse_map(geometry)
                 sv[0,:] = sun_P[time,:]
                 dirrr_map = self.direct_map(sv, geometry)
                 self.dir_map[:,time] = dirrr_map 
                 self.diff_map[:,time] = difff_map
-        else:
-            self.diff_map = self.diffuse_map(self.geometry,
-                                             scheme=scheme,
-                                             MF=MF,
-                                             n_small_suns=n_small_suns)
+        else:  # no sun tracking
+            self.diff_map = self.diffuse_map(self.geometry)
             self.dir_map = self.direct_map(sun_P, self.geometry)
 
     def self_intercept(self,SourcePoints,intercept_points,id_rays_stopped,tol = 0.01):
@@ -488,45 +485,33 @@ class Ray_casting_scene:
         return np.unique(id_rays_stopped[delta>tol])
         
  
-    def diffuse_map(self, geometry, scheme='Reinhart', MF=1, n_small_suns=180):
+    def diffuse_map(self, geometry):
         """
-        Public method, compute the diffuse light at the point sources defined in the input mesh using
-         the approximation of a isotropic half sphere sky. 
-         The method uses the mesh and the geometry set at the initialization of the instance
-        
-        Parameters:
-            n_small_suns (int): number of sources consider in the sky for the diffuse light computation
-            higher number will provide a better accuracy but heavier computation
+        Public method, compute the diffuse light at the point sources defined
+        in the input mesh (see this class' constructor) under a
+        discrete_sky sky model.
+
+        :param geometry: geometry set at the initialization of the instance
+        :return: diffuse_mask: diffuse map binary mask for each source point and discrete_sky element. Shape: (n_sky_patches, n_source_points)
 
-        Returns:
-           Diffu (np.array 1 x n):  Providing a vector with the fraction ([0-1]) of diffuse light 
-                                   for each of the "n" source points defined in the mesh
         """
-      
-    
+
         #Get direction of ray to reach the small suns and compute the sky view of each point
-        if scheme.lower() == 'reinhart':
-            sky = ReinhartSky(MF=MF)
-            pTarget = np.column_stack([sky.reinhart_patches.x,
-                                       sky.reinhart_patches.y,
-                                       sky.reinhart_patches.z])
-            n_small_suns = len(sky.reinhart_patches)
-        elif scheme.lower() == 'fibonacci':
-            pTarget = fibonacci_half_sphere(n_small_suns)
-        else:
-            raise NotImplementedError('Unrecognized sky discretization scheme')
-        
+        pTarget = np.column_stack([self.discrete_sky.x,
+                                   self.discrete_sky.y,
+                                   self.discrete_sky.z])
+        n_sky_elements = len(self.discrete_sky)
         
-        #Creation of the source points array (Nx3) with N = len(Source) * len(n_small_suns)
+        #Creation of the source points array (Nx3) with N = len(Source) * len(n_sky_elements)
         SourcePoints = np.repeat(np.column_stack((
                                                   self.sourcepoints[:,0],
                                                   self.sourcepoints[:,1],
                                                   self.sourcepoints[:,2]
                                                  )),
-                                      n_small_suns,
+                                      n_sky_elements,
                                       axis=0)
         
-        #Creation of the target points array (Nx3) with N = len(Source) * len(n_small_suns)
+        #Creation of the target points array (Nx3) with N = len(Source) * len(n_sky_elements)
         TargetPoints = np.tile(pTarget,[self.n_sourcepoints,1])
         
         #Computation of the ray interception of the N rays
@@ -544,7 +529,7 @@ class Ray_casting_scene:
         SourceID = np.repeat(np.linspace(0,
                                          self.n_sourcepoints-1,
                                          self.n_sourcepoints),
-                             n_small_suns,
+                             n_sky_elements,
                              axis=0)
         
         #Touched provide a list with the sourceID of the intercept ray
@@ -555,24 +540,17 @@ class Ray_casting_scene:
 
         # Compute the cos(zenith angle) of all the small suns for the normalization
         _, _zenith_angle_all = cf.get_zenith_angle_from_cart(TargetPoints)
-        _cos_zenith_angle_all = np.cos(_zenith_angle_all).reshape(self.n_sourcepoints, n_small_suns)
+        _cos_zenith_angle_all = np.cos(_zenith_angle_all).reshape(self.n_sourcepoints, n_sky_elements)
 
         # Compute the cos(zenith angle) of the small suns that do NOT contribute to the diffuse map
         # (i.e. rays that were intercepted)
 
-        _mask = np.ones(_cos_zenith_angle_all.size, bool)
-        _mask[id_rays_stopped_filtred] = 0
-        _cos_zenith_angle_blocked = _cos_zenith_angle_all.copy()
-        _mask = _mask.reshape(self.n_sourcepoints, n_small_suns)
-        _cos_zenith_angle_blocked[_mask] = 0
-
-        #Creation of the empty matrix of sky view
-        Diffu = np.ones(self.n_sourcepoints, dtype=np.float16)
+        diffuse_mask = np.ones(_cos_zenith_angle_all.size, bool)
+        diffuse_mask[id_rays_stopped_filtred] = 0
 
-        #Computation of the sky view by removing the fraction of intercepted ray at each location
-        Diffu[unique.astype("int")] = 1 - np.sum(_cos_zenith_angle_blocked[unique.astype("int"), :], axis=1)/np.sum(_cos_zenith_angle_all[unique.astype("int"), :], axis=1)
+        diffuse_mask = diffuse_mask.reshape(self.n_sourcepoints, n_sky_elements)
 
-        return Diffu
+        return diffuse_mask
 
     def get_direct_map_by_flag(self,Flags):
         """
@@ -721,10 +699,8 @@ class Ray_casting_scene:
                 n = 6
                 
   
-            #dfWeather = light_data[year].tz_localize('Etc/GMT+0').reset_index()
             dfWeather = light_data[year].tz_localize(None).reset_index()
             dfWeather['doy'] = dfWeather['index'].dt.dayofyear
-     #       dfWeather['hour'] = dfWeather['index'].dt.hour
             dfWeather['second'] = pd.to_timedelta(dfWeather['index'].dt.time.astype(str)).dt.total_seconds()
 
             # Jointure entre les données météos et les données d'ombrage en vue de déterminée la date/index liée aux données d'ombrage
@@ -733,34 +709,58 @@ class Ray_casting_scene:
             irradianceMap_direct = {}
             irradianceMap_diffus = {}
             
-            # Boucle sur les n jours de l'annee afin de calculer l'irradiation journaliere
+            # Loop over days of year (doy) to compute daily irradiance.
             doy = dfWeatherMerged['index'].dt.dayofyear.unique()
-            doy.sort()
+            doy.sort()  # doy = [1, 2, 3, ..., 365]
             for day in doy:
                 
-                #Creation de sous dataframe comprenant les donnees du jour numero "doy"
+                #Creation de sous dataframe comprenant les donnees du jour numero "day"
                 sublight_df = dfWeatherMerged.loc[dfWeatherMerged['index'].dt.dayofyear==day,:]
                 df_subShade = dfShade.loc[sublight_df['RefDate'].unique()[0]==dfShade['RefDate']]
                 df_subShade_merged = df_subShade.join(sublight_df[['second','BHI','GHI','DHI']].set_index('second'),on='second',how='left',rsuffix='_',lsuffix="__")
                 
                 #df_subShade_merged = df_subShade.join(sublight_df[['hour','BHI','GHI','DHI']].set_index('hour'),on='hour',how='left',rsuffix='_',lsuffix="__")
                 #Jointure sur l'instant de la journee la plus proche sur base des secondes écoulées depuis le debut de la journee
-                df_subShade_merged = pd.merge_asof(df_subShade,sublight_df[['second','BHI','GHI','DHI','doy']].set_index('second'),on=['second'],direction='nearest',suffixes=('_x','_y')).sort_values('second')
+                df_subShade_merged = pd.merge_asof(df_subShade,sublight_df[['second','BHI','GHI','DHI','doy', 'CIE Sky Type']].set_index('second'),on=['second'],direction='nearest',suffixes=('_x','_y')).sort_values('second')
 
                 #Calcul de l'irradiation directe et diffuse dans ce jour et injection de la donnee dans le dictionnaire lie
                 irradianceMap_direct[day] = np.sum(self.dir_map[:,list(df_subShade.index)] * df_subShade_merged['BHI'].to_numpy(),axis=1)*10**-6*60*60/n
-                if type(self.geometry) != list:
-                    irradianceMap_diffus[day] = np.sum(df_subShade_merged['DHI'].to_numpy())*self.diff_map*10**-6*60*60/n
-                else:
-                    print('day')
+                if type(self.geometry) != list:  # no sun tracking
+                    irradianceMap_diffus[day] = self.compute_daily_diff_irradiation(df_subShade_merged.dropna())
+                else:  # sun tracking -> in this case, the mask changes at each time step
+                    print(f'{day=}')
                     irradianceMap_diffus[day] = np.sum(df_subShade_merged['DHI'].to_numpy()*self.diff_map[:,list(df_subShade.index)], axis=1)*10**-6*60*60/n
-            
+
             #Conversion des dictionnaires en matrice numpy et ajout dans l attribut ad-hoc
-            self.daily_irr_spat[year] = pd.DataFrame.from_dict(irradianceMap_diffus).to_numpy()    + pd.DataFrame.from_dict(irradianceMap_direct).to_numpy() 
-            self.daily_dir_irr_spat[year] = pd.DataFrame.from_dict(irradianceMap_direct).to_numpy()    
+            self.daily_irr_spat[year] = (pd.DataFrame.from_dict(irradianceMap_diffus).to_numpy()
+                                         + pd.DataFrame.from_dict(irradianceMap_direct).to_numpy())
+            self.daily_dir_irr_spat[year] = pd.DataFrame.from_dict(irradianceMap_direct).to_numpy()
             self.daily_diff_irr_spat[year] = pd.DataFrame.from_dict(irradianceMap_diffus).to_numpy()     
 
-    
+    def compute_daily_diff_irradiation(self, df):
+        # loop over each timestep in the day's dataframes
+        for i in range(len(df)):
+            unweighted_radiance = CIEStandardSky(self.discrete_sky,
+                                                 df['azimuth'].iloc[i],
+                                                 df['elevation'].iloc[i],
+                                                 df['CIE Sky Type'].iloc[i]).rel_radiance_distribution
+
+            # sky patch area and cos(z) weighting
+            weighted_radiance = unweighted_radiance * self.discrete_sky['Normalized surf area'] * self.discrete_sky['cos(z)']
+
+            # normalize to relative radiance contributions
+            rel_radiance_contrib = weighted_radiance/np.sum(weighted_radiance)
+
+            # multiply by the mask
+            if type(self.geometry) == list:  # sun tracking
+                raise NotImplementedError('Sun tracking is temporarily broken')
+                # shaded_radiance_contrib = rel_radiance_contrib * self.diff_map[:,time]
+            else:
+                shaded_radiance_contrib = rel_radiance_contrib.values * self.diff_map
+
+            diff_irradiance_map = (shaded_radiance_contrib * df['DHI'].iloc[i]).sum(axis=1)
+
+            return diff_irradiance_map
 
 def show_light_map(light_matrix, msh_grid, PV_central, lim_min, lim_max, lgd_title):
     
-- 
GitLab


From 873c6ad428a5affc0ca79625ff07098838dd4c8d Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Tue, 22 Apr 2025 16:57:16 +0200
Subject: [PATCH 25/30] [feature] Instantiate discrete sky in the example
 script and adjust inputs of Ray_casting_scene instantiation

n/a
---
 example.py | 16 ++++++++++------
 1 file changed, 10 insertions(+), 6 deletions(-)

diff --git a/example.py b/example.py
index 5c384b1..ebd86c6 100644
--- a/example.py
+++ b/example.py
@@ -17,6 +17,7 @@ from MODULES.PHOTOVOLTAICS.configuration import PV_Configuration_3D
 from MODULES.ENVIRONMENT.light import Sun_positions_sampled, Sun_positions, Light
 from MODULES.ENVIRONMENT.light import show_light_map2, Ray_casting_scene
 from MODULES.ENVIRONMENT.mesh import Mesh
+from MODULES.ENVIRONMENT.sky_model import ReinhartSky
 from MODULES.PHOTOVOLTAICS.production import PV_Production
 from MODULES.CROPS.run_crop_simulations import run_crop_simu
 
@@ -62,7 +63,7 @@ Sun_positions_complete = Sun_positions(Loc_1['Latitude'],
 # Instantiation of the 3D PV central
 PV_1_3Dconfig = PV_Configuration_3D(PV_params_dict,
                                     Sun_positions_samp.solar_vector,
-                                    visualization=True)  # !!!! Problem with rotation angle that are negative
+                                    visualization=False)  # !!!! Problem with rotation angle that are negative
 
 # Initiation of the object containing points of interest to compute light
 M = Mesh()
@@ -76,15 +77,18 @@ M.add_plane_ground_regular_meshes(Loc_1['Xmin_InterestZone'],
                                   Loc_1['dY_InterestZone'],
                                   flag="crop")
 
+# Discrete sky model
+discrete_sky = ReinhartSky(MF=Loc_1['MF']).reinhart_patches
+
 # Computation of sun and light data
 Light_instance = Light(WD.nyears_data, Sun_positions_complete)
 
 # Iniation and run of light ray casting model (direct and diffuse) with points of interest and scene
-L = Ray_casting_scene(mesh=M, geometry=PV_1_3Dconfig.PV_central_PD)
-L.get_light_maps(Sun_positions_samp.solar_vector,
-                 scheme=Loc_1['SkyDiscretizationScheme'],
-                 MF=Loc_1['MF'],
-                 n_small_suns=Loc_1['FibonacciSamples'])
+L = Ray_casting_scene(mesh=M,
+                      geometry=PV_1_3Dconfig.PV_central_PD,
+                      discrete_sky=discrete_sky)
+
+L.get_light_maps(Sun_positions_samp.solar_vector)
 
 # Integration of irradiation along days
 L.get_daily_irradiation_map(Sun_positions_samp.SP,
-- 
GitLab


From 62ba07aedccc4c3b53dbbe11fe7a170ebc51e409 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Tue, 22 Apr 2025 16:57:53 +0200
Subject: [PATCH 26/30] [fix] Remove unused input of Ray_casting_scene

sun position precision level is not used
---
 MODULES/ENVIRONMENT/light.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/MODULES/ENVIRONMENT/light.py b/MODULES/ENVIRONMENT/light.py
index 7f942cb..207fdde 100644
--- a/MODULES/ENVIRONMENT/light.py
+++ b/MODULES/ENVIRONMENT/light.py
@@ -435,7 +435,7 @@ class Ray_casting_scene:
     The class is initiate with a mesh containing the source points and a geometry
     '''
     #The class light shade scene init with a geometry (pyvista.polydata) and a mesh instance
-    def __init__(self, mesh, geometry, discrete_sky, sun_pos_prec_lvl):
+    def __init__(self, mesh, geometry, discrete_sky):
         self.mesh = mesh
         self.sourcepoints = self.mesh.get_sourcepoints()
 
@@ -444,8 +444,7 @@ class Ray_casting_scene:
         self.sources_flag_dict = mesh.get_sources_flag_dict()
 
         self.discrete_sky = discrete_sky
-        self.sun_pos_prec_lvl = sun_pos_prec_lvl
-        
+
     def get_light_maps(self, sun_P):
     
         if type(self.geometry) == list:  # sun tracking
-- 
GitLab


From 07fb0d97a6fd3055697a8351ed4ccccd64890a81 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Tue, 22 Apr 2025 16:59:02 +0200
Subject: [PATCH 27/30] [feature] Add field DiffuseSkyType to config file

see Definition
---
 INPUTS/SCENARIOS/Example1_loc.yaml | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/INPUTS/SCENARIOS/Example1_loc.yaml b/INPUTS/SCENARIOS/Example1_loc.yaml
index 78d08b5..dee079f 100644
--- a/INPUTS/SCENARIOS/Example1_loc.yaml
+++ b/INPUTS/SCENARIOS/Example1_loc.yaml
@@ -71,6 +71,14 @@ FibonacciSamples:
    Type: integer
    Limit: [1,100000]
    Definition: Samples in the Fibonacci sky subdivision scheme.
+DiffuseSkyType:
+   Value: From weather data
+   Type: string
+   Possibilities: [From weather data, Uniform]
+   Definition: Determine how to compute sky diffuse radiance distribution. 
+               "From weather data" uses a simplified implementation of Igawa (2014) 
+               to derive the sky type from weather data.
+               "Uniform" assumes uniform radiance distribution (faster but less accurate).
 Xmin_InterestZone:
    Value: -2
    Type: float
-- 
GitLab


From 0182820b0d9937817b14148d15055c88794c4bb0 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Thu, 24 Apr 2025 09:14:05 +0200
Subject: [PATCH 28/30] [feature] Try a derivation of absolute radiance
 distribution

Attempt in the script part of sky_model, no impact on the PASE codebase.
---
 MODULES/ENVIRONMENT/sky_model.py | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/MODULES/ENVIRONMENT/sky_model.py b/MODULES/ENVIRONMENT/sky_model.py
index d0cd70d..ac9d75f 100644
--- a/MODULES/ENVIRONMENT/sky_model.py
+++ b/MODULES/ENVIRONMENT/sky_model.py
@@ -498,4 +498,12 @@ if __name__ == '__main__':
             sky.scatter_3D(text_id=text_id)
 
     sky.reinhart_patches['Relative radiance distribution'] = CIEStandardSky(sky.reinhart_patches,
-                               38.02, 147.67, 15).rel_radiance_distribution
+                               sun_az=147.67, sun_el=38.02, sky_type=15).rel_radiance_distribution
+
+    DHI = 100  # [W/m²] arbitrary value for example's sake
+    integral_rel_sky_radiance = (sky.reinhart_patches['Relative radiance distribution']
+                                 * sky.reinhart_patches['solid_angle_sr']).sum()
+    zenith_radiance = DHI/integral_rel_sky_radiance
+    sky.reinhart_patches['Absolute radiance distribution'] = (sky.reinhart_patches['Relative radiance distribution']
+                                                              * zenith_radiance)
+
-- 
GitLab


From 4ac94c6a1d23aa9abfb07c371f2568e97303d893 Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Thu, 24 Apr 2025 14:10:12 +0200
Subject: [PATCH 29/30] [feature] Normalize the displayed diffuse map

Normalize the diffuse map to reflect the proportion of sky patches seen by each ground control point (so that the scalar values lie between 0 and 1 for respectively no obstacles at all and totally obstructed sky).
---
 example.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/example.py b/example.py
index ebd86c6..50796fb 100644
--- a/example.py
+++ b/example.py
@@ -104,8 +104,11 @@ show_light_map2(L.sourcepoints[:,:-1],
                 PV_1_3Dconfig.PV_central_PD,
                 "Direct map [-]")
 
+# Display the Sky visibility map (= fraction of sky patches viewed by each control point)
+# It is purely derived from the geometry of the central
+# (i.e. computed once, or for each tracking position if tracking is activated)
 show_light_map2(L.sourcepoints[:,:-1], 
-                np.array(L.diff_map,dtype=np.float32), 
+                np.array(L.diff_map.sum(axis=1)/len(discrete_sky),dtype=np.float32),
                 PV_1_3Dconfig.PV_central_PD,
                 "Sky visibility map [-]")
 
-- 
GitLab


From 1080f9acdefcc548dfa6509c83782efae8b6e75e Mon Sep 17 00:00:00 2001
From: Arnaud Bouvry <abouvry@uliege.be>
Date: Thu, 24 Apr 2025 15:36:49 +0200
Subject: [PATCH 30/30] [feature] Add csv files used to compute sky radiance
 distribution

They were ignored because of the .gitignore file
---
 INPUTS/CIE_standard_skies.csv    |  16 +
 INPUTS/Igawa-5_sky_types_lut.csv | 599 +++++++++++++++++++++++++++++++
 2 files changed, 615 insertions(+)
 create mode 100644 INPUTS/CIE_standard_skies.csv
 create mode 100644 INPUTS/Igawa-5_sky_types_lut.csv

diff --git a/INPUTS/CIE_standard_skies.csv b/INPUTS/CIE_standard_skies.csv
new file mode 100644
index 0000000..463d374
--- /dev/null
+++ b/INPUTS/CIE_standard_skies.csv
@@ -0,0 +1,16 @@
+Type,Gradation,Indikatrix,a,b,c,d,e,Description of luminance distribution
+1,1,1,4.0 ,-0.70 ,0,-1.0 ,0.00 ,"CIE Standard Overcast Sky, alternative form. Steep luminance gradation towards zenith, azimuthal uniformity"
+2,1,2,4.0 ,-0.70 ,2,-1.5 ,0.15 ,"Overcast, with steep luminance gradation and slight brightening towards the sun"
+3,2,1,1.1 ,-0.8 ,0,-1.0 ,0.00 ,"Overcast, moderately graded with azimuthal uniformity"
+4,2,2,1.1 ,-0.8 ,2,-1.5 ,0.15 ,"Overcast, moderately graded and slight brightening towards the sun"
+5,3,1,0.0 ,-1.0 ,0,-1.0 ,0.00 ,"Sky of uniform luminance"
+6 ,3,2,0.0 ,-1.0 ,2,-1.5 ,0.15 ,"Partly cloudy sky, no gradation towards zenith, slight brightening towards the sun"
+7,3,3,0.0 ,-1.0 ,5,-2.5 ,0.30 ,"Partly cloudy sky, no gradation towards zenith, brighter circumsolar region"
+8,3,4,0.0 ,-1.0 ,10,-3.0 ,0.45 ,"Partly cloudy sky, no gradation towards zenith, distinct solar corona"
+9,4,2,-1.0 ,-0.55 ,2,-1.5 ,0.15 ,"Partly cloudy, with the obscured sun"
+10,4,3,-1.0 ,-0.55 ,5,-2.5 ,0.30 ,"Partly cloudy, with brighter circumsolar region"
+11,4,4,-1.0 ,-0.55 ,10,-3.0 ,0.45 ,"White-blue sky with distinct solar corona"
+12,5,4,-1.0 ,-0.32 ,10,-3.0 ,0.45 ,"CIE Standard Clear Sky, low illuminance turbidity"
+13 ,5,5,-1.0 ,-0.32 ,16,-3.0 ,0.30 ,"CIE Standard Clear Sky, polluted atmosphere"
+14 ,6,5,-1.0 ,-0.15 ,16,-3.0 ,0.30 ,"Cloudless turbid sky with broad solar corona"
+15,6,6,-1.0 ,-0.15 ,24,-2.8 ,0.15 ,"White-blue turbid sky with broad solar corona"
diff --git a/INPUTS/Igawa-5_sky_types_lut.csv b/INPUTS/Igawa-5_sky_types_lut.csv
new file mode 100644
index 0000000..7a0f0c9
--- /dev/null
+++ b/INPUTS/Igawa-5_sky_types_lut.csv
@@ -0,0 +1,599 @@
+;Kc;Cle;Sky category;CIE Sky Type
+0;0.01;0.0;O;1
+1;0.05;0.0;O;1
+2;0.1;0.0;O;1
+3;0.15;0.0;O;1
+4;0.2;0.0;O;1
+5;0.25;0.0;O;1
+6;0.3;0.0;O;1
+7;0.35;0.0;O;1
+8;0.4;0.0;O;1
+9;0.45;0.0;IO;4
+10;0.5;0.0;IO;4
+11;0.55;0.0;IO;4
+12;0.6;0.0;IO;4
+13;0.65;0.0;IO;4
+14;0.7;0.0;IO;4
+15;0.75;0.0;;5
+16;0.8;0.0;;5
+17;0.85;0.0;;5
+18;0.9;0.0;;5
+19;0.95;0.0;;5
+20;1.0;0.0;;5
+21;1.05;0.0;;5
+22;1.1;0.0;;5
+23;1.15;0.0;;5
+24;1.2;0.0;;5
+25;1.25;0.0;;5
+26;0.01;0.05;O;1
+27;0.05;0.05;O;1
+28;0.1;0.05;O;1
+29;0.15;0.05;O;1
+30;0.2;0.05;IO;4
+31;0.25;0.05;IO;4
+32;0.3;0.05;IO;4
+33;0.35;0.05;IO;4
+34;0.4;0.05;IO;4
+35;0.45;0.05;IO;4
+36;0.5;0.05;IO;4
+37;0.55;0.05;I;7
+38;0.6;0.05;I;7
+39;0.65;0.05;I;7
+40;0.7;0.05;I;7
+41;0.75;0.05;I;7
+42;0.8;0.05;I;7
+43;0.85;0.05;I;7
+44;0.9;0.05;;5
+45;0.95;0.05;;5
+46;1.0;0.05;;5
+47;1.05;0.05;;5
+48;1.1;0.05;;5
+49;1.15;0.05;;5
+50;1.2;0.05;;5
+51;1.25;0.05;;5
+52;0.01;0.1;O;1
+53;0.05;0.1;O;1
+54;0.1;0.1;IO;4
+55;0.15;0.1;IO;4
+56;0.2;0.1;IO;4
+57;0.25;0.1;IO;4
+58;0.3;0.1;IO;4
+59;0.35;0.1;IO;4
+60;0.4;0.1;IO;4
+61;0.45;0.1;I;7
+62;0.5;0.1;I;7
+63;0.55;0.1;I;7
+64;0.6;0.1;I;7
+65;0.65;0.1;I;7
+66;0.7;0.1;I;7
+67;0.75;0.1;I;7
+68;0.8;0.1;I;7
+69;0.85;0.1;I;7
+70;0.9;0.1;I;7
+71;0.95;0.1;;5
+72;1.0;0.1;;5
+73;1.05;0.1;;5
+74;1.1;0.1;;5
+75;1.15;0.1;;5
+76;1.2;0.1;;5
+77;1.25;0.1;;5
+78;0.01;0.15;O;1
+79;0.05;0.15;IO;4
+80;0.1;0.15;IO;4
+81;0.15;0.15;IO;4
+82;0.2;0.15;IO;4
+83;0.25;0.15;IO;4
+84;0.3;0.15;IO;4
+85;0.35;0.15;I;7
+86;0.4;0.15;I;7
+87;0.45;0.15;I;7
+88;0.5;0.15;I;7
+89;0.55;0.15;I;7
+90;0.6;0.15;I;7
+91;0.65;0.15;I;7
+92;0.7;0.15;I;7
+93;0.75;0.15;I;7
+94;0.8;0.15;I;7
+95;0.85;0.15;I;7
+96;0.9;0.15;I;7
+97;0.95;0.15;;5
+98;1.0;0.15;;5
+99;1.05;0.15;;5
+100;1.1;0.15;;5
+101;1.15;0.15;;5
+102;1.2;0.15;;5
+103;1.25;0.15;;5
+104;0.01;0.2;IO;4
+105;0.05;0.2;IO;4
+106;0.1;0.2;IO;4
+107;0.15;0.2;IO;4
+108;0.2;0.2;IO;4
+109;0.25;0.2;IO;4
+110;0.3;0.2;I;7
+111;0.35;0.2;I;7
+112;0.4;0.2;I;7
+113;0.45;0.2;I;7
+114;0.5;0.2;I;7
+115;0.55;0.2;I;7
+116;0.6;0.2;I;7
+117;0.65;0.2;I;7
+118;0.7;0.2;I;7
+119;0.75;0.2;I;7
+120;0.8;0.2;I;7
+121;0.85;0.2;I;7
+122;0.9;0.2;I;7
+123;0.95;0.2;I;7
+124;1.0;0.2;;5
+125;1.05;0.2;;5
+126;1.1;0.2;;5
+127;1.15;0.2;;5
+128;1.2;0.2;;5
+129;1.25;0.2;;5
+130;0.01;0.25;;5
+131;0.05;0.25;;5
+132;0.1;0.25;;5
+133;0.15;0.25;;5
+134;0.2;0.25;;5
+135;0.25;0.25;IO;4
+136;0.3;0.25;I;7
+137;0.35;0.25;I;7
+138;0.4;0.25;I;7
+139;0.45;0.25;I;7
+140;0.5;0.25;I;7
+141;0.55;0.25;I;7
+142;0.6;0.25;I;7
+143;0.65;0.25;I;7
+144;0.7;0.25;I;7
+145;0.75;0.25;I;7
+146;0.8;0.25;I;7
+147;0.85;0.25;I;7
+148;0.9;0.25;I;7
+149;0.95;0.25;I;7
+150;1.0;0.25;I;7
+151;1.05;0.25;I;7
+152;1.1;0.25;I;7
+153;1.15;0.25;;5
+154;1.2;0.25;;5
+155;1.25;0.25;;5
+156;0.01;0.3;;5
+157;0.05;0.3;;5
+158;0.1;0.3;;5
+159;0.15;0.3;;5
+160;0.2;0.3;;5
+161;0.25;0.3;I;7
+162;0.3;0.3;I;7
+163;0.35;0.3;I;7
+164;0.4;0.3;I;7
+165;0.45;0.3;I;7
+166;0.5;0.3;I;7
+167;0.55;0.3;I;7
+168;0.6;0.3;I;7
+169;0.65;0.3;I;7
+170;0.7;0.3;I;7
+171;0.75;0.3;I;7
+172;0.8;0.3;I;7
+173;0.85;0.3;I;7
+174;0.9;0.3;I;7
+175;0.95;0.3;I;7
+176;1.0;0.3;I;7
+177;1.05;0.3;I;7
+178;1.1;0.3;I;7
+179;1.15;0.3;I;7
+180;1.2;0.3;;5
+181;1.25;0.3;;5
+182;0.01;0.35;;5
+183;0.05;0.35;;5
+184;0.1;0.35;;5
+185;0.15;0.35;;5
+186;0.2;0.35;;5
+187;0.25;0.35;;5
+188;0.3;0.35;I;7
+189;0.35;0.35;I;7
+190;0.4;0.35;I;7
+191;0.45;0.35;I;7
+192;0.5;0.35;I;7
+193;0.55;0.35;I;7
+194;0.6;0.35;I;7
+195;0.65;0.35;I;7
+196;0.7;0.35;I;7
+197;0.75;0.35;I;7
+198;0.8;0.35;I;7
+199;0.85;0.35;I;7
+200;0.9;0.35;I;7
+201;0.95;0.35;I;7
+202;1.0;0.35;I;7
+203;1.05;0.35;I;7
+204;1.1;0.35;I;7
+205;1.15;0.35;I;7
+206;1.2;0.35;;5
+207;1.25;0.35;;5
+208;0.01;0.4;;5
+209;0.05;0.4;;5
+210;0.1;0.4;;5
+211;0.15;0.4;;5
+212;0.2;0.4;;5
+213;0.25;0.4;;5
+214;0.3;0.4;I;7
+215;0.35;0.4;I;7
+216;0.4;0.4;I;7
+217;0.45;0.4;I;7
+218;0.5;0.4;I;7
+219;0.55;0.4;I;7
+220;0.6;0.4;I;7
+221;0.65;0.4;I;7
+222;0.7;0.4;I;7
+223;0.75;0.4;I;7
+224;0.8;0.4;I;7
+225;0.85;0.4;I;7
+226;0.9;0.4;I;7
+227;0.95;0.4;I;7
+228;1.0;0.4;I;7
+229;1.05;0.4;I;7
+230;1.1;0.4;I;7
+231;1.15;0.4;I;7
+232;1.2;0.4;I;7
+233;1.25;0.4;;5
+234;0.01;0.45;;5
+235;0.05;0.45;;5
+236;0.1;0.45;;5
+237;0.15;0.45;;5
+238;0.2;0.45;;5
+239;0.25;0.45;;5
+240;0.3;0.45;;5
+241;0.35;0.45;I;7
+242;0.4;0.45;I;7
+243;0.45;0.45;I;7
+244;0.5;0.45;I;7
+245;0.55;0.45;I;7
+246;0.6;0.45;I;7
+247;0.65;0.45;I;7
+248;0.7;0.45;I;7
+249;0.75;0.45;I;7
+250;0.8;0.45;I;7
+251;0.85;0.45;I;7
+252;0.9;0.45;I;7
+253;0.95;0.45;I;7
+254;1.0;0.45;I;7
+255;1.05;0.45;I;7
+256;1.1;0.45;I;7
+257;1.15;0.45;I;7
+258;1.2;0.45;I;7
+259;1.25;0.45;I;7
+260;0.01;0.5;;5
+261;0.05;0.5;;5
+262;0.1;0.5;;5
+263;0.15;0.5;;5
+264;0.2;0.5;;5
+265;0.25;0.5;;5
+266;0.3;0.5;;5
+267;0.35;0.5;;5
+268;0.4;0.5;;5
+269;0.45;0.5;I;7
+270;0.5;0.5;I;7
+271;0.55;0.5;I;7
+272;0.6;0.5;I;7
+273;0.65;0.5;I;7
+274;0.7;0.5;I;7
+275;0.75;0.5;I;7
+276;0.8;0.5;I;7
+277;0.85;0.5;I;7
+278;0.9;0.5;I;7
+279;0.95;0.5;IC;11
+280;1.0;0.5;IC;11
+281;1.05;0.5;IC;11
+282;1.1;0.5;I;7
+283;1.15;0.5;I;7
+284;1.2;0.5;I;7
+285;1.25;0.5;I;7
+286;0.01;0.55;;5
+287;0.05;0.55;;5
+288;0.1;0.55;;5
+289;0.15;0.55;;5
+290;0.2;0.55;;5
+291;0.25;0.55;;5
+292;0.3;0.55;;5
+293;0.35;0.55;;5
+294;0.4;0.55;;5
+295;0.45;0.55;;5
+296;0.5;0.55;I;7
+297;0.55;0.55;I;7
+298;0.6;0.55;I;7
+299;0.65;0.55;I;7
+300;0.7;0.55;I;7
+301;0.75;0.55;I;7
+302;0.8;0.55;I;7
+303;0.85;0.55;IC;11
+304;0.9;0.55;IC;11
+305;0.95;0.55;IC;11
+306;1.0;0.55;IC;11
+307;1.05;0.55;IC;11
+308;1.1;0.55;IC;11
+309;1.15;0.55;IC;11
+310;1.2;0.55;I;7
+311;1.25;0.55;I;7
+312;0.01;0.6;;5
+313;0.05;0.6;;5
+314;0.1;0.6;;5
+315;0.15;0.6;;5
+316;0.2;0.6;;5
+317;0.25;0.6;;5
+318;0.3;0.6;;5
+319;0.35;0.6;;5
+320;0.4;0.6;;5
+321;0.45;0.6;;5
+322;0.5;0.6;I;7
+323;0.55;0.6;I;7
+324;0.6;0.6;I;7
+325;0.65;0.6;I;7
+326;0.7;0.6;I;7
+327;0.75;0.6;I;7
+328;0.8;0.6;I;7
+329;0.85;0.6;IC;11
+330;0.9;0.6;IC;11
+331;0.95;0.6;IC;11
+332;1.0;0.6;IC;11
+333;1.05;0.6;IC;11
+334;1.1;0.6;IC;11
+335;1.15;0.6;IC;11
+336;1.2;0.6;I;7
+337;1.25;0.6;I;7
+338;0.01;0.65;;5
+339;0.05;0.65;;5
+340;0.1;0.65;;5
+341;0.15;0.65;;5
+342;0.2;0.65;;5
+343;0.25;0.65;;5
+344;0.3;0.65;;5
+345;0.35;0.65;;5
+346;0.4;0.65;;5
+347;0.45;0.65;;5
+348;0.5;0.65;;5
+349;0.55;0.65;I;7
+350;0.6;0.65;I;7
+351;0.65;0.65;I;7
+352;0.7;0.65;I;7
+353;0.75;0.65;I;7
+354;0.8;0.65;IC;11
+355;0.85;0.65;IC;11
+356;0.9;0.65;IC;11
+357;0.95;0.65;IC;11
+358;1.0;0.65;IC;11
+359;1.05;0.65;IC;11
+360;1.1;0.65;IC;11
+361;1.15;0.65;IC;11
+362;1.2;0.65;IC;11
+363;1.25;0.65;I;7
+364;0.01;0.7;;5
+365;0.05;0.7;;5
+366;0.1;0.7;;5
+367;0.15;0.7;;5
+368;0.2;0.7;;5
+369;0.25;0.7;;5
+370;0.3;0.7;;5
+371;0.35;0.7;;5
+372;0.4;0.7;;5
+373;0.45;0.7;;5
+374;0.5;0.7;;5
+375;0.55;0.7;I;7
+376;0.6;0.7;I;7
+377;0.65;0.7;I;7
+378;0.7;0.7;I;7
+379;0.75;0.7;IC;11
+380;0.8;0.7;IC;11
+381;0.85;0.7;IC;11
+382;0.9;0.7;IC;11
+383;0.95;0.7;IC;11
+384;1.0;0.7;IC;11
+385;1.05;0.7;IC;11
+386;1.1;0.7;IC;11
+387;1.15;0.7;IC;11
+388;1.2;0.7;IC;11
+389;1.25;0.7;IC;11
+390;0.01;0.75;;5
+391;0.05;0.75;;5
+392;0.1;0.75;;5
+393;0.15;0.75;;5
+394;0.2;0.75;;5
+395;0.25;0.75;;5
+396;0.3;0.75;;5
+397;0.35;0.75;;5
+398;0.4;0.75;;5
+399;0.45;0.75;;5
+400;0.5;0.75;;5
+401;0.55;0.75;;5
+402;0.6;0.75;I;7
+403;0.65;0.75;I;7
+404;0.7;0.75;I;7
+405;0.75;0.75;IC;11
+406;0.8;0.75;IC;11
+407;0.85;0.75;IC;11
+408;0.9;0.75;IC;11
+409;0.95;0.75;C;13
+410;1.0;0.75;C;13
+411;1.05;0.75;C;13
+412;1.1;0.75;IC;11
+413;1.15;0.75;IC;11
+414;1.2;0.75;IC;11
+415;1.25;0.75;IC;11
+416;0.01;0.8;;5
+417;0.05;0.8;;5
+418;0.1;0.8;;5
+419;0.15;0.8;;5
+420;0.2;0.8;;5
+421;0.25;0.8;;5
+422;0.3;0.8;;5
+423;0.35;0.8;;5
+424;0.4;0.8;;5
+425;0.45;0.8;;5
+426;0.5;0.8;;5
+427;0.55;0.8;;5
+428;0.6;0.8;I;7
+429;0.65;0.8;I;7
+430;0.7;0.8;I;7
+431;0.75;0.8;IC;11
+432;0.8;0.8;IC;11
+433;0.85;0.8;IC;11
+434;0.9;0.8;C;13
+435;0.95;0.8;C;13
+436;1.0;0.8;C;13
+437;1.05;0.8;C;13
+438;1.1;0.8;C;13
+439;1.15;0.8;IC;11
+440;1.2;0.8;IC;11
+441;1.25;0.8;IC;11
+442;0.01;0.85;;5
+443;0.05;0.85;;5
+444;0.1;0.85;;5
+445;0.15;0.85;;5
+446;0.2;0.85;;5
+447;0.25;0.85;;5
+448;0.3;0.85;;5
+449;0.35;0.85;;5
+450;0.4;0.85;;5
+451;0.45;0.85;;5
+452;0.5;0.85;;5
+453;0.55;0.85;;5
+454;0.6;0.85;;5
+455;0.65;0.85;I;7
+456;0.7;0.85;I;7
+457;0.75;0.85;IC;11
+458;0.8;0.85;IC;11
+459;0.85;0.85;IC;11
+460;0.9;0.85;C;13
+461;0.95;0.85;C;13
+462;1.0;0.85;C;13
+463;1.05;0.85;C;13
+464;1.1;0.85;C;13
+465;1.15;0.85;IC;11
+466;1.2;0.85;IC;11
+467;1.25;0.85;IC;11
+468;0.01;0.9;;5
+469;0.05;0.9;;5
+470;0.1;0.9;;5
+471;0.15;0.9;;5
+472;0.2;0.9;;5
+473;0.25;0.9;;5
+474;0.3;0.9;;5
+475;0.35;0.9;;5
+476;0.4;0.9;;5
+477;0.45;0.9;;5
+478;0.5;0.9;;5
+479;0.55;0.9;;5
+480;0.6;0.9;;5
+481;0.65;0.9;;5
+482;0.7;0.9;I;7
+483;0.75;0.9;IC;11
+484;0.8;0.9;IC;11
+485;0.85;0.9;IC;11
+486;0.9;0.9;C;13
+487;0.95;0.9;C;13
+488;1.0;0.9;C;13
+489;1.05;0.9;C;13
+490;1.1;0.9;C;13
+491;1.15;0.9;IC;11
+492;1.2;0.9;IC;11
+493;1.25;0.9;IC;11
+494;0.01;0.95;;5
+495;0.05;0.95;;5
+496;0.1;0.95;;5
+497;0.15;0.95;;5
+498;0.2;0.95;;5
+499;0.25;0.95;;5
+500;0.3;0.95;;5
+501;0.35;0.95;;5
+502;0.4;0.95;;5
+503;0.45;0.95;;5
+504;0.5;0.95;;5
+505;0.55;0.95;;5
+506;0.6;0.95;;5
+507;0.65;0.95;;5
+508;0.7;0.95;I;7
+509;0.75;0.95;IC;11
+510;0.8;0.95;IC;11
+511;0.85;0.95;IC;11
+512;0.9;0.95;C;13
+513;0.95;0.95;C;13
+514;1.0;0.95;C;13
+515;1.05;0.95;C;13
+516;1.1;0.95;C;13
+517;1.15;0.95;IC;11
+518;1.2;0.95;IC;11
+519;1.25;0.95;;5
+520;0.01;1.0;;5
+521;0.05;1.0;;5
+522;0.1;1.0;;5
+523;0.15;1.0;;5
+524;0.2;1.0;;5
+525;0.25;1.0;;5
+526;0.3;1.0;;5
+527;0.35;1.0;;5
+528;0.4;1.0;;5
+529;0.45;1.0;;5
+530;0.5;1.0;;5
+531;0.55;1.0;;5
+532;0.6;1.0;;5
+533;0.65;1.0;;5
+534;0.7;1.0;;5
+535;0.75;1.0;IC;11
+536;0.8;1.0;IC;11
+537;0.85;1.0;IC;11
+538;0.9;1.0;C;13
+539;0.95;1.0;C;13
+540;1.0;1.0;C;13
+541;1.05;1.0;C;13
+542;1.1;1.0;C;13
+543;1.15;1.0;IC;11
+544;1.2;1.0;;5
+545;1.25;1.0;;5
+546;0.01;1.05;;5
+547;0.05;1.05;;5
+548;0.1;1.05;;5
+549;0.15;1.05;;5
+550;0.2;1.05;;5
+551;0.25;1.05;;5
+552;0.3;1.05;;5
+553;0.35;1.05;;5
+554;0.4;1.05;;5
+555;0.45;1.05;;5
+556;0.5;1.05;;5
+557;0.55;1.05;;5
+558;0.6;1.05;;5
+559;0.65;1.05;;5
+560;0.7;1.05;;5
+561;0.75;1.05;IC;11
+562;0.8;1.05;IC;11
+563;0.85;1.05;IC;11
+564;0.9;1.05;C;13
+565;0.95;1.05;C;13
+566;1.0;1.05;C;13
+567;1.05;1.05;C;13
+568;1.1;1.05;C;13
+569;1.15;1.05;IC;11
+570;1.2;1.05;;5
+571;1.25;1.05;;5
+572;0.01;1.1;;5
+573;0.05;1.1;;5
+574;0.1;1.1;;5
+575;0.15;1.1;;5
+576;0.2;1.1;;5
+577;0.25;1.1;;5
+578;0.3;1.1;;5
+579;0.35;1.1;;5
+580;0.4;1.1;;5
+581;0.45;1.1;;5
+582;0.5;1.1;;5
+583;0.55;1.1;;5
+584;0.6;1.1;;5
+585;0.65;1.1;;5
+586;0.7;1.1;;5
+587;0.75;1.1;;5
+588;0.8;1.1;IC;11
+589;0.85;1.1;IC;11
+590;0.9;1.1;C;13
+591;0.95;1.1;C;13
+592;1.0;1.1;C;13
+593;1.05;1.1;C;13
+594;1.1;1.1;C;13
+595;1.15;1.1;IC;11
+596;1.2;1.1;;5
+597;1.25;1.1;;5
-- 
GitLab