Remarkably, the reconstruction's mean error (or total volume error) is exceptionally low, i.e. by several orders of magnitude lower than the absolute mean error for each sampling point. The higher mean absolute error is also a result of the somewhat “noisy” input data. As the samples are the results of a stochastic raytracing process, they are subject to a certain degree of statistical errors (see 5.1.3). However, since the (signed) mean error (or total volume error) is not an explicit optimisation criterion, the

excellent match regarding this measure indicates that the spherical harmonics can capture the geometric nature of the underlying process rather than only reconstruct the individual samples. The high data compression rate also supports this interpretation: the information of each case, represented by 1440 sampling points, can be reconstructed by summation of less than 90 spherical harmonics, as shown in the first row of Table 11.

The polar plot representations for the sampling data as well as for the final functional SIOP, are depicted in Figure 81. It can be seen how the final shape of the SIOP matches the sampled data. The polar plots presented in Figure 82 reveal the directional dependencies of the different shading states. The unshaded state is open to a broader angular range with a peak at normal incidence. Note that the ideal case of a simple planar absorber, governed only by the cosine law, would result in a perfectly spherical SIOP. However, due to the shading effects of the window reveal and the angular characteristics of the glazing, the SIOP of the unshaded window detail shows a rather balloon-like shape with smaller values for higher angles of incidence.

<figure><img src="image_86.png" /><figcaption>Figure 81 Polar plot representation of the sampled data (right) vs. the fitted SIOP (left); both blended (center); (example transmitted power for case 30° tilted slats)</figcaption></figure>

<figure><img src="image_87.png" /><figcaption>Figure 82 Polar plot representation of generated SIOPs for shaded window (unshaded; shaded 0°; shaded 30°; shaded 60°)</figcaption></figure>

Once the shading is closed, maintaining horizontally oriented slats, the SIOP is squeezed towards a lentil-like shape, as predominantly light from relatively horizontal directions can pass the blinds. Compared to the unshaded state, the hemispherical mean value of transmittance is reduced to 57% for the horizontal (0°) slat angle, as can be seen in the polar plot by the decrease in size. Tilting the slats further reduces the hemispherical mean value to 24% (for 30° tilt) and 6%, respectively (for 60° tilt). It can also be seen how the direction of maximum transmission is tilted along with the slat angle.
## 6.2. SIOP evaluation and incidence profiles

Evaluating the SIOP means to determine the quantity of interest $q$ represented by the SIOP for a given solar incidence profile. Based on the definition of the SIOP only a simple integration over the hemisphere is required for this purpose. The solar incidence profile will usually contain radiance information but could also contain time-integrated values for the radiance of unit $Jm^{-2}sr^{-1}$. In this case, the evaluation of the SIOP would also yield time-integrated quantities (e.g. energy instead of power). The equations derived below are valid for any type of incidence profile; however, for the sake of simplicity, the physical quantity radiance $L(\theta, \varphi)$ [W/(m².sr)] is used for the incidence profile. The evaluation of the target quantity $q$ can then be written as:

$$q = \iint_{\text{hemi}} (\text{SIOP}_q \cdot I) d\Omega = \int_{\theta=0}^{\pi/2} \int_{\varphi=0}^{2\pi} \text{SIOP}_q(\theta, \varphi) \cdot L(\theta, \varphi) \sin(\theta) d\varphi d\theta \quad (109)$$

In order to perform the integration, the incidence profile in world coordinates has to be rotated into the local SIOP coordinate system according to the element's actual orientation (e.g. a window) first (or vice versa).

In practice, a sampling method is used to evaluate the integral. As equally-spaced, discrete directions are used, the solid angle can be pulled out of the integration, which results in an elementwise multiplication. Consequently, it is possible to state the discrete integration over the entire hemisphere as:

$$q = \iint_{\text{hemi}} (\text{SIOP}_q \cdot L) d\Omega = \frac{2\pi}{N} \sum_{i=0}^{N-1} \text{SIOP}_q(i) \cdot L(i) \quad (110)$$

In more general terms, the integration over any segment of the hemisphere can be formulated as:

$$q = \iint_{\text{seg}} (\text{SIOP}_q \cdot L) d\Omega = \frac{\Omega_{\text{seg}}}{N} \sum_{j=0}^{M-1} \text{SIOP}_q(j) \cdot L(j) \quad (111)$$

Again, it is required that the directions $j$ are distributed equally spaced in the relevant segment covering the evaluated solid angle $\Omega_{seg}$.

Rarely the angular radiance distribution on the entire hemisphere is available. Instead, a generic incidence profile can be generated based on the commonly used direct and diffuse irradiance information. In such cases, the angular incidence profile is not explicitly calculated but implicitly contained in the evaluation process. In commonly used irradiance models, the diffuse radiance is assumed constant on specific angular segments; e.g. the most simple isotropic model assumes constant radiance values for the diffuse radiation from the sky and the ground-reflected radiation. In this case, the radiance is constant in the according integration areas but generally a function of time: $I(t, j) = I_{0,seg}(t)$. This allows to additionally pull the radiance out of the integration:

$$q(t) = \Omega_{\text{seg}} \cdot L_{0,\text{seg}}(t) \left( \frac{1}{N} \sum_{j=0}^{M-1} \text{SIOP}_q(j) \right) \quad (112)$$

Since the term in the brackets reflects the average value $\bar{S}_{q,seg}$ of the SIOP$_q$ on a specific segment of the hemisphere, the evaluation can be written as:

$$q(t) = \Omega_{\text{seg}} \cdot \bar{S}_{q,\text{seg}} \cdot L_{0,\text{seg}}(t) \quad (113)$$

Utilizing this formulation for the evaluation has a significant advantage: in order to achieve high precision, the evaluation of the SIOP should be performed with a large number of directional samples, which requires quite some calculation time as trigonometric functions are involved. However, if the
radiance value is a constant for a specific fixed segment of the hemisphere, the value of $\bar{S}_{q,seg}$ can be precalculated before the evaluation of the time series. Since this calculation must be performed only once, it is no longer time-critical. Consequently, any appropriate number of samples and any appropriate distribution of directions can be used for the numerical determination of the constants $\bar{S}_{q,seg}$. The subsequent evaluation of irradiance time series can now be performed very fast, as it merely consists of a couple of multiplications for each diffuse component. Only for the direct solar component, an evaluation of the SIOP has to be performed on every time step. However, this can now be performed for the precise location of the sun and, therefore, requires only a single evaluation. The process is explained in detail in the following sections.

### 6.2.1. Evaluation with angularly resolved measured data

For the case that measured data is available with high resolution for the entire hemisphere, the full evaluation according to equations (109) or (110) can be performed using a Monte Carlo method (see section 5.1). In practice, such angular-resolved data is more likely to be available with a lower resolution (e.g. based on the Klems segmentation). In such cases, the SIOP can be evaluated by applying equation (113) for each segment that represents a locally constant (modelled or measured) radiance value. The necessary coefficients, i.e. the solid angle $\Omega_{seg}$ and the SIOP average value for the segment $\bar{S}_{q,seg}$ can again be precomputed based on the given segmentation. This allows a fast evaluation afterwards.

### 6.2.2. Evaluation using a generic isotropic incidence profile

In the simple but commonly used isotropic sky model (see Hottel et al. (1942) or Loutzenhiser et al. (2007)), the solar irradiance is modelled based on three commonly used, measured (or modelled) components direct normal incidence DNI, global horizontal irradiance GHI and diffuse horizontal irradiance DHI. The model assumes that the radiance over the sky dome is constant in all directions. The model further assumes the radiance resulting from ground-reflected solar radiation is also a constant and can be modelled as a fraction of the global horizontal irradiance GHI. The relevant coefficient that determines this fraction is referred to as albedo value $\alpha$. It is usually a constant representing the ground's bi-hemispherical reflectance value.

Since the model is based on the assumption of perfectly diffuse Lambert behaviour (see section 4.3.3), the radiances for the sky and ground component of the horizontally measured irradiance quantities can be derived as:

$$L_{sky}(t) = \frac{DHI(t)}{\pi} \qquad L_{gnd}(t) = \frac{albedo \cdot GHI(t)}{\pi} \qquad (114)$$

Based on these definitions, the time-series $q(t)$ for any quantity of interest represented as SIOP can easily be evaluated by applying equation (113). In addition to the two diffuse components, the direct normal irradiance is considered in a third term that evaluates the SIOP for the solar angles $\theta(t), \varphi(t)$ (relative to the SIOP orientation):

$$q_{dir}(t) = SIOP_q(\theta(t), \varphi(t)) \cdot DNI(t) \qquad (115)$$

Considering the two diffuse and the direct component, the evaluation can be written as:

$$q(t) = \Omega_{sky} \cdot \bar{S}_{sky} \cdot L_{sky}(t) + \Omega_{gnd} \cdot \bar{S}_{gnd} \cdot L_{gnd}(t) + SIOP_q(\theta(t), \varphi(t)) \cdot DNI(t) \qquad (116)$$
The coefficients $\Omega$ and $\bar{S}$, relevant for the diffuse components, can again be precomputed, provided that the SIOP reference system is first rotated into the correct world position.

In the following, it is demonstrated how the RadiCal approach complies with the standard method by calculating the total irradiance incident on a vertically tilted plane using the SIOP approach. Assuming a planar, totally absorbing surface, the cosine law leads to an absorption SIOP that has the shape of a perfect sphere, extending to a value of 1.0 in the normal direction. It is crucial for further derivations and is hereinafter referred to as Lambert-SIOP. Its polar-plot representation is depicted in Figure 83. Following geometrical considerations, it is obvious that the pre-computed coefficients for the determination of the diffuse components are:

$$ \bar{S}_{sky} = 0.5, \bar{S}_{gnd} = 0.5, \Omega_{sky} = \pi sr, \Omega_{gnd} = \pi sr $$

Substituting these coefficients into equation (116), the time-dependent total irradiance on the plane $E(t)$ can be written as:

$$ E(t) = \frac{\pi}{2} \cdot L_{sky}(t) + \frac{\pi}{2} \cdot L_{gnd}(t) + SIOP(\theta(t), \varphi(t)) \cdot DNI(t) \quad (117) $$

By substituting the above radiances using eq. (114), and by utilizing the fact that the Lambertian, sphere-shaped SIOP (see Figure 83) can also be expressed simply by $cos(\gamma)$, with $\gamma$ being the angle of incidence relative to the plane, equation (117) can be written as:

$$ E(t) = \frac{DHI(t)}{2} + \frac{\alpha GHI(t)}{2} + cos(\gamma(t)) \cdot DNI(t) \quad \text{Q.E.D.} \quad (118) $$

This is the standard equation used to determine the irradiance based on the isotropic sky model for a plane tilted by 90°, which has been derived consistently using the SIOP approach. Of course, the same could be shown for any other tilt angle of the plane, as the coefficients for the sky and ground term determined on a tilted SIOP will change accordingly.

<figure><img src="image_88.png" /><figcaption>Figure 83 Lambert-SIOP - sphere shaped SIOP representing absorption by a Lambertian plane (polar plot representation)</figcaption></figure>
### 6.2.3. Evaluation using a generic Perez incidence profile

#### 6.2.3.1. The anisotropic Perez sky model

While the isotropic sky model is still frequently used, the assumption of constant radiance over the entire sky dome does not properly reflect the real-world irradiance profiles, leading to erroneous results. The actual radiance distribution over the sky dome and its dynamic diurnal and seasonal variation is complex. It depends on specific local meteorological conditions, in particular, of course, on the weather and cloud cover. Therefore, it is impossible to accurately describe a local irradiance situation by applying simple analytic models. Nonetheless, the isotropic model can be 'upgraded' by

adding additional features. In order to achieve higher accuracy, Perez et al. (1988; 1990) have identified two dominant non-isotropic components that can be modelled in an analytic way: circumsolar radiation and horizon brightening (or darkening). To be able to take local variations into account, the models are based on free parameters that can specifically be defined based on local conditions. The Perez model can, therefore, be considered a semi-empirical approach.

While the direct normal irradiance

and ground reflected irradiance is treated in the same way as in the isotropic sky model, the diffuse sky irradiance received by a tilted surface is now determined by (Perez et al., 1988; Perez et al., 1990):

<figure><img src="image_89.png" /><figcaption>Figure 84 The five irradiance components considered by the Perez model</figcaption></figure>

$$E_{diff,sky} = DHI \cdot \left[ (1 - F_1) \left( \frac{1 + \cos\beta}{2} \right) + F_1 \frac{a}{b} + F_2 \sin\beta \right] \quad (119)$$
with $a = \max(0, \cos \gamma)$ and $b = \max(\cos 85^\circ, \cos \theta_z)$

$\beta$ ... tilt angle of plane (0°=horiz. plane)
$F_1$, $F_2$ ... "Perez coefficients"

$\gamma$ ... incidence angle
$\theta_z$ ... solar zenith angle

For the application in RadiCal, it is essential to identify and separate the components of the Perez model: The first term in the square brackets of equation (119) represents the isotropic diffuse irradiance of the sky, excluding the circumsolar component. The second term, led by the coefficient $F_1$, reflects the circumsolar component of the model, whereas the last term, led by $F_2$, represents the horizon component of the model.

The coefficients $F_1$ and $F_2$ are, in turn, determined based on the irradiance components DHI, DNI, the solar zenith angle $\theta_z$, as well as on the local pressure (for air mass calculation) and date, i.e. day of the year (for calculation of extraterrestrial radiation). Additionally, empirical parameters $f_{ij}$, specific to the location, are required in the calculation process. The implementation performed here follows the often applied approach to use global values as provided in table 6 in the publication of Perez et al. (1990). The relevant parameters $f_{ij}$ are provided in tabular form depending on the clearness parameter $\epsilon$, which depends on diffuse/direct irradiance ratio and the zenith angle.
### 6.2.3.2. SIOP evaluation using the anisotropic Perez sky model

The necessary details for the implementation of the Perez model are available in the original publications (Perez et al., 1988; Perez et al., 1990) or, e.g., on the relevant page of the *PV Performance Modeling Collaborative of Sandia National Laboratories* (Sandia National Laboratories, 2018b). The application of the model in *RadiCal* basically follows these standards with two variations: first, the determination of the coefficients $f_{ij}$ based on the clearness parameter $\epsilon$ is not performed based on the discrete interval approach, as suggested. Instead, a linear interpolation method is applied (see implementation, section 6.2.5). Second, the three components forming the diffuse sky irradiance (eq. (119)) are calculated and applied separately. Based on Perez model radiance and irradiance components are defined subsequently used for the *SIOP* evaluation.

By following the principle of a Lambert receiver or by comparing with the equations in the previous section, the constant sky radiance value $L_{sky,Perez}$ can be derived from the Perez model as:

$$L_{sky,Perez} = \frac{DHI}{\pi} \cdot (1 - F_1) \qquad (120)$$

Note that this definition still considers the entire sky dome, i.e. the upper hemisphere. The horizon-band introduced by Perez does not reduce its extension, as the horizon radiance values are superimposed on the sky radiance value. Since these horizon radiance values can be positive and negative, they can either increase or reduce the original sky radiance values.

For the ground-reflected component, commonly, the standard approach for the isotropic model is used. It was introduced in the previous section:

$$L_{gnd} = \frac{albedo \cdot GHI}{\pi} \qquad (114)$$

The Perez model assumes a horizon band that extends to an angle of 6.5° above the horizon. To be able to perform an evaluation analogue to the two diffuse components just introduced, a constant radiance is assumed for this region. The corresponding value for the related radiance can be found by inverting the Perez model, i.e. demanding that the irradiance evaluation of the *Lambert-SIOP*, representing a vertical plane, according to equation (113), has to equal the horizon irradiance component of the Perez model:

$$E_{horizon} = DHI \cdot F_2 \sin 90^{\circ} = \Omega_{0,horizon} \cdot \bar{S}_{0,horizon} \cdot L_{0,horizon} \qquad (121)$$

Consequently, the constant radiance value of the horizon $L_{0,horizon}$ is obtained as:

$$L_{0,horizon} = \frac{DHI \cdot F_2}{\Omega_{0,horizon} \bar{S}_{0,horizon}} = \frac{DHI \cdot F_2}{0.3556 \cdot 0.6353} = DHI \cdot F_2 \cdot 4.426 \qquad (122)$$

The solid angle and average value of the *Lambert-SIOP* used represent constants. The value for the solid angle $\Omega_{0,horizon}$ is achieved simply by integration over the 6.5 degrees extended horizon band, considering only the hemisphere in front of the plane. The average value $\bar{S}_{0,horizon}$ is obtained by evaluating the *SIOP* for this same region. The determined constants ensure that the irradiance values resulting from the horizon brightening (darkening) align with the original definitions of the Perez model. The actual solid angles and average values are used in the evaluation process to perform a *SIOP*-based evaluation of the Perez model for a plane tilted at angles other than 90° (see below). However, the constants contained in equation (122) always remain unchanged.

Next, the second anisotropic component introduced by Perez *et al.* is considered. In their definition, the circumsolar component represents a part of the diffuse sky radiation. This is reasonable since it
scales with the diffuse irradiance *DHI*. However, for its evaluation, the sky's clearness, which is in turn linked to the direct normal component *DNI*, is also relevant. Further, to calculate the irradiance resulting from the circumsolar component, Perez introduced the parameter *a*, which by definition is nothing but the non-negative cosine of the incidence angle (*max(0, cos γ)*). Hence, it is reasonable to consider the circumsolar component as directed and to handle it in the same way as the *DNI* component. Consequently, it is not represented as a radiance term [W/(m².sr)] but as an irradiance component [W/m²]. The related term can easily be found by replacing the parameter *a* with the *SIOP* for the incidence direction, as the *SIOP* operator contains the angular profile in a more general way:

$$E_{circ,Perez} = SIOP_q(\theta, \varphi) \cdot DHI \cdot F_1 \cdot \frac{1}{b} \qquad (123)$$

Alternatively, the *SIOP* operator could also be evaluated on the entire circumsolar disc as defined in Perez' model. In future work, it will be studied if this improves the accuracy. However, the simple approach based on a singular *SIOP* evaluation for the circumsolar and direct component is pursued for all evaluations performed in the present work for two reasons. First, it is more in line with the commonly performed evaluation of the Perez model. Second, the single evaluation of the *SIOP* for the incidence angle will generally be a good estimator for the average value of the circumsolar disk centred around it.

Finally, only the direct component is missing. In analogy to the circumsolar term just derived, the irradiance resulting from direct normal irradiance *DNI* is stated as:

$$E_{dir} = SIOP_q(\theta, \varphi) \cdot DNI \qquad (124)$$

The general form for the evaluation of a *SIOP* using the Perez sky model can now be written as:

$$\begin{aligned} E(t) = & \Omega_{sky} \cdot \bar{S}_{sky} \cdot L_{sky,Perez}(t) \\ & + \Omega_{gnd} \cdot \bar{S}_{gnd} \cdot L_{gnd}(t) \\ & + \Omega_{horizon} \cdot \bar{S}_{horizon} \cdot L_{horizon,Perez}(t) \\ & + E_{circ,Perez} \\ & + E_{dir} \end{aligned} \qquad (125)$$
### 6.2.3.3. Irradiance on the vertical plane - consistency with standard approach

In analogy to the derivation performed for the isotropic model in the last section, the consistency of the Perez evaluation for the SIOP is demonstrated by explicitly calculating the irradiance on a vertically oriented plane.

The solid angles for this orientation represent quarter spheres for the sky and ground component: $\Omega_{sky} = \Omega_{gnd} = \pi$. The 6.5° extended horizon results in: $\Omega_{horizon} \approx 0.3556$. The average values over the related SIOP segments can be determined as $\bar{S}_{sky} = \bar{S}_{gnd} = 0.5$ and $\bar{S}_{horizon} \approx 0.6353$.

Substituting these values into the previous equation:

$$ 
\begin{aligned}
E(t) = & \pi \cdot 0.5 \cdot L_{sky,Perez}(t) \\
& + \pi \cdot 0.5 \cdot L_{gnd}(t) \\
& + 0.3556 \cdot 0.6353 \cdot L_{horizon,Perez}(t) \\
& + SIOP_q(\theta(t), \varphi(t)) \cdot \left[ DNI(t) + DHI(t) \cdot \frac{F_1(t)}{b(\theta(t))} \right]
\end{aligned} 
\quad (126) $$

For the irradiance on a plane, the evaluation of the Lambert-SIOP for the solar position is again substituted by $cos(\gamma)$. Further, equations (120) to (124) are substituted into the equation:

$$ 
\begin{aligned}
E(t) = & \pi \cdot 0.5 \cdot \frac{DHI(t)}{\pi} \cdot (1 - F_1(t)) \\
& + \pi \cdot 0.5 \cdot \frac{albedo \cdot GHI(t)}{\pi} \\
& + 0.3556 \cdot 0.6353 \cdot \frac{DHI}{0.3556} \cdot \frac{F_2(t)}{0.6353} \\
& + cos(\gamma(t)) \cdot DNI(t) \\
& + DHI(t) \cdot \frac{F_1(t)}{b(\theta(t))} \cdot cos(\gamma(t))
\end{aligned} 
\quad (127) $$

$$ = (1 - F_1(t)) DHI(t) \cdot \frac{1}{2} \triangleq \text{Perez sky on vert. plane} $$

$$ + F_1(t) \cdot DHI(t) \cdot \frac{a(\gamma(t))}{b(\theta(t))} \triangleq \text{Perez circum-solar on vert. plane} $$

$$ + F_2(t) \cdot DHI(t) \triangleq \text{Perez horizon on vert. plane} $$

$$ + cos(\gamma(t)) \cdot DNI(t) \triangleq \text{std. direct normal on vert. plane} $$

$$ + \frac{albedo \cdot GHI(t)}{2} \triangleq \text{std. ground reflection on vert. plane} $$

This, again, demonstrates that the SIOP approach transitions into the standard model for a simple planar geometry. The evaluation of equation (125) for any other orientation simply requires the (pre)calculation of the solid angles and SIOP averages applying the actual tilt angle.
### 6.2.3.4. Perez evaluation for a general tilt angle

In order to perform the *SIOP* evaluation for the Perez model according to equation (127), the three solid-angle and the three *SIOP* average values have to be provided. All values are specific for the tilt-angle and do not depend on the azimuthal, i.e. compass orientation. The solid angles represent the fractional area of the particular regions (sky, ground, horizon band) visible to an observer on the tilted plane. Following this consideration, it is evident that the areas representing the sky and ground are spherical lunes. Hence the solid angles of the visible sky ($\Omega_{sky}$) and ground ($\Omega_{gnd}$) are simply linear functions ranging from 0 to $2\pi$ (see Figure 85). The solid-angle of the visible horizon band ($\Omega_{horizon}$) follows a slightly more complex function. For most tilt angles, it is close to 0.3556sr as one half of the horizon will be visible. However, there is a slight slope, as additional wedges around the axis of rotation have to be considered. When the downward-pointing horizontal orientation is reached, the solid angle will diminish to zero as the horizon band moves out of view. In contrast, for the upright horizontal plane, the value reaches 0.7113sr, as the entire 6.5° extended band is visible.

While the solid angles used in the evaluation depend on the tilt angle but are generally valid for any *SIOPs*, the corresponding *SIOP* average values for the specific regions naturally depend on the actual *SIOP* shape. For reference and comparison with the standard model, the *Lambert-SIOP*, representing the irradiance on a plane, is once more evaluated. The *SIOP* averages as functions of tilt angle are determined by numerical integration and depicted in Figure 86. Note that the average for the horizon band ($\bar{S}_{horizon}$) reaches its maximum when tilted slightly upwards (red arrow) and does not reach zero for the horizontal orientation (red line). In order to analyse the impact that the radiance of each component has on the resulting irradiance on the plane, the solid angles are multiplied by the relevant *SIOP* averages for each component. The results are depicted in Figure 87. As can be seen, the products for the sky and ground components result in cosine- and sine-shaped functions of the

<figure><img src="image_90.png" /><figcaption>Figure 85 Relevant solid-angles as function of tilt-angle (0=vertical plane)</figcaption></figure>

<figure><img src="image_91.png" /><figcaption>Figure 86 Cosine-SIOP averages as function of tilt-angle</figcaption></figure>

<figure><img src="image_92.png" /><figcaption>Figure 87 Products of solid-angles and SIOP averages as function of tilt-angle</figcaption></figure>
tilt angle. This corresponds perfectly to the standard Perez model, eq. (119). However, the irradiance resulting from the horizon radiance, determined by the product $\Omega_{horizon} \cdot \bar{S}_{horizon}$, follows a slightly different shape than the $\cos\beta$ term considered in equation (119). Firstly, the maximum is offset by few degrees to the upside, which is reasonable based on geometrical considerations. Second, the impact of the horizon band radiance does not vanish for the (upright) horizontal orientation. From a geometric point of view, this is also reasonable. As stated above, the visible solid angle fraction of the horizon band increases to double the value for a horizontal plane. Even if the absorbed irradiance converges to zero for incoming light at grazing angles (according to Lambert's cosine law), there remains a non-zero contribution of the 6.5° extended horizon band on a horizontal plane. This is an interesting deviation of the SIOP-based Perez evaluation compared to the standard model. Even if the irradiance contribution of the horizon band might be limited for a planar surface, SIOPs can generally be used to model complex structures that could, intentionally or not, absorb or redirect light coming from these directions. For such cases, it is essential to accurately consider the impact of the radiance resulting from the horizon band.

<figure><img src="image_93.png" /><figcaption>Figure 88 Perez evaluation of SIOP orientation in world-coordinate for different tilt angles (0°, 30°, 90° vs. vertical) and compass orientation 220° (SSW)</figcaption></figure>
### 6.2.4. General approach for real-world environments

The evaluations discussed in the previous two sections assume that the view to the sky, horizon and ground is unobstructed. This represents a base case that is fundamentally important but rarely relevant for practical applications.

Following theoretical considerations, as well as the results of the full system validation performed in this work (Chapter 8), it is evident that the constant albedo approach leads to significant inaccuracies. It assumes a ground-level horizontal plane with constant Lambert reflectance (section 4.3.3). In real environments, all three assumptions generally prove to be insufficient. First, relevant surfaces in the environment generally have varying orientations and rise significantly above ground level, e.g. facades, mountains or tilted roofs. Consequently, irradiance values vary considerably for each surface, depending on its orientation and shading by other objects. Second, many surfaces exhibit significant amounts of specular reflection, e.g. metallic roofs. Finally, the reflectance will change depending on the climate and weather, e.g. vegetation surfaces change their colour and precipitation like rain or especially snow will significantly alter the reflectance of surfaces.

For practical application, individual radiance components reflected off all visible surfaces in the target's view have to be considered. In order to do this, the geometry of these objects is mapped onto the SIOP, i.e. a segmentation of the SIOP is performed in a way that the radiance of each segment is approximately constant.

In such cases, the SIOP can be evaluated by applying equation (113) for each segment representing a locally constant (modelled or measured) radiance value. The necessary coefficients, i.e. the solid angle $\Omega_i$ and the SIOP average value for the segment $\bar{S}_i$ can again be precomputed based on the given segmentation. Finally, the solid angles and SIOP average values are calculated for the remaining visible regions of the sky and horizon band. These regions additionally determine if the sun is visible at its current location, i.e. if direct irradiance must be considered.

If the Perez sky model is applied for this approach, the total irradiance can be written as:

$$ 
\begin{aligned}
E(t) = & \sum_{i=1}^{N} \Omega_i \cdot \bar{S}_i \cdot L_i(t) \\
& + \Omega'_{\text{sky}} \cdot \bar{S}'_{\text{sky}} \cdot L_{\text{sky,Perez}}(t) \\
& + \Omega'_{\text{horizon}} \cdot \bar{S}'_{\text{horizon}} \cdot L_{\text{horizon,Perez}}(t) \\
& + \Omega'_{\text{circ}} \cdot \bar{S}'_{\text{circ}} \cdot L_{\text{circ,Perez}}(t) \\
& + \delta_{\text{sun}} E_{\text{dir}}
\end{aligned} 
\quad (128) $$

$\Omega_i$ ... solid angle of segment i

$\bar{S}_i$ ... SIOP average of segment i

$L_i(t)$ ... locally constant radiance of segment i

$\Omega'_{\text{sky}}, \Omega'_{\text{horizon}}, \Omega'_{\text{circ}}$ ... solid angles of remaining unobstructed regions

$\bar{S}'_{\text{sky}}, \bar{S}'_{\text{horizon}}, \bar{S}'_{\text{circ}}$ ... SIOP averages of remaining unobstructed regions

$\delta_{\text{sun}}$ ... delta-function: 1...sun visible, 0...view on sun obstructed
<figure><img src="image_94.png" /><figcaption>Figure 89 projection of radiance patches on located and oriented SIOP (here the Lambert-SIOP for irradiance on plane is depicted)</figcaption></figure>

While the segmentation and evaluation of the *SIOP* based on the locally visible objects is a straight-forward task, assigning suitable radiance values $L_i(t)$ to the different patches is more challenging. The time-resolved radiance values for each visible surface can be determined at various levels of detail.

Some potential approaches, ordered by increasing complexity, are:

* consideration of sky, horizon and ground radiance only, but with an arbitrarily shaped horizon line

* consideration of patches with varying albedo values (diffuse reflection, horizontal plane)

* consideration of plane irradiance on patches (diffuse reflection, oriented plane)

* consideration of plane irradiance on patches (diffuse reflection, oriented plane, shadowing)

* consideration of plane irradiance and directed reflection (direct reflection, diffuse reflection, oriented plane, shadowing)

* elaborate calculation of radiances using *RadiCal* raytracing and an adapted *SIOP* approach for the patches.

An efficient approach is currently being implemented and tested but cannot be detailed here as this would go beyond the scope of the PhD. This work focuses on the raytracing algorithm, physical optics, *SIOP* definition, and their validation. Future publications will, in particular, cover the handling of spatially resolved incidence profiles.

However, it is important to stress that the *SIOP* concept does not impose any restrictions regarding the evaluation of a detailed radiance environment. On the contrary, *SIOPs* in particular allow and facilitate the implementation of such angularly resolved approaches. The significant advantage of the *SIOP* is that radiance distributions on an arbitrarily shaped segmentation can be evaluated very accurately and efficiently.
### 6.2.5. Implementation Perez module

The program code relevant for the Perez evaluation are contained in the, already introduced, classes `TSunposIrradiance` and `TRCSIOP`. The latter class contains the algorithm for calculating the solid angles and `SIOP` average values. Both tasks are performed in the function `precalcHemiMeansStdSYS`:

```c
procedure TRCSIOPSamples.precalcHemiMeansStdSYS(const relPrecTarget: Double);
const NperSample = 1000;
var
  vSum, v2Sum: Double;
  sCnt: Integer;
  sdev, sem: double;
  mn: Double;
  costh: double;
  phi: double;
  c, ulc: Integer;
  vd: Tvect3d;
  fc: Double;
  oldAzi: Double;
  fsHit: Boolean;
  missed: Integer;
  hitN: Integer;
  hSum: Double;
  nrm: Double;

const sin6_5deg = 0.11320321; // height of Perez horizon band (6.5°)
const horizon360bandSolidAngleDivPI = 0.226406428;
// = 0.711276769 sr / pi = 2*PI() * SIN(6,5/180*PI()) / pi = sin(6.5°)*2

procedure AddSampleStatsCalcSEM;
  Begin
    if hitN=0 then mn:=0 else mn:=mn/hitN;
    Inc(scnt);
    v2sum:=v2sum+sqr(mn); vSum:=vSum+mn;
    hSum:=hSum+(hitN/NPerSample);
    sDev:=sqrt( abs(v2sum/sCnt - sqr(vSum/sCnt) ) ); // std. dev of sample means
    SEM:=sdev/sqrt(sCnt); // standard error of sample means
    if mn=0 then SEM:=0 else SEM:=abs(SEM/mn); // relative SEM
  end;

begin
  oldAzi:=azi;
  initSIOPorientation(tilt, SOUTH); // make sure orientatin is set
                                   // temporarily point SIOP to SOUTH
  for ulc:=0 to 2 do
    Begin
      sCnt:=0;
      vSum:=0;
      hSum:=0;
      v2Sum:=0;
      Repeat
        mn:=0;
        missed:=0;
        For c:=1 to NPerSample do
          Begin
            if ulc=0 then costh:=+random else // upper hemi/sky
              if ulc=1 then costh:=-random // lower hemi/gnd
              else costh:=+random*sin6_5deg; // 6.5° horizon band Perez
            // equally sample the hemisphere or horizon band
            vd.z:=costh; // z is uniformly distributed
            fc:=sqrt(1-sqr(costh));
            // phi:=random*2*pi; // x,y unused slower trigonometric version
            // vd.x:=sin phi*fc; vd.y:=cos phi*fc;

            Repeat // x,y faster version
              vd.x:=2*random-1; // rejection sampling circle in unit box
              vd.y:=2*random-1;
              nrm:=sqr(vd.x)+sqr(vd.y);
            until (nrm>0) and (nrm<=1);
            nrm:=sqrt(nrm);
            vd.x:=vd.x/nrm*fc;
            vd.y:=vd.y/nrm*fc;

            mn:=mn+evalSIOPWorld(vd, fsHit);
            if fsHit=false then Inc(missed);
          End;

          hitN:=(NPerSample-missed);

          AddSampleStatsCalcSEM; // calculate realtive standard error of means
      until ((sCnt>500) and (SEM<relPrecTarget)) or (sCnt>100000);
```
```javascript
// integrate until precision reached or em.brake
if ulc=0 then ssiop.meanUpperHemi:=vSum/sCnt else
if ulc=1 then ssiop.meanLowerHemi:=vSum/sCnt
else Begin
    ssiop.meanHorizBand:=vSum/sCnt;
    // numerically get solid angle for vis. horizon
    ssiop.solidAngleDivPiHorizBand:=(hSum/sCnt)*horizon360bandSolidAngleDivPI;
end;

// analytic solid angle for visible sky and ground lunes:
SSIOP orientation (Tilt, oldAzi);

// analytic solid angle for visible sky and ground lunes:
SSIOP orientation (Tilt, oldAzi);

end;
```

The function samples the average SIOP values for the relevant sections of the sphere (sky, ground, horizon band) until a given precision target is met. The central limit theorem is again applied to abort the sampling at the earliest possible iteration (see section 5.1.3). In order to sample the sphere's sections uniformly, again rejection sampling, referred to as 'Malley's method' (Malley, 1988) is applied (compare section 4.11.2). This algorithm is approx. 30% faster than the standard code, which involves the evaluation of sine and cosine functions. Note that all solid angles are kept as dimensionless fractions (of quarter-spheres) by dividing them by $\pi sr$.

The implementation of the Perez algorithm is contained in the class TSunposIrradiance.

The different irradiance/radiance components are provided by the function getPerezIrradiance:

```javascript
function TSunposIrradiance.getPerezIrradiance(const dt: Double; const DNI_, DHI_,
    GHI_: Double; const locAlbedo_, locPressHPa_, locLat_, locLon_, locAlt_: Double;
    const refPlaneV: Tvec3d; const refPlaneCosBeta, refPlaneSinBeta: Double;
    const UsePerezInterpolation: boolean; sunOpts: TSG2optParams;
    var irrDat: TIrradianceData): Boolean;

var
    y_, m_, d_: Word;
    sp: TsolarPos;
    F1, F2, b: Double;

begin
    DecodeDate(dt, y_, m_, d_);
    sp:=calcSunPos(JulianD(y_, m_, d_, Frac(dt), +1), loclat_, loclon_, localt_, sunOpts);
    with irrDat do
        Begin
            DHI:=DHI_;
            GHI:=GHI_;
            DNI:=DNI_;
            Az:=sp.alpha;
            El:=sp.gamma;
            ThetaZ:=pi/2-sp.gamma;
            V:=VRotOrigin(Vect(0,-1,0), XAxis, -El); // point vector to solar pos
            V:=VRotOrigin irrDat.V, ZAxis, -Az+pi);
            refPlanecosDNI:=VDot(V, refPlaneV); // equals coin of angle of incidence
            if (GHI<=0) or (DHI<=0) then // no irradiance measured, return 0 for components
                Begin
                    circumI:=0; diffSkyI:=0; diffHorizonI:=0; diffGndI:=0;
                    refPlaneTot:=0; refPlaneSky:=0; refPlaneGnd:=0; RefPlaneDNI:=0;
                    Exit;
                end;
            if UsePerezInterpolation=false then getPerezF1F2(dt, locPressHpa_, DHI, DNI, thetaZ, F1, F2, b)
                else getPerezF1F2interpol(dt, locPressHpa_, DHI, DNI, thetaZ, F1, F2, b);

            // Perez components for SIOP evaluation
            circumI:=DNI + F1*(DHI/b); // sum of DNI and Perez circumsolar forms directed component
            diffSkyI:=DHI*(1-F1); // diffuse Sky irr. ex circumsolar (result = radiance*Pi)
            diffHorizonI:=DHI*F2*(pi/(0.35564*0.63530)); // 6.5° horizon average (radiance*Pi)
            diffGndI:=GHI*locAlbedo_; // diffuse ground irr. ex circumsolar (result = radiance*Pi)

            // additionally components for on-plane irradiance, for reference
            if refPlanecosDNI<=0 then refPlaneDNI:=0
                else refPlaneDNI:=DNI*refPlanecosDNI + DHI*F1*(refPlanecosDNI/b);
            refPlaneSky:=DHI*(1-F1)*(1+refPlaneCosBeta)/2 + F2*refPlaneSinBeta;
            refPlaneGnd:=GHI*locAlbedo_*(1-refPlaneCosBeta)/2;
            refPlaneTot:=refPlaneDNI + refPlaneSky + refPlaneGnd;
        end;
end;
```
The four irradiance components `circumI`, `diffSkyI`, `diffHorizonI`, `diffGndI` are calculated as described in the previous section. However, the radiance values of the diffuse components (sky, ground, horizon) are already multiplied with a factor of $\pi sr$. Consequently, they are already in units of irradiance and can directly be multiplied with the dimensionless solid angles, as specified on the previous page.

The Perez parameters F1 and F2, which determine the circumsolar and horizon irradiance, can optionally be calculated based on the standard algorithm or based on the interpolated version suggested here. Only the function `getPerezF1F2interpol` is presented here:

```cpp
procedure TSunposIrradiance.getPerezF1F2interpol(const tim, locPressHPa: Double; const DHI, DNI, thetaz: Double; var F1, F2, b: Double);

Const
  f11 : array[1..8] of Double = (-0.008, 0.13, 0.33, 0.568, 0.873, 1.132, 1.06, 0.678);
  f12 : array[1..8] of Double = (0.588, 0.683, 0.487, 0.187, -0.392, -1.237, -1.6, -0.327);
  f13 : array[1..8] of Double = (-0.062, -0.151, -0.221, -0.295, -0.362, -0.412, -0.359, -0.25);
  f21 : array[1..8] of Double = (-0.06, -0.019, 0.055, 0.109, 0.226, 0.288, 0.264, 0.156);
  f22 : array[1..8] of Double = (0.072, 0.066, -0.064, -0.152, -0.462, -0.823, -1.127, -1.377);
  f23 : array[1..8] of Double = (-0.022, -0.029, -0.026, -0.014, 0.001, 0.056, 0.131, 0.251);

Const epsLoLimit : array[1..8] of Double =
  (1.0000, 1.0650, 1.2300, 1.5000, 1.9500, 2.8000, 4.5000, 6.2000);
Const epsCenter : array[1..8] of Double = (1.0325, 1.1475, 1.365, 1.725, 2.375, 3.65, 5.35, 7.05
  (=6.2+0.85 equal slope of previous two intervals, this is somewhat random));
const cos85deg = 0.08715574;

var eps : Double;
  kappath3 : Double;
  costhetaz : Double;
  delta : Double;
  c : Integer;
  c1, c2 : Integer;
  wt1, wt2 : Double;

begin
  kappath3:=1.041*thetaz*thetaz*thetaz;
  eps:=(DHI+DNI)/(DHI+1E-8) + kappath3)/(1+kappath3);
  costhetaz:=cos(thetaz);
  b:=max(cos85deg, costhetaz);
  delta:=(DHI*airMassKastenSandia(costhetaz, thetaz, locPressHPa))/getSeasonalExaSolarConst(tim);
  c:=1; while (c<8)and(eps>epsLoLimit[c+1]) do Inc(C);

  if (eps<epsCenter[c])or(c=8) then
    Begin // no interpolation for upper or lower end
      F1:=max(0, f11[c] + f12[c]*delta + thetaz*f13[c]);
      F2:= f21[c] + f22[c]*delta + thetaz*f23[c];
    end else
    Begin
      if eps<epsCenter[c] then Begin c1:=c-1; c2:=c; End else
        Begin c1:=c; c2:=c+1; End;
      wt2:=(eps-epsCenter[c1])/(epsCenter[c2]-epsCenter[c1]);
      wt1:=1-wt2;
    F1:=max(0, f11[c1] + f12[c1]*delta + thetaz*f13[c1])*wt1
        + max(0, f11[c2] + f12[c2]*delta + thetaz*f13[c2])*wt2;
    F2:=( f21[c1] + f22[c1]*delta + thetaz*f23[c1])*wt1
        + ( f21[c2] + f22[c2]*delta + thetaz*f23[c2])*wt2;
    end;
  end;
```
As depicted in Figure 90, a piecewise linear interpolation is applied in between the centres of the clearness $\varepsilon$ intervals defined in the standard algorithm (Perez et al., 1990; Sandia National Laboratories, 2018b). The standard Perez model's bin number representing integer values ranging from 1 to 8 is interpreted as a real number. Linear interpolation is used to replace the discrete eight values for the $f_{xx}$ parameters provided by Perez et al. with continuous values. The determination of the clearness index $\varepsilon$ function `getPerezF1F2interpol` relies on the values of the air-mass and extraterrestrial radiation. In order to determine the extraterrestrial

radiation, the commonly used Spencer approximation (Spencer, 1971; Sandia National Laboratories, 2018a) was implemented. For performance reasons, the constants are precomputed on a daily base. The commonly used model of Kasten (1964; Kasten et al., 1989) was implemented to determine the air-mass. Both implementations are provided below:

<figure><img src="image_95.png" /><figcaption>Figure 90 Perez clearness interpolation</figcaption></figure>

```javascript
function TSunposIrradiance.airMassKastenSandia(const costhetaz_, thetaz, locPressHpa: Double): Double;

const amLimitD = 89.9;
var thetazD: Double;
    cosThetaz: Double;
begin
    thetazD:=thetaz/pi*180;
    if thetazD>amLimitD then
        Begin thetaZD:=amLimitD;costhetaz:=0;End
    else cosThetaz:=costhetaz_;

    result := 1/( costhetaz+0.50572*power(96.07995-thetazD,-1.6364) );
    result:=result*locPressHpa/1013.25;
end;
```

```pascal
procedure TSunposIrradiance.preCalcSeasonalSolarConst;
var b: Double;
    rfctSqr: Double;
    c: integer;
begin
    For c:=1 to 366 do
    Begin
        b:=2*pi*(c/365);
        if b>2*pi then b:=2*pi;
        rfctSqr:=1.00011 + 0.034221*cos(b) + 0.00128*sin(b) + 0.000719*cos(2*b)
            + 0.000077 * sin(2*b); // this equals (R/Ravg)^2 R...sun-earth distance
        seasonalSolarConst[c]:=1367*rfctSqr; // save values in precalc. array
    end;
end;
```

The presented algorithms provide all the required elements for evaluating the SIOP for any given irradiance profile. The actual evaluation, representing the implementation of equation (125), is now merely a matter of a few multiplications. The full code required to finally evaluate the irradiance considering the diffuse and direct components can simply be written as:

```pascal
siopI.Dir:=siop.evalSIOPWorld(sun.V,frontSideHit)*circumI;
siopI.Sky:=siop.ssiop.solidAngleDivPIupperHemi*siop.ssiop.meanUpperHemi*diffSkyI
    + siop.ssiop.solidAngleDivPIhorizBand*siop.ssiop.meanHorizBand*diffHorizI;
siopI.Re:=siop.ssiop.solidAngleDivPIlowerHemi*siop.ssiop.meanLowerHemi*GHI*locAlbedo;
siopI.Tot:=siopI.Dir+siopI.Sky+siopI.Re;
```
## 6.2.6. Testing and validation

While many consistency and plausibility tests were carried out during the implementation of the relevant modules, an explicit two-level test for evaluating the Perez module was performed. In the first step, the validity of the Perez module, providing the individual radiance and irradiance components, was tested by calculating the time series for two different days and orientations. Afterwards, the results of the SIOP-based evaluation for horizontal and vertical orientation were computed. All results were compared against irradiance values calculated with the building performance simulation software IDA ICE (Equa S.A., 2021). The program features the Perez irradiance model and allows separate export of the irradiance components (diffuse sky, ground and direct) on different vertical planes and the horizontal plane ("roof").

Measured high-resolution (1 second time steps) irradiance values of the local weather station for two days were used as input data. The first day (14.6.2021) had a clear sky, while the second (17.7.2021) was mostly overcast. The location of the weather station is the rooftop of the AEE INTEC laboratory in Gleisdorf (latitude: 47.1097°, longitude: 15.7096°, altitude: 369 m). The albedo value was set to 10% in both calculations.

### 6.2.6.1. Validation of Perez module - comparison of diffuse components

First, the correct implementation of the Perez algorithm was tested. For this purpose, the diffuse sky components, and additionally the ground component, were compared. According to the Perez model, the diffuse sky irradiance depends on several parameters in a complex way, including direct normal irradiance and day-of-year. Consequently, comparing computed diffuse irradiance time series for different orientations represents a comprehensive and robust validation of the algorithms involved.

Since the RadiCal implementation contains an additional, optional interpolation feature, two calculations were tested. The results of the validation are shown in Figure 91.

<table><caption>Table 12 Error statistics, comparison of Perez diffuse sky irradiance RadiCal vs. IDA ICE</caption><thead><tr><td>case</td><td>R<sup>2</sup></td><td>mean error (bias) [W/m<sup>2</sup>]</td><td>mean abs. error [W/m<sup>2</sup>]</td><td>max abs. error [W/m<sup>2</sup>]</td></tr></thead><tbody><tr><td>south, 14.6.2021</td><td>0.99987</td><td>0.06</td><td>0.18</td><td>3.74</td></tr><tr><td>south, 14.6.2021 with interpolation</td><td>0.99955</td><td>-0.03</td><td>0.47</td><td>3.28</td></tr><tr><td>horiz., 17.7.2021</td><td>0.99863</td><td>0.16</td><td>2.29</td><td>30.6</td></tr><tr><td>horiz., 17.7.2021 with interpolation</td><td>0.99864</td><td>0.15</td><td>2.29</td><td>30.6</td></tr></tbody></table>

The comparison of the diffuse sky irradiance time series shows an excellent agreement between the RadiCal and IDA ICE calculations. The coefficients of determinations and errors are listed in Table 12. As can be seen, the calculation results are mostly identical. The singular, relatively high maximum error for the horizontal is likely caused by the different handling of time frames in connection with the strong short-term fluctuations. The irradiance data measured in one-second intervals have to be averaged to minute intervals. Since obviously different averaging methods or different reference times are applied here, strong short-term irradiance fluctuations lead to a deviation in the results. The small values for the mean error indicate this, as the fluctuations average out when mean values are calculated.
<figure><img src="image_96.png" /><figcaption>Figure 91 Comparison of Perez diffuse sky irradiance
top row: south orient., 14.6.2021, left: std. method, right: interpolation
bottom row: horizontal orient., 17.7.2021, left: std. method, right: interpolation</figcaption></figure>

It is also apparent that the interpolation feature of *RadiCal* has no significant impact on the comparison regarding average errors or correlation. Nevertheless, the feature leads to noticeable smoothing. Albeit the quantitative effect is small, this represents an appropriate improvement of the calculation, as it leads to a more realistic dynamic behaviour. It removes the transient, step-like changes, benefiting dynamic thermal simulations. Figure 92 shows a zoomed detail of the comparison. It can be seen that the standard Perez model (grey) leads to step-like changes in the diffuse irradiance, whereas the interpolation feature leads to an unbiased, smoothed curve.

<figure><img src="image_97.png" /><figcaption>Figure 92 Comparison of Perez diffuse sky irradiance.</figcaption></figure>
### 6.2.6.2. Validation of SIOP evaluation vs. IDA ICE calculation

In this next step, the entire SIOP evaluation process is validated. Again, the Lambert-SIOP (see Figure 83) representing a planar absorber is used. Also, the same days as for the previous steps are used for the evaluation. The tests are performed for three vertical orientations (south, west, north) and the horizontal plane.

The results of the evaluation are presented in Figure 93, Figure 94 and Table 13.

<figure><img src="image_98.png" /><figcaption>Figure 93 Comparison SIOP evaluation vs. IDA ICE (Perez model)
clear sky; 14.06.2021</figcaption></figure>
<figure><img src="image_99.png" /><figcaption>Figure 94 Comparison SIOP evaluation vs. IDA ICE (Perez model)
mostly overcast; 17.07.2021</figcaption></figure>

<table><caption>Table 13 Error statistics, comparison SIOP evaluation for plane vs. IDA ICE for irradiance on plane</caption><thead><tr><td>case</td><td>R<sup>2</sup></td><td>mean error (bias) [W/m<sup>2</sup>]</td><td>mean abs. error [W/m<sup>2</sup>]</td><td>max abs. error [W/m<sup>2</sup>]</td></tr></thead><tbody><tr><td>horizontal, 14.6.2021</td><td>0.99997</td><td>2.4</td><td>2.5</td><td>13.3</td></tr><tr><td>south, 14.6.2021</td><td>0.99993</td><td>2.7</td><td>2.8</td><td>9.1</td></tr><tr><td>north, 14.6.2021</td><td>0.99370</td><td>-1.3</td><td>3.0</td><td>48.4</td></tr><tr><td>west, 14.6.2021</td><td>0.99942</td><td>0.0</td><td>1.6</td><td>56.4</td></tr><tr><td>horizontal, 17.7.2021</td><td>0.99084</td><td>-1.8</td><td>7.6</td><td>138.1</td></tr><tr><td>south, 17.7.2021</td><td>0.99117</td><td>-0.3</td><td>2.9</td><td>72.9</td></tr><tr><td>north, 17.7.2021</td><td>0.99862</td><td>-0.8</td><td>1.1</td><td>23.4</td></tr><tr><td>west, 17.7.2021</td><td>0.98529</td><td>-1.1</td><td>7.3</td><td>203.8</td></tr></tbody></table>

It can be seen in the charts of Figure 93 and particularly Figure 94 that there are noise-like deviations at times with strongly fluctuating irradiance values. As mentioned in the previous section, these deviations are the results of different averaging processes or different time-frames used to generate minute interval data of the original measurement data with one-second time steps. This is confirmed by the mean errors (Table 13), as they show no significant bias. Apart from the fluctuations, the calculated irradiance values show an almost perfect match (see Figure 93 and Figure 94 - dashed red curve overlapping grey curve). The only significant systematic deviation can be found in the case north, 14.6.2021 (see Figure 93, bottom left). A systematic deviation of approximately 5 W/m²
occurs in the morning and -9 W/m² in the evening. Surprisingly, the direct normal irradiance components differ at these times. The results of the *RadiCal* implementation were cross-checked by comparison with the standard Perez model, where no error was detected. The reason for the deviating calculation of IDA ICE remains unclear. As the solar positions calculated were found to be almost identical, only a minor tilt of the northern plane could potentially explain the occurring deviation.

### 6.2.6.3. Validation of *Lambert-SIOP* evaluation vs. standard on-plane irradiance model

In a further validation step, the *SIOP* evaluation results for the *Lambert-SIOP* were compared with on-plane irradiance calculations using the Perez method. For the latter, the *RadiCal* Perez module was used, as this module has already been validated against the IDA ICE calculations in the first step (see above). As both calculations can now be performed using the same time frames, no sampling-related noise is affecting the comparison. The results of both calculations were found to be virtually identical, showing average deviations in the range of a few thousandths Watt per square metre. Consequently, the results are not presented in graphical form. The comparison of the *Lambert-SIOP* evaluation vs. the standard Perez algorithm seems trivial. Still, it is a robust and comprehensive proof of the implementation's validity, as the two calculation paths differ significantly. While the standard calculation is based on trigonometric functions of the tilt angle and solar angle, the *SIOP* evaluation is based on numerical integrations and vector calculus. The perfect match of the resulting irradiance on planes with different orientations confirms, in particular, the validity of these evaluation steps:

* generation of *SIOP*

* evaluation of *SIOP* (for singular direction)

* integration on *SIOP* over segments (for average values)

* calculation of solid angles of segments (sky, horizon, ground)

* determination of radiance values for segments

* two-axis rotation into the world coordinate system
# 7. A full spectral polarisation Monte-Carlo Raytracer

## 7.1. Tracing a ray - overview

In the field of computer graphics, many different implementations of Monte-Carlo (MC) raytracing can be found. A good overview is provided by Jansen *et al.* (2003). The methods apply different advanced methods of MC sampling known as, e.g. Metropolis sampling, importance sampling or Russian roulette. All methods aim at providing realistic renderings of objects with the highest possible efficiency, i.e. the shortest possible calculation time. In order to achieve this, the methods cast rays (light particles) from the eye-point towards the viewed area. Since the light is traced against its natural direction, these methods are referred to as backward raytracing, while tracing the light in its natural direction is referred to as forward raytracing (unfortunately, some authors with computer graphics background use the terms with interchanged meanings).

Backward raytracing is very efficient on the first hand, as it guarantees that every cast ray hits an object in the visible area first and therefore provides valuable information for the rendering. The efficiency of the methods is then determined after this first hit when the light incident to this first collision point has to be sampled efficiently. The light incident to this point will generally show high spatial fluctuations, as usually distinct, angularly limited light sources, or reflections of these, are present. The various MC sampling strategies aim to solve this integration of the incoming light as efficiently as possible. Special requirements must be met in this sampling process, as the visual appearance of the object, and not the energy of the light, is the relevant measure.

Reversing directions and using forward raytracing for image rendering is possible in theory but not from a practical point of view, as it is simply too inefficient. Forward raytracing, e.g. of an outdoor environment, would mean that every ray of light initiates from the sun's direction or some from the sky dome if scattering is considered. The rays are then scattered on objects, and only a very small fraction of them would finally reach the virtual camera lens to provide imaging information. However, most rays would either be absorbed or scattered in other directions. The overwhelming majority of cast rays would therefore be useless for the rendering process.

However, a similar reasoning applies in reverse in the application presented here. The focus of this application is determining the paths that the energy (or power) carried by the incident light particles are taking. In this context, any light particle contributes equally relevant information, regardless of whether it is absorbed or transmitted through the object. Hence, an unbiased sampling strategy fulfils the requirements of this energetic analysis best. In the context of rendering and backward raytracing, this unbiased strategy would cause the appearance of some surfaces to be determined by very few samples only, leading to noise-like colour or brightness fluctuations in the image. In the energetic context, however, the number of samples absorbed by a surface is always related to the density of the incidence light. Therefore, the relative accuracy of the results is determined by the total of light particles cast and can, therefore, even be high for a surface absorbing no particles.

Following this line of reasoning, an important design choice is that no recursive path splitting is applied during the raytracing process. This means that whenever light is scattered into different directions or simultaneously reflected, absorbed and transmitted, always one path is selected and followed further on, maintaining the full energy of the light particle. Since the probability of these paths is determined based on the underlying probabilities, the law of large numbers ensures that the
results converge to the expected value. Splitting the path into N sub-paths would lead to a recursive division, resulting in a large number of branched paths. On average, all of these paths would be computationally as expensive as the single path; however, most would carry little energy and, therefore, contribute little relevant information. Avoiding path splitting is, therefore, significantly more efficient for the intended application.

## 7.2. Rays and light particles

In a not entirely accurate – analogy to quantum mechanics, the term “light particle” is used alongside the term “ray” below. It is helpful to have the concept of light particles in mind when scattering events are considered in raytracing, while the concept of rays is useful to describe linear propagation. The term *light particle* does, however, not refer to a specific physical entity but to the sampling entity used in the implementation, as the ray is sequentially traced in the form of a sample from its generation over scattering events to absorption. The term *photon* is avoided, as it relates to the quantum dynamic nature of light, which can, of course, not be considered here in full. The energies of photons further depend on their wavelengths, whereas the introduced *light particles* are defined as being equienergetic, i.e. each sample represents the same constant flux density fraction, regardless of its wavelength. Beyond that, a Stokes vector is assigned to each light particle to describe its polarisation state. As discussed above (sections 4.2.4 and 4.8.1), the Stokes vector describes the polarisation state of a multitude of photons, or wave trains, in a statistical nature. Hence, the light particles introduced here can be best described as a fraction of a parallel light beam consisting of many photons with the same wavelength and a specified constant radiant flux density.

In the implementation, the following type definitions are used to describe the ray/light particle:

```csharp
Type TCollision = record
    cp : TVect3d; // collision point
    cn : TVect3d; // surface normal at collision point
    cxv, cyv : TVect3d; // local x,y vector at collision
    cu, cv : Double; // barycentric collision coordinates
    sprod : Double; // dot product at collision
    oNo : UInt32; // object index of hit object
    fNo : UInt32; // face index of hit object
    mNo : UInt32; // material index of hit object
    enteredTransparentMediumAlpha : Double; // abs. coeff of mat. hit
    lsiR : TLSIevent; // type of scattering event
    scatterCnt : UInt32; // count of scattering events suffered
End;
```

```csharp
Type TStokesRay = record
    rp, rdir : TVect3d; // starting point and direction
    pState : TStokesVect; // polarisation state
    spol : TVect3d; // direction of S-Polaritation unit vect,
                    perpend. to rdir (P-pol not explicit)
    wl : Double; // wavelength in m
    curColl : TCollision; // current collision
    lastColl : TCollision; // recent collision
    collCn : Integer; // count of scatterings suffered
    rinvDir : TVect3d; // required for fast box-test
    xDirSign, yDirSign, zDirSign : Integer; // required for fast box-test
    0:positive, 1:negative
End;
```

In the implementation, a light particle or ray sample is defined by the type `TStokesRay` in the global module `RC_Vects`. It includes the current position and direction (`rp`, `rdir`), the polarisation state of Stokes vector (`pState`) and reference axis (`spol`) and the wavelength (`wl`).

It carries information on the current (`curColl`) and most recent scattering event (`lastColl`), as well as an ordinal counter that tracks the collisions (`collCn`). It further contains four variables that are used to increase the performance of the collision detection algorithm.
The two encapsulated `Tcollision` instances contain information regarding the location, orientation and character of the two most recent object collisions representing scattering events (reflection or transmission).

## 7.3. 3D model and geometry modules

### 7.3.1. OBJ file import

In order to maximize the potential application of the software, an import module allowing the import of standard CAD files was implemented. There are various common file formats used to store and deploy three-dimensional CAD models. Because of its simplicity and widespread use as an exchange format, the OBJ format (The Library of Congress (USA), no date; Bourke, 2022) has been chosen for implementation.

The format is clear-text-based, so it can be viewed with any text file editor. The following sample shows the data of a simple triangulated cube with a dimension of one located at the origin (see Figure 95).

<figure><img src="image_100.png" /><figcaption>Figure 95 Triangulated cube with normal vectors</figcaption></figure>

<table><tr><td>o Cube</td><td>// object name</td></tr><tr><td>mtllib testcube2.mtl</td><td>// reference to material definition file</td></tr><tr><td>v 1.000000 1.000000 -1.000000</td><td>// list of all object vertices</td></tr><tr><td>v 1.000000 0.000000 -1.000000</td><td></td></tr><tr><td>v 1.000000 1.000000 0.000000</td><td></td></tr><tr><td>v 1.000000 0.000000 0.000000</td><td></td></tr><tr><td>v 0.000000 1.000000 -1.000000</td><td></td></tr><tr><td>v 0.000000 0.000000 -1.000000</td><td></td></tr><tr><td>v 0.000000 1.000000 0.000000</td><td></td></tr><tr><td>v 0.000000 0.000000 0.000000</td><td></td></tr><tr><td>vn -0.5773 0.5773 -0.5773</td><td>// list of normal vectors for vertices</td></tr><tr><td>vn 0.5773 0.5773 0.5773</td><td></td></tr><tr><td>vn 0.5773 0.5773 -0.5773</td><td></td></tr><tr><td>vn -0.5773 -0.5773 0.5773</td><td></td></tr><tr><td>vn 0.5773 -0.5773 0.5773</td><td></td></tr><tr><td>vn -0.5773 -0.5773 -0.5773</td><td></td></tr><tr><td>vn 0.5773 -0.5773 -0.5773</td><td></td></tr><tr><td>usemtl Material</td><td>// material identifier for following faces</td></tr><tr><td>s 1</td><td>// smoothing activated for following faces</td></tr><tr><td>f 5//1 3//2 1//3</td><td>// face definitions using</td></tr><tr><td>f 3//2 8//4 4//5</td><td>// Vertex-IDs and Normal-vector IDs</td></tr><tr><td>f 7//6 6//7 8//4</td><td></td></tr><tr><td>f 2//8 8//4 6//7</td><td></td></tr><tr><td>f 1//3 4//5 2//8</td><td></td></tr><tr><td>f 5//1 2//8 6//7</td><td></td></tr><tr><td>f 5//1 7//6 3//2</td><td></td></tr><tr><td>f 3//2 7//6 8//4</td><td></td></tr><tr><td>f 7//6 5//1 6//7</td><td></td></tr><tr><td>f 2//8 4//5 8//4</td><td></td></tr><tr><td>f 1//3 3//2 4//5</td><td></td></tr><tr><td>f 5//1 1//3 2//8</td><td></td></tr></table>

Firstly a list of vertices ("v") and normal vectors ("vn") is contained. In the bottom section, the indices of these elements are referenced ("f") to define 12 triangular surfaces. While the format allows additional features, see e.g. (Bourke, 2022), the data contained in this simple example covers all types that are required for importing geometry to the *RadiCal* method.

Beyond the features contained in the above example, only the functionality to import texture coordinates is additionally implemented. The texture feature can potentially be used to define the directional system for non-isotropic surface roughness. However, this feature is currently not fully implemented in the raytracer.
### 7.3.1.1. Implementation

The `OBJ` import is contained in class `TImportOBJ` of the module `rc_OBJECT`.

```python
Type TSubDivLimitsEntry = record nm : String; lim : Double; End;
Type TSubDivLimitsArray = array of TSubDivLimitsEntry;

Type
{
    TImportOBJ
}

TImportOBJ = class(TObject)
    private
        subDivMaxSize : Double;
        subDivMaxSizes : TSubDivLimitsArray;
    public
        constructor Create;
    function ReadOBJ(const fn: String;
        var geo: TRCGeometry;
        var mlib: TRCMatLib;
        var errList: TStringList;
        const SmoothNormalMinAngleDevDeg, SmoothNormalMaxAngleDevDeg: Double;
        const subDivideLimit: TSubDivLimitsArray): Integer;
        destructor Destroy; override;
end;
```

The class contains solely the function `ReadOBJ` which performs the import. All required additional routines are contained within this function. The function returns the imported geometry (`geo`), extracted material list (`mlib`) and a list with errors or warnings (`errList`) that can potentially occur during the import.

Apart from the filename (`fn`) the function has two input parameters that are used to control the function of the smooth shading feature (see 7.6). Additionally, an optional array can be passed to the function triggering an automatic subdivision of the geometry during the import. This array provides the relevant material identifiers (as strings) and the subdivision limit that should be applied during the import for the specified material layer.

### 7.3.2. Internal geometry definitions

Internally, the imported geometry is stored in a hierarchical structure:

The master node (`TRCGeometry`) contains a list of objects (`TRCGeoObject`), and each object contains a list of faces (`TRCface`). Along with the objects, the boundary boxes (relevant to collision detection) are included at each level.

```python
Type TRCFace = record // face-class, either triangle or quad (parallelogram)
    v0 : TVect3d; // origin vertex
    edge1, edge2 : TVect3d; // vectors for edges
    snv : TVect3d; // surface normal vector = z -axis (for flat surfaces)
    lxv, lyv : TVect3d; // X,Y axis of aligned local system
    is4P : boolean; // true:quad, false: triangle
    frontID : Integer; // id of frontside material
    backID : Integer; // ID of backside material
    sn0, sn1, sn2 : TVect3d; // optional surface normal vector
    // for smooth shading (triangles only)
end;
```

```python
Type TRCGeoObject = record // object class
    objName : String; // object identifier
    visible : Boolean; // visibilty flag
    faceCnt : Integer; // face count
    face : array of TRCFace; // face list
    faceBB : TBoundBoxArray; // face bdry. boxes
end;
```

```python
Type TRCGeometry = record // geometry class
    objCnt : Integer; // object count
    obj : array of TRCGeoObject; // object list
    ObjBB : TBoundBoxArray; // obj bdry. boxes
end;
```
As shown in the class definitions above, the 3D data contained in the OBJ file is stored in a hierar-
chical structure that additionally contains pre-computed data to increase the performance of the
raytracing algorithm.

The three vertices provided by the
OBJ face definition are converted into
a face origin (**v0**) and two vectors
(**edge1**, **edge2**). Based on the edge
vectors, a local reference system is pre-
computed for every face (**lxv**, **lyv**,
**snv**). The x-axis of the reference sys-
tem is aligned to the **edge1** vector,
whereas the z-axis is aligned to the
surface normal vector of the face. If
the imported model contains normal
vectors defined on the vertices (**s0**,
**s1**, **s2**) the smooth shading function
can be used to derive surface smooth
normal vectors throughout the face
(see 7.6). If no vertex normals are contained in the file, a single constant surface normal vector (**snv**)
is generated by applying the cross-product of **edge1** and **edge2**. The direction of the surface normal
vector is used to distinguish the front and backside of the surface. The definition follows the conven-
tion of computer graphics: a counter-clockwise sequence of the vertices defines the front side of the
surface, while the clockwise order defines its back side. Based on the ordinal variables **frontID** and
**backID**, each side can be assigned a different surface material.

<figure><img src="image_101.png" /><figcaption>Figure 96 Vector definitions contained in the TRCface class</figcaption></figure>
## 7.4. Collision detection

The objective of the collision detection algorithm is to determine the intersection of a ray with the geometric model. From a geometrical point of view, this means that the closest intersection point of the line representing the ray and any face of the model has to be determined (see Figure 97). In order to solve this efficiently, two components are required: first, an algorithm to calculate the intersection point between any face and the ray, i.e. a line, is required. Second, a partitioning algorithm is needed to reduce the number of necessary intersection tests.

### 7.4.1. Line triangle intersection

Only the line-triangle intersection algorithm is applied in the implementation, as already triangulated geometry is used. Models that do contain quad faces

(consisting of four points in a plane) can trivially be split into two triangles for the test. There are many different approaches for calculating the intersection point. In the current implementation, the Müller-Trumbore algorithm (Möller and Trumbore, 1997) is applied. The test, on the one hand, has a high algorithmic efficiency and, on the other hand, can already provide essential information regarding the location of the intersection point relative to the ray origin point and in the coordinates of the face hit.

As can be seen in the implementation below, the algorithm is quite short if vector notation is used.

<figure><img src="image_102.png" /><figcaption>Figure 97 Collision detection to determine intersection point closest to ray origin ($\vec{P}_0$) in direction $\vec{d}$</figcaption></figure>

```javascript
function TFaceTree.RayHitT(const rayP0, rayDir: TVect3d; const Tri: TRCface; var ipoint: TVect3d; var u,v,t: Double): Boolean;
    Const Eps = 0.0000001; // MÖLLER-TROMBORE test
    var h,s,q: TVect3d;
        a,f,u,v: Double;
    begin
        h:=VCrossProdraydir,tri.edge2);
        a:=Vdot(tri.edge1,h);
        if (a>-eps)and(a<eps) then Begin Result:=false;Exit;End; // ray is parallel->leave
        f:=1/a; // scaling parameter
        s:=rayp0-tri.v0;
        u:=f*Vdot(s,h);
        if (u<0)or(u>1) then Begin Result:=false;Exit;End; // hit outside tri.area->leave
        q:=VCrossProd(s,tri.edge1);
        v:=f*Vdot(raydir,q);
        if (v<0)or((u+v)>1) then Begin Result:=false;Exit;End; // hit outside tri.area->leave
        t:=f*Vdot(tri.edge2,q);
        if t>eps then // distance origin->intersection point is positive
            Begin // hit provide barycentric coordinates and intersection point
                u_:=u;
                v_:=v;
                ipoint:=rayp0+t*raydir;
                Result:=true;
            End else Result:=false;
        end;
```

The parameters ray origin point (rayP0) and direction (rayDir) as well as the face (Tri) are passed as input to the function. The function is of boolean type and returns false if the triangle is missed. Otherwise, it returns true and passes the intersection point in world coordinates (ipoint) as well as
in barycentric coordinates of the triangle forming the face (u_, v_). Beyond that, it outputs the parameter t that contains the distance between the ray origin and the intersection point. This information is crucial to determine the closest collision point along the ray's path and to avoid unnecessary calculations for further distances. As can be seen in the code, a part of the efficiency of the algorithms results from the integrated early termination tests. At various points, the execution is aborted if the current calculation results no longer allow an intersection.

### 7.4.2. Spatial partitioning – search tree

The triangle intersection test could simply be applied to all faces of the model to finally determine the closest intersection point based on the returned distance parameters t. However, this is a very inefficient approach that would require vast computational time. In a much more efficient approach, a spatial partition is applied to the model first. The partition is used to create a search tree for the faces that allows identification of potential target faces more easily - or conversely expressed, the exclusion of non-relevant regions from the test.

The partitioning can be performed based on various algorithms. The most common ones are called *uniform subdivision* (Cleary and Wyvill, 1988), *kd-tree* (Ooi, 1987) and *octree* (Meagher, 1980). The latter two show significantly higher performance. While the actual performance

depends on the specifics of the geometric model, a slightly higher performance can be expected of the kd-tree for many models (Szirmay-Kalos et al., 2002). However, since the octree method is less complex to implement and the performance difference is limited, the octree method was implemented here. It can easily be replaced or extended with other methods later on.

The implemented octree method is based on the initial ideas by Meagher (1980) and many further implementations. The algorithm was, however, implemented from scratch for the *RadiCal* method, applying today's available modern programming features such as object-oriented programming and dynamic memory allocation. The method begins with a bounding box containing the entire model. This root box is then recursively subdivided into eight equally sized sub-boxes, and the number of faces contained in each sub-box is analysed. The subdivision is carried on recursively if the face count lies above a certain threshold. Otherwise, the subdivision is aborted, and a list with indices of the contained faces is created for the last element. As commonly applied for search trees, all elements containing face lists are referred to as “leafs”, while all elements that refer to child boxes are referred to as “nodes” (see Figure 109). Both elements, contain the coordinates of the associated bounding box.

<figure><img src="image_103.png" /><figcaption>Figure 98 Octree partition of imported CAD model to increase efficiency of collision detection</figcaption></figure>
<figure><img src="image_104.png" /><figcaption>Figure 99 Octree hierarchy</figcaption></figure>

There are mainly two advantages resulting from this partitioning. First, the nested hierarchy effi-
ciently excludes large parts of the model if the bounding boxes of higher-level nodes do not intersect
with the ray. Second, the required box-line intersection tests, which are always performed before the
final face intersection tests, can be performed very fast, as the rectangular bounding boxes are aligned
to the axis of the coordinate system. More considerations regarding the algorithm can be found in
Meagher’s original publication (1980).

A customized implementation has been written for use in the raytracer. It is contained in the class
TFaceTree of the module rc_facetree.

Type TOFListEntry = record
    objNo : Integer;
    faceNo : Integer;
End;
Type TOFList = array of TOFListEntry;
Type TFaceTreeNode = record
    bb : Txyzboundbox;
    child : array [0..7] of Integer;
end;
Type TFaceTreeLeaf = record
    bb : Txyzboundbox;
    oflCnt : Integer;
    ofl : TOFList;
end;
type TFaceTree = class(SharpObject) { TFaceTree }
private
    geo : TVGeo;
    maxLeafSize : Integer;
    minSqrDim : Double;
    nodeCnt : Integer;
    node : array of TFaceTreeNode;
    leafCnt : Integer;
    leaf : array of TFaceTreeLeaf;
    function AddLeaf: Integer;
    function RayHitT(const rayP0, rayDir: TVect3d; const Tri: TRCface;
        var ipoint: TVect3d; var u_,v_,t: Double): Boolean;
    procedure SetupNode(const nno : Integer; const bb : txyzboundbox);
    function CreateList(const bb: txyzboundbox; var ofl: TOFList;
        const abortOnMax: Boolean): Integer;
    function AddNode : Integer;
    procedure truncNodeLeafSize;
public
    constructor Create;
    procedure createTree (var geo_ : TVGeo; const bb0: txyzboundbox;
        const maxLeafSize: Integer; const minDim: Double;
        const isoScale, minDimOfVolume: boolean);
    function treeStats : String;
    procedure save(const fn : String);
    function readFromFile (var geo_ : TVgeo; const fn: String): Boolean;
    procedure dumpOcttree(const fn : String);
    procedure dumpOcttree2Objects (const fn: String);
    function nextRayCollRecursive (var ray: TStokesRay): Boolean;
```cpp
destructor Destroy; override;
end;
```

There are essentially only two functions relevant for the application of the module: `createTree` and `nextRayCollRecursive`. The first one is used to establish the octree. For this purpose, the geometry (`geo_`) and root bounding box (`bb0`) are passed to the module. In order to control the tree generation, the maximum number of faces per leaf is provided (`maxLeafSize`). During the implementation, two additional features were implemented and found to be useful for increasing the performance on some geometries: `minDim_` allows defining a secondary recursion abortion criterion based on the volume of the bounding box. Finally, the `isoScale` parameters can be used to transform the root box into a cube. If this option is not chosen, the proportions of the initial bounding box will be maintained throughout the recursion. Since the tree generation can take up several minutes, the functions `saveToFile` and `readFromFile` can be used to store and reuse an already generated hierarchy.

Finally, the function `nextRayCollRecursive` is the key function that performs the collision tests in the raytracing process. In the exceptional case of not detecting any collisions, it returns `false`; otherwise, it updates the `currColl` variable in the `ray` parameter with specific information regarding the detected collision. As shown in section 7.3.1, the collision structure contains comprehensive information for the face hit. Most of these parameters are provided by `nextRayCollRecursive` already. This increases the performance of the raytracing, as calculations results from the triangle intersection can be used for this purpose. Specifically, the following parameters of `currColl` are provided:

```cpp
cp : TVect3d; // collision point
cn : TVect3d; // surface normal at collision point
cxv,cyv : TVect3d; // local x,y axis at collision point
cu,cv : Double; // barycentric (triangle) collision coordinates
sprod : Double; // dot product at collision (incidence angle)
oNo : UInt32; // object index of hit object
fNo : UInt32; // face index of hit object
mNo : UInt32; // material index of hit object
```
## 7.5. Periodic planes

In currently used approaches, both in the standard methods as well as in the BSDF-based methods, the surrounding objects of glazings, e.g. frames or rails for the shadings, are usually neglected. This can, on the one hand, lead to significant inaccuracies, but on the other hand, provides results of a more general nature, i.e. the results can later be used to model glazings of any size.

In order to be able to perform such calculations efficiently with the *RadiCal* method, a periodic-plane feature is implemented. Up to two sets of parallel planes can be included in the CAD drawing and defined as periodic planes. Each plane has a front and back side, defined by its surface normal vectors. Any rays intersecting these planes from the front side are relocated to the corresponding, opposite plane. All translations are performed parallel to the specified axis, i.e. two coordinates are held constant. Rays that reach the plane coming from the backside are not affected.

This concept allows for efficient modelling of periodic elements, as only a single element (or any multiple) has to be modelled, provided it is surrounded by periodic planes that correctly reflect the periodic lengths of the structure. E.g. for modelling Venetian blinds (see Figure 100), the vertically aligned planes (Z+, Z-) have to be placed precisely at a distance matching the slat distances, while any distance can be chosen for the horizontally aligned periodic planes (X+, X-).

<figure><img src="image_105.png" /><figcaption>Figure 100 Principle of ray relocation at periodic planes</figcaption></figure>

The concept of periodic planes can either be used to model “quasi-homogeneous” surfaces, e.g. glazings shaded by Venetian blinds, as shown in the example in Figure 100, or - more generally - to model structures that are composed of periodic structures, e.g. curtain wall façades or photovoltaic arrays.
## 7.6. Smooth shading

In order to be able to describe curved surfaces efficiently, a smooth shading feature is implemented in the *RadiCal* method. In the context of the *RadiCal* raytracer, the term smooth shading is somewhat misleading, as the surfaces are not actually shaded. However, the term mostly refers to the same principles when used in computer graphics. Technically, smooth shading is an interpolation method for surface normal vectors on curved surfaces. The idea dates back to an algorithm introduced

by Henri Gouraud (1971) and is fairly simple: it assumes that the surface normal vectors on a smooth surface change continuously but are explicitly defined at the vertices of the surface mesh only. Linear interpolation is applied to derive the surface normal vector at any point of the given normal vectors of the surrounding vertices. A threshold value usually determines if smooth shading is applied on a face. It is defined as a maximum angle (dot product of surface normal vectors) that determines if the face is part of a curved surface or represents an edge of the object. In the latter case, no smooth shading is applied.

The integration into the *RadiCal* raytracing algorithm is simple, as the implemented collision detection algorithm (see above) already provides barycentric coordinates for the collision point inside a triangular face. Hence, the linearly interpolated surface normal vector $\vec{n}_p$ at the collision point P (see Figure 101) can be found by:

$$ \begin{aligned} \vec{n}'_p &= (1 - u - v) \cdot \vec{n}_a + u \cdot \vec{n}_b + v \cdot \vec{n}_c \\ \vec{n}_p &= \frac{1}{|\vec{n}'_p|} \vec{n}'_p \end{aligned} \quad (129) $$

$\vec{n}_a, \vec{n}_b, \vec{n}_c$ ... surface normal vectors at vertices A, B, C
u, v ... barycentric coordinates in triangle ABC

The smooth shading feature is a simple but essential feature of the *RadiCal* method. It significantly improves the raytracer's performance (or accuracy), as considerably fewer vertices are required to accurately model the reflection of curved surfaces (e.g. the lamellae of shades).

<figure><img src="image_106.png" /><figcaption>Figure 101 Smooth shading / linear interpolation of surface normal vectors</figcaption></figure>

<figure><img src="image_107.png" /><figcaption>Figure 102 Smooth shading example – partially silver coated ETFE membrane cushion
(left: wireframe model; right: RadiCal rendering)</figcaption></figure>
## 7.7. Material-Interface class

Material interface objects are assigned to the surfaces of the imported 3D model. The assignment is based on the material identifier already present in the CAD file. However, the front side and the back side of the surface are generally assigned differently (see 7.3.2). The internal faces of opaque objects can be omitted, as no ray should impact there. This allows raising exception warnings if inner faces are still accessible, which would then indicate an error in the geometry.

The class definition of the `TRCMaterialInterface` is contained in the module `TRCMaterial`.

```python
Uses rc_Fresnel, rc_refractivityIndex, rc_MuMatMap, rc_subsurfaceRho;

Type TRCMaterialInterface = class(TObject) { TRCMaterialInterface }
private
public
    name : String;                 // name of mat-interface
    Substrate : TRIbaseclass;       // "frontside" medium RI
    Medium : TRIbaseClass;         // "backside" medium RI
    ThinFilmCnt : Integer;         // number of present thin-films
    ThinFilm : TEightLayerRIBaseClassList; // thin-film refractive indices
    ThinFilmD : TEightLayerThicknessList; // thin-film thicknesses (meter)
    Rcorrelated : Boolean;         // use correlated microfacet roughness
    Ralpha_x, Ralpha_y : Double;    // roughness alphas for microfacets
    RhoDiff : TsubSurfaceRhoBaseClass; // diffuse subsattering instance
    transparent : Boolean;         // calculate also MuMat for Transmission
    MMM : TMuMatMap;                // Müller matrix calc. module
    constructor Create(setname : String);
    procedure SetupMMM(const calcBiCu: Boolean);
    function InitMatofNKdatFile(const fn: String; const 165oolean165e: 165oolean): Boolean;
    function InitMatofConstNK(const n0,k0 : Double; const 165oolean165e: 165oolean) : Boolean;
    destructor Destroy; override;
end;
```

<figure><img src="image_108.png" /><figcaption>Figure 103 Overview of main components of TMaterialInterface</figcaption></figure>

As can be seen in the class definition, the material interface comprises the three main elements that allow the definition of a generic material based on the models discussed in Chapter 4 and depicted in the overview in Figure 103. The fundamental elements of the class are:

1. Fresnel module: the necessary parameters for the Fresnel approach are covered by the refractive index classes (`TRIbaseclass`) (section 4.6) of the two materials and any existing thin-films, as well as the `müller matrix` component `TMuMatMap`.

2. Subsurface scattering: the `RhoDiff` encapsulates the desired model and relevant parameters to perform the diffuse sub-subsurface reflection (section 4.11).

3. Surface roughness: the `microfacet` roughness model (section 4.10) is simply defined by the parameters `Rcorrelated`, `Ralpha_x` and `Ralpha_y`.
The method SetupMMM initialises the Müller matrix calculations (see section 4.9). For materials with high absorption, the Boolean parameter *transparent* can be set to *true* in order to absorb any rays transmitted below the surface immediately. Otherwise, the ray would be traced through the object and set to “absorbed” only at the next collision point. Therefore, defining the material as transparent increases the efficiency of the algorithm. By utilizing these models, the TRCMaterialInterface is able to model a large variety of opaque or transparent materials with different roughness and subsurface scattering behaviours. A typical scattering profile for an opaque material, representing a rough painted surface, that combines surface-roughness and subsurface scattering is depicted in Figure 104.

<figure><img src="image_109.png" /><figcaption>Figure 104 Typical scattering on an opaque, slightly rough painted surface ($\alpha$ = 0.1, $n$ = 1.5, $\rho_{ss}$ = 0.8)
left: polar plot showing BRDF (bidirectional reflectivity), right: polar plot showing radiance of reflection</figcaption></figure>

As can be seen, the parameterized material model is able to generate typical BRDF functions, including specular and diffuse components. The Fresnel reflection on the microfacet surface results in a polarisation of the specularly reflected beam, while the diffuse subsurface reflection is unpolarised (see also section 4.3.3).

The material interface class is primarily used to initialise the associated surface scattering TLSISroughPol module (see below). So far, the options provided by the combined model have been sufficiently versatile to model all relevant materials (coated surfaces, painted surfaces, metallic surface, coated and uncoated glasses). If necessary, alternative material interface definitions can be added in future versions of the raytracer.
## 7.8. Scattering classes (TLightSurfaceInteractionStokes)

The implementation of the raytracer follows a strict object-oriented programming hierarchy. This means that an abstract base class defines the essential functionality that has to be provided, whereas the descendants of that base class carry the actual implementations based on different physical models.

In order to perform the scattering of the traced ray at material interfaces, the abstract (=empty) base class `TLightSurfaceInteractionStokes` is introduced. Its straightforward definition is:

```typescript
Type TLSIevent = (lsiReflected, lsiAbsorbed, lsiTransmitted);

Type TLightSurfaceInteractionStokes = class (TObject) { TLightSurfaceInteractionStokes }

private
    LSIid: Integer;

public
    constructor Create;
    function init(params : TDoubleArray) : Boolean; virtual; abstract;
    function Process (var ray : TStokesRay ) : TLSIevent; virtual; abstract;
    function dumpParams : String; virtual; abstract;
    destructor Destroy; override;
end;
```

As can be seen in the definition above, apart from an initialisation function, each descendant class of `TLightSurfaceInteractionStokes` essentially needs to implement a `Process` method. The method receives the current ray object, which, in turn, contains information on the current and previous collision point (see section 7.3.1). The `Process` function performs all evaluations and samplings processes that are required to determine the outcome of the current scattering event.

Depending on the three potential outcomes and result values of the function, the following manipulations of the current ray have to be performed by any implementation of the `Process` function:

* **lsiAbsorbed**. The ray is absorbed, and the raytracing can be terminated. No further calculations of the `Process` function have to be performed, as the collision object already contains all information required for further processing (e.g. material etc.).

* **lsiReflected**. The ray has been reflected off the surface. The ray origin point (rp) is set to the collision point, and the outgoing direction (rdir) is provided. Further, a new polarisation state defined by the Stokes vector (pState) and polarisation reference direction vector (spol), have to be set.

* **lsiTransmitted**. All information regarding the refracted ray's new origin, direction and polarisation has to be set. Beyond that, the entered medium's absorption coefficient has to be provided (enteredTransparentMediumAlpha).
The currently implemented descendant classes of `TLightSurfaceInteractionStokes` are listed in Table 14. For reference, each class name begins with `LSIS` for *Light Surface Interaction Stokes*, as an early implementation of the raytracer did not consider polarisation.

<table><caption>Table 14 Currently implemented descendant classes of `TLightSurfaceInteractionStokes`</caption><thead><tr><td>class name</td><td>description</td></tr></thead><tbody><tr><td>TLSISidealBlack</td><td>Absorbs any incident ray. Mainly used for the transmittance measurement surfaces.</td></tr><tr><td>TLSISidealTransparent</td><td>Transmits any incident ray with unaltered direction and polarisation state. Mainly used for transmittance measurements.</td></tr><tr><td>TLSISsimpleLambert</td><td>Implementation of simple, ideally diffuse Lambert surface. Returns unpolarised radiation with cosine weighted directions. Controlled by a wavelength-independent reflectance parameter.</td></tr><tr><td>TLSISsimpleSchlickFresnel</td><td>Implementation of the simple, in computer graphics often used, Schlick approximation (Schlick, 1994) that describes the angular dependence of the reflectance. The implementation allows an additional diffuse component but is wavelength independent.</td></tr><tr><td>TLSISsimpleGlass</td><td>A wavelength-independent implementation that allows to describe the angular reflectance based on polynomial fits.</td></tr><tr><td>TLSISroughPol</td><td>The key scattering class for opaque and transparent materials covering Fresnel reflection, microfacet roughness and subsurface scattering considering spectral and angular dependencies. (see below)</td></tr></tbody></table>

Only the classes `TLSISroughPol`, `TLSISidealBlack` and `TLSISidealTransparent` are relevant for all validation and application cases covered in chapters 8 and 9. While the first is used to model any physical materials, the latter two are used to define surfaces used for measurement purposes (see section 7.10.4). The other classes were implemented for testing and additional validation purposes during implementation only.

Following the logic of object-oriented programming, further scattering classes can be implemented straightforwardly, inheriting definitions either from the base class or from its descendants.

## 7.9. The TLSISroughPol scattering class

The scattering class `TLSISroughPol` is currently the key class used to simulate scattering events. It was used to model all materials relevant to this work. It represents a generic surface model that is able to model many different types of material interfaces ranging from transparent glass materials to opaque, coated surfaces. The method is usually initialised by providing the material interface data in the form of a `TRCMaterialInterface` object (see section 0) and specifying the reference direction of incidence (frontside or backside). As discussed above, the model essentially comprises a Fresnel module, a microfacet roughness module and a subsurface scattering module.
<figure><img src="image_110.png" /><figcaption>Figure 105 Process flow diagram of TLSISrouhPol showing the processing of a scattering event</figcaption></figure>

An overview of how scattering events are processed by the **Process** method of the **TLSISrouhPol** class is depicted in Figure 105. The key processing steps, as indicated in the diagram, are:

(1) Using the parameters `enteredTransparentMediumAlpha` of the current ray (see section 7.3.1) and the distance between the current and the previous collision point, the probability for absorption is determined using equations (25) and (26). A potential absorption is sampled based on the integrated attenuation calculated along the path. In the case of absorption, the process halts (A1).
(2) The direction of the incoming ray and the reference vector for the polarisation state is rotated to the local (macroscopic) coordinate systems (see section 7.3.2) of the face hit. In the local system, the z-axis represents the surface normal vector $\vec{n}$.

(3) If the roughness parameter of the surface is greater than zero, the orientation of a hit micro-facet is sampled using the function `sampleGGXVNDF`. Based on the direction of the incoming ray and the surface roughness parameter, the function returns a random surface normal vector $\vec{m}$ of the microfacet (see section 4.10.2). If the roughness parameter is zero, an ideally flat surface is assumed, applying $\vec{m} \equiv \vec{n}$.

(4) To correctly perform the Müller matrix operations, the Stokes vector coordinates have to be transformed into the reference system of the microfacet (with $\vec{m}$ as its z-axis). The Stokes vector transformation is performed by determining the rotation angle and then applying equation (67).

(5) The locally aligned Stokes vector of the incoming ray $\vec{S}$ is multiplied with the Müller matrices for reflection and transmission (equations (64) and (65)) in order to determine the Stokes vectors of the reflected ($\vec{S}_R$) and transmitted ($\vec{S}_T$) ray.

(6) The total irradiance of the reflected and transmitted ray is contained in the first coordinate ($S_0$) of the related Stokes vectors. The parameters represent the Fresnel coefficients. Following the MC approach, they are used as probabilities to sample the outcome of the scattering at the microfacet:

$$ \text{probability for reflection: } R = \vec{S}_R[0] $$

$$ \text{probability for transmission: } T = \vec{S}_T[0] $$

$$ \text{probability for absorption: } A = 1 - R - T $$

The absorption is determined based on energy conservation. It reflects the absorption occurring at the material interface and will, therefore, only have a non-zero value if thin-films are present.

(7) Using a uniformly distributed random number and the determined probabilities R, T and A, the outcome of the scattering event is determined.

(8) The direction of the reflected ray is determined by applying the law of reflection, considering the incidence direction and the orientation of the current microfacet.

(9) Since no path splitting is applied (see section 7.1), the reflection path is followed with the full initial intensity. Hence, the Stokes vector of the outgoing ray is defined as the normalised Stokes vector of the reflected ray (equation (62)).

(10) By calling the shadowing function for the chosen, correlated or uncorrelated surface roughness model (see section 4.10), the probability for an unobstructed path in the outgoing direction is determined. If the sampling determines that a further microfacet obstructs the outgoing path, the process returns to step (3) and determines a new microfacet based on the current direction. If the path is unobstructed, the process exits with a reflection event (R1). In this case, the reflected beam is polarised.

(11) Based on the current wavelength and the applied subsurface reflection model (see section 4.11.2), it is determined if the ray is scattered back directly below the surface.

(12) The subsurface-scattered ray is assumed to be unpolarised and diffusely reflected. The outgoing direction is sampled accordingly, and the scattering process is terminated (R2). Further
reflections or refractions on the microfacet structure are not considered, as it is assumed that
the occurring subsurface scattering is the final result of multiple scattering events.

(13) If the material is defined as opaque by settings its transparent property to false, any ray penetrating below the surface is absorbed, and the raytracing terminated (A3).

(14) The new direction of the transmitted ray is determined, considering refraction at the micro-facet. The angle of refraction is a by-product of the Müller-matrix generation in step (5). It is used to generate the direction vector of the outgoing ray.

(15) In the rather exceptional case that the microfacet refracts the light towards the incidence medium again, it is assumed that after one or more further scattering events, the ray will finally leaves the surface again. Therefore, this case is treated like the subsurface scattering event.

(16) In analogy to subsurface-scattering (12), the outwards refracted ray is assumed to leave the surface unpolarised with a direction following the ideal diffuse distribution.

(17) Since no path splitting is applied (see section 7.1), the transmission path is followed with the full initial intensity. Hence, the Stokes vector of the outgoing ray is defined as the normalised Stokes vector of the reflected ray (equation (62)).

As can be seen in the diagram, the *Process* function can terminate in seven different ways:

(A1) The ray is (was) absorbed by attenuation in the medium on its path towards the current collision point.

(A2) The ray is absorbed within the thin-flim stack applied on the surface.

(A3) The ray is absorbed after being transmitted below the surface.

(R1) The ray is reflected on the flat surface or a microfacet. The ray is polarised.

(R2) The ray is reflected on multiple scattering events below the surface. The ray is fully depolarised and diffusely reflected.

(R3) The ray is refracted towards the incidence medium within the microfacet surface. The ray is fully depolarised and diffusely reflected.

(T) The ray is transmitted, i.e. refracted into the second medium. The ray is polarised.

The function exits and returns the result value `lsiAbsorbed`, `lsiReflected` or `lsiTransmitted`. If the ray is reflected or transmitted, the polarisation, direction and origin of the `TStokesRay` object that was passed to the function are updated accordingly.
# 7.10. Directional scanning process

## 7.10.1. Overview

In the directional scanning process, the scattering characteristic of the target object is determined for various target quantities. For the applications addressed by this work, the target quantities will mostly be the power absorbed in various elements of the object as well as the power transmitted through the target object. Based on the principle of the linearity of light (see section 4.2.2), it is possible to capture the optical behaviour of the target object using parallel beams of light coming from different directions. This angularly resolved information can later be superimposed to determine the target quantities for any given radiance profile (see section 6.2).

Technically, the scanning process can be characterised as a Monte-Carlo forward raytracing process. This means that random rays, or light particles, are cast according to their actual direction of propagation. Following the Monte-Carlo approach, the wavelengths of the cast rays as well as all physical models governing the absorption and scattering of light are modelled as stochastic processes, i.e. random samples representing the underlying physics are used throughout to model all processes. In line with the Monte-Carlo approach, the raytracing process is iterated for a large number of samples, and the resulting mean values of the target quantities and their standard deviations are tracked. By applying the central limit theorem, see equation (91), these values can be used to determine the required number of iterations for achieving the results with specified precision.

<figure><img src="image_111.png" /><figcaption>Figure 106 Schematic overview of the directional scanning process</figcaption></figure>

An overview of the process is depicted in Figure 106. Essentially, the process involves the following steps:

<figure><img src="image_112.png" /></figure>

① A distinct direction of the incidence beam is determined using the Fibonacci hemisphere (section 5.1). A large number of rays/light particles are defined. Their wavelengths are determined based on the global radiation spectrum (see section 4.4), or any other given spectral density function. Their initial polarisation state is generally defined as unpolarised, but can also be set to a specific polarisation state.
② The equienergetic light particles are cast from a caster object to impact the target with a specified and uniform flux density (see below).

③ The light particles are scattered, i.e. reflected or absorbed, on the surfaces of the target object. All scattering processes follow Monte-Carlo implementations of physical-optics-based models (chapter 2). The particles always maintain their full initial energy (or flux density), or, in other words, the path of the light ray is not branched into different subpaths. Instead, at any scattering event that partially absorbs or reflects a ray, one potential outcome is chosen following the probability distributions of the underlying physical processes.

④ Since equienergetic particles are used throughout the process, the powers absorbed, reflected or transmitted by specific surfaces of the target objects can be determined merely by counting the particles. These values provide essential information for later modelling of the energy balance of the analysed object.

⑤ A catcher object finally captures any rays that are transmitted through the target object. The catcher object can be seen as an ideally black surface that absorbs any light particles. However, due to its virtual nature, the power absorbed by the catcher is considered a transmitted power. Alternatively, a set of catcher surfaces can be used to capture the transmitted power's directional distribution (e.g. by again using the Fibonacci Hemisphere or a Klems segmentation).

## 7.10.2. Caster/Catcher concept

### 7.10.2.1. The caster concept

In the raytracing scanning process, the target is illuminated by parallel light beams coming from all sampling directions. In the raytracing approach, parallel light beams are represented by rays from randomly evenly distributed starting points with the same directional vector. The starting point distribution has to fulfil two criteria: first, the origin points must be uniformly distributed on a plane perpendicular to the direction with a predefined density. Second, the location of the starting points should guarantee that all relevant surfaces of the target are hit, but only a few, or in the best case, none, should miss the targeted area.

In order to achieve these goals efficiently, the caster concept was conceived. The caster represents an object that wraps around the targeted area of the target object. The caster object can be of random geometry as long as it fulfils certain criteria:

* the entire target area has to be covered

* the target area, together with the caster object has to form a closed object

* the caster object has to be convex so that it will not block itself, i.e. rays cast from any faces of the caster object that potentially hit the target area will not be obstructed by other faces of the caster.

* for identification, the object should be named “CASTER” in the CAD model

* the surface normals of the faces that constitute the caster have to point towards the target, i.e. towards the inside of the caster object.

Some potential geometries of caster objects are depicted in Figure 107.
<figure><img src="image_113.png" /><figcaption>Figure 107 Appropriate geometries for caster objects</figcaption></figure>

### 7.10.2.2. Generating starting points on the caster surface

As discussed above, the task is to generate randomly distributed starting points on the caster so that the associated flux density perpendicular to this direction is uniform and matches a defined value. Following these criteria and in the context of the MC approach, this transforms into two stochastic tasks:

1. face-sampling: randomly determine a face of the caster object

2. point-sampling: randomly choose a starting point on the face

In order to pick a random face efficiently, some auxiliary variables are determined before the distribution is generated. First, the probability of selecting a specific face $i$, with surface normal $\vec{n}_i$ and area $A_i$ is calculated. The probability for this face is determined by the radiant flux $\Phi_i$ (or radiant power) that the face contributes to achieving a constant flux density $E$ in the specific direction $\vec{v}$:

$$ \Phi_i = E \cdot A_i \cdot \max(\vec{n}_i \cdot \vec{v}, 0) \quad (130) $$

As $E$ is constant for all faces, the probability of selecting a face is proportional to the weighted projected area, or in this case more accurately, the projecting area of the face. The maximum function ensures that only faces pointing towards the direction of the light are considered. Using equation (130), the radiant fluxes are calculated for all faces forming the caster. The results are cumulated successively in an array. In the actual implementation, however, already the proportional quantity, density of rays per watt and square metre, denoted as $P$, is considered:

```pascal
procedure TVGeo.InitCaster(PperWm2: Double; Dir: TVect3d);
var fc: Integer;
    d: Double;
begin
  With caster, target do
  Begin
    rdPperWm2:=PperWm2;
    rdDir:=VNormalize(Dir);
    rdTotalP:=0;
    SetLength rdFaceP,length rdFaceId); // init arrays rays per face
    SetLength rdFacePcum,length rdFaceId); // init array cum. rays per face
    For fc:=0 to high rdFaceId do // iterate through all faces of caster
      with obj[rdObjId].face[rdFaceId[fc]] do
      Begin
        d:=Vdot (rdDir,snv); // dot product
        if d>0 then rdFaceP[fc]:=(rdFaceSz[fc]*rdPperWm2*d) // rays per face
        else rdFaceP[fc]:=0;
        rdTotalP:=rdTotalP+rdFaceP[fc]; // cumulate rays per face
        rdFacePcum[fc]:=rdTotalP;
      end;
    end;
  end;
end;
```
The array rdFacePcum is part of the caster variable and now contains real values indicating the cumulated assigned numbers of rays per face resulting from the chosen direction and chosen numbers of rays per watt and square meter. The array would represent the cumulated distribution function if the results were normalised. The variable rdTotalP is identical to rdFacePcum of the last face and indicates the number of rays that have to be cast in total from the caster object to achieve a flux density of 1 W/m². A random real number t in the range of [0,rdTotalP] can now be used to randomly determine a face of the caster to initiate a ray. This is done by finding the corresponding bin number c of the array rdFacePcum so that rdFacePcum[c]<t< rdFacePcum[c].

Once the face is determined, the second task is to pick a ran- domly distributed point on its surface. Since a face can represent either a quad or a triangle, two different sampling algorithms are required. The algorithm for the quad is trivial. Based on two uniform random variables u and v and the face vector defini- tions (see Figure 96), the coordinates of the random point ⃗ p can be determined by:

$$
\vec{p} = \vec{v}_0 + u \cdot \overrightarrow{edge1} + v \cdot \overrightarrow{edge2} \quad (131)
$$

In order to uniformly sample triangular faces, a short derivation
is required. *u* and *v* are again the coordinates of the point; how-
ever, they are no longer independent of each other. Firstly, the
*u* coordinate is considered. If the points should be distributed
evenly, it is evident that the probability distribution function
(*PDF*) of *u* is linearly decreasing to zero along the edge1 axis as
the reachable points decrease towards the tip of the triangle (see
Figure 109). Adding a factor of 2 for normalisation, the *PDF*
function for *u* is obtained, see equation (132). By integration,
the cumulated distribution function (*CDF*) is determined, see
equation (133). This function is inverted in order to be able to
determine *u* based on a uniformly distributed random number
ξ₁, equation (134). Once the *u* coordinate is determined, the *v*
coordinate is uniformly distributed (αξ₂), but needs to be scaled
down according to its location regarding the *u* coordinate, equa-
tion (135).

$$
PDF(u) = 2(1 - u) \tag{132}
$$

$$
CDF(u) = 2u - u^2 \tag{133}
$$

$$
u = CDF^{-1}(\xi_1) = 1 - \sqrt{1 - \xi_1} \quad (134)
$$

$$
v = \xi_2 (1 - u) \tag{135}
$$

Based on these definitions for *u* and *v*, a uniformly sampled point on the triangle can be found by using equation (131) again. Both sampling algorithms have been tested by using different caster geometries and analysing the uniform distribution of the rays when projected on a screen.

<figure><img src="image_114.png" /><figcaption>Figure 108 Uniform sampling
of quad face</figcaption></figure>

<figure><img src="image_115.png" /><figcaption>Figure 109 Uniform sampling
of triangular face</figcaption></figure>
### 7.10.2.3. The catcher object

The concept is completed by a second complementary object that is used to catch any rays not absorbed by the target object. The requirements for this object are, however, less complex. Actually, the only criterion is that it should wrap the geometry so that it can absorb any rays that were not absorbed by the target or the caster object. Of course, the catcher object should not penetrate or interfere with the target object or caster in a way that it captures rays that are still relevant for the raytracing process. Again, the face normals of the catcher object in the CAD model should be oriented inwards. The distinct face orientation allows raising warnings in case the catcher object is hit from the outside, which would then indicate an error in the geometry or raytracing process.

## 7.10.3. Incidence beam - starting a ray

The direction of the incidence beam is determined based on the chosen method: the N discrete sampling directions can either be provided by a Fibonacci hemisphere (see section 5.3) or, alternatively, by the Klems segmentation. On every iteration, the given direction is provided to a method called `StartRayDir`, which initialises the remaining parameters that are required to start a ray (light particle). The following steps are required to launch a single ray:

1. Determine starting point

This is done based on the algorithms explained in the preceding chapter (see 7.10.2). Firstly, a subsurface of the caster object is determined, and then a random point on the face is picked.

2. Initialize wavelength

In the standard process, the function `getSampleWLat` (see 4.4.2) is called to generate a random wavelength sample based on the global radiation spectrum. Alternatively, for special applications, any other spectral distribution could be considered. Beyond that, the scanning process could also be performed for monochromatic samples, particularly for three wavelengths representing the (R, G, B) values for lighting applications.

3. Initialize polarisation state

For most applications, the polarisation state will be set to the unpolarised Stokes vector (see Table 6). Again, alternative definitions for other application fields are possible (e.g. for performing virtual ellipsometric measurements, see, e.g. validation examples in section 4.9.3). Further, the polarisation reference system has to be initialised by providing an axis for unpolarisation. In the case of unpolarised radiation, the axis can be randomly chosen but has to be defined perpendicular to the direction of light.

4. Initialize absorption coefficient

The volume absorption coefficient has to be provided for the first path of the light particle after its creation. Later this information will be provided by the light-surface interaction modules. Since no atmospheric absorption is considered and the rays will usually be launched in the air medium, the value is initially set to zero.

5. Initialize collision trackers

The trackers for the current and recent collision points have to be initialised. While the recent collision is defined as unset, the current collision point references the relevant surface of the caster object. This is necessary to exclude the starting surface from the first collision detection cycle. Otherwise, the starting face could be erroneously detected as the first collision.
The implementation of the ray initialisation routine is provided in full in the code snippet below.
The initialisation steps are carried out in the same order as specified above.

```pascal
procedure TRCScan.StartRayDIR;
  var t: Double;
    c: Integer;
    u, v: Double;
begin
  With geo.radi do
  Begin
    t:=random*rdTotalP;

    // speed up and avoid overflow/range check, check last bin first
    if (t>=rdFacePCum[high(rdFacePCum)-1]) then c:=high(rdFacePCum)
    else Begin // find bin=starting face on caster object (radiator)
      c:=0; while (t>=rdFacePCum[c]) do Inc(c);
    end;

    with geo.target.obj[rdObjId].face[rdfaceid[c]] do
    Begin // uniformly sample starting face
      u:=random; v:=random; // two random numbers
      if is4p=false then // triangle sampling
        Begin u:=(1-sqrt(u)); v:=v*(1-u); End; // v=0..u;
      ray.rp:=v0+ (u*edge1) + (v*edge2);
      ray.rdir:=rdDir;
    end;

    with ray do Begin
      // init coll. counter
      wl:=GR.getSampleWLat(random); // get sample of global rad. specttrum
      pState:=unpolstokes; // start with unpolarised stokes vector
      // chose initial reference axis for Stokes vector
      if abs(rdir.z)<0.99 then spol:=VNormalize(VCrossProd(ZAxis,rdir))
        else spol:=VNormalize(VCrossProd(YAxis,rdir));
      curColl.sprod:=NoColl; // init dot prod.
      curColl.oNo:=MaxUInt32; // init hit object no.
      curColl.cp:=ray.rp; // set ray origin
      curColl.enteredTransparentMediumAlpha:=0;
      lastColl.sprod:=NoColl; // int dot product
      lastColl.oNo:=rdObjId; // point on caster, to avoid self-collision on start
      lastColl.Fno:=rdfaceid[c]; // point on caster face, to avoid self-collision on start
      lastColl.enteredTransparentMediumAlpha:=0; // no absorbtion on first path
      CollN:=0;
    end;
  end;
end;
```

### 7.10.4. Measurement of target quantities

In order to provide a versatile method to perform measurements in the raytracing scan, measurement functions can be assigned based on the surface's material assignment. This allows a flexible approach, as the user can freely define measurement regions on any object by assigning specific materials. While these materials can still be based on identical physical properties, the separation allows for performing individual measurements. Currently, only energetically relevant measurement functions are implemented. However, the evaluation of other quantities, such as degree of polarisation, angle of polarisation or spectral properties, can easily be added in future.

The intended measurements are defined in the initialisation of the scanning process using the function `addPowerMeasurement`. The function tracks the absorbed, transmitted or reflected power for either a single material or cumulated for two different materials. The data is stored as an array for all N directions of the scanning process. Later on, this sampled data is used to generate the *SIOP* (see section 6.1.2).

Each measurement is initialised by providing the following arguments for the `addPowerMeasurement` function:

* `mId[1,2]`
The material identifier is used to assign the measurement function to a specific material. If the value -1 is passed for the second material, only one material is considered in the measurement.

* `lsiType[1,2]`

Specifies the “light-surface-interaction” type that should be measured. The options `lsi-Absorbed`, `lsiTransmitted` or `lsiReflected` are used in order to determine which power should be tracked on the material.

* `inciType[1,2]`

This parameter is used to specify the incidence direction for the measurements. If `itFront` or `itBack` is selected, only hits from either the frontside or backside will be considered. If `itBoth` is chosen, all hits of the surfaces are tracked.

* `referenceArea`

A reference area can be provided to convert the measured powers into irradiance (flux density) values [W/m²]. If no reference area is provided, the results are considered to be of type power [W].

* `filename`

Finally, a filename for the sampled measured data is provided. The measured data, along with specifications regarding the scanning process, are stored in the specified file.

The function can best be explained with the help of an example. Below is a typical example for the initialisation of measurements used to evaluate the performance of a shaded window detail.

```
{1} addPowerMeasurement(mLib.MatByName('ABSORBER_INT'), -1, lsiAbsorbed, lsiAbsorbed,
                         masterName+'-ABS_INT', itFront,itFront,1);
{2} addPowerMeasurement(mLib.MatByName('PROBE_CC'),-1,lsiTransmitted,lsiTransmitted,
                         masterName+'-PROBE_CC',itFront,itFront,0.25*0.25);
{3} addPowerMeasurement(mLib.MatByName('SLAT'), -1, lsiAbsorbed, lsiAbsorbed,
                         masterName+'-SLATS', itBOTH,itBOTH,1);
{4} addPowerMeasurement(mLib.MatByName('PANE1_e'), mLib.MatByName('PANE1_i'),
                         lsiAbsorbed, lsiAbsorbed, masterName+'-PANE1abs', itBOTH,itBOTH,1);
{5} addPowerMeasurement(mLib.MatByName('PANE2_e'), mLib.MatByName('PANE2_i'),
                         lsiAbsorbed, lsiAbsorbed, masterName+'-PANE2abs', itBOTH,itBOTH,1);
{6} addPowerMeasurement(mLib.MatByName('PANE3_e'), mLib.MatByName('PANE3_i'),
                         lsiAbsorbed, lsiAbsorbed, masterName+'-PANE3abs', itBOTH,itBOTH,1);
```

The sequence defines the measurement of the total transmitted power (1), the irradiance transmitted through a specific probe area with dimension 25 cm x 25 cm (2), the power absorbed by the slats of the shading (3), and the (bidirectionally) absorbed powers in the three panes of the window (4,5,6). A schematic overview of the six measurements is depicted in Figure 115.

As can be seen in the initialisation code, the type `lsiTransmitted` is used for the measurement in the

<figure><img src="image_116.png" /><figcaption>Figure 110 Schematic overview of typical measurement configuration</figcaption></figure>
probe area, as the material of the probe area is defined as `TLSISidealTransparent` (see Table 14).
However, the type `lsiAbsorbed` is used for the internal absorber object (catcher), as its material is defined as `TLSISidealBlack`. Both measurements, however, are actually used to determine transmitted radiation. This approach ensures that the probe can be used for measurements restricted to a specific region but does not distort the results of the internal absorber. Furthermore, both types are limited to frontside incidence using the `itFront` specifier, as no incidence from the backside is expected.

In order to determine the absorption of the three window panes and slates, the specifier `lsiAbsorbed` is used. To further consider the absorption of any incidence direction, the specifier `itBoth` is chosen. While the slats are defined using a single material, two materials are used for each glass pane. In the example, the absorptions on each side of the panes are cumulated into one measurement by referencing the materials of both sides within a single measurement (`PANE1_e` and `PANE1_i` etc.).

While the provided example is kept simple, the material-based measurement concept can be used to perform refined measurements in specific areas, e.g. distinct parts of the window frame or segmented or stratified regions of the glazing.
### 7.10.5. Precision and convergence

In order to provide the desired accuracy for the raytracing results, the number of iterations, i.e. cast rays, is controlled and adapted for each direction individually.

The statistical accuracy of the result can be determined by tracking the standard deviation of the result and applying the central limit theorem (see section 5.1.3). For this purpose, one relevant measurement (e.g. total transmitted power) has to be selected, and a precision target for this quantity has to be specified. The precision target is defined as a threshold value below which the determined standard error of the mean $\sigma_{SEM}$ of the selected quantity must fall. It can be provided either as a target for the absolute error or as a target for the relative error. For most applications using absolute precision targets is more appropriate, as usually a suitable reference value for deriving the precision target, e.g. the power of the incidence radiation, is available. The latter approach, providing relative targets, can be problematic for some applications if the target quantity is marginal or zero for some directions.

As discussed in section 5.1.3 the error estimators are derived from sampled means values. In order to achieve this, mean values of cycles with constant numbers of iterations are tracked. In the example depicted in Figure 111, each sampling cycle consists of 10,000 single raytracing processes. After each cycle, the mean value of all cycle results ($q$) and the standard deviation ($s$) are computed. Finally, equation (91) is used to evaluate the standard error of the mean $\sigma_{SEM}$, which is then compared against the precision target. The iteration is terminated if it falls below this threshold for four consecutive cycles. In the example provided in Figure 111, the precision target of 0.0025 was reached after 230 cycles, equivalent to 2.300.000 million single raytracing processes. A false, early termination can be triggered during the first few cycles when the statistical evaluation is still unreliable. In order to prevent this, the error tracking is suspended for an initial period with a constant number of iterations.

<figure><img src="image_117.png" /><figcaption>Figure 111 Typical convergence during raytracing top: target quantity with confidence intervals bottom: standard deviation and standard error of the mean</figcaption></figure>
## 7.11. Rendering and backward raytracing

Unlike tracing rays in their natural direction (forward raytracing), as applied for the *RadiCal* method, it is also possible to perform backward raytracing. As stated above, most raytracing applications actually perform backward raytracing, as their purpose is the generation of (more or less) photorealistic images of three-dimensional objects (renderings).

Following the fundamental principle of the reversibility of light propagation, referred to as Helmholtz or Stokes-Helmholtz reversibility, it is always possible to reverse the path of light, applying the same models to describe any refractions or reflections along this path. This does include the detailed modelling based on polarisation states, as applied here. Therefore, the *RadiCal* raytracer can be used to generate images of the target objects if rays originating from a specific eye point are directed towards the target object.

Unlike in common raytracing tools, no explicit light sources are included in the model. Instead, a high-contrast spherical environment provides illumination for the rendering. This is because considering light sources with a singular location (point light source) or a singular directional profile (sunlight) requires different sampling techniques. In such cases, biased sampling or correlated sampling strategies (e.g. Metropolis sampling) have to be applied to ensure that the relatively few relevant directions having high intensities are sampled adequately; see, e.g., Veach *et al.* (1997).

However, the forward raytracing approach implemented in the *RadiCal* method relies on unbiased, uncorrelated sampling, which is more efficient for the chosen forward raytracing approach. Since the main purpose of generating renderings with the *RadiCal* method is validating the applied algorithms, it would not be reasonable to implement alternative sampling strategies for backward raytracing. Nevertheless, the implementation has proven that unbiased sampling can still be used for rendering purposes if singular light sources are avoided.

Instead, illumination is provided by a panoramic image in HDRI format. HDRI stands for *high dynamic range image*. Unlike in regular image formats, the RGB values are not stored as either 8 or 16-bit integer values but as 32-bit floating point values. While the colour channels with integer values are

<figure><img src="image_118.png" /><figcaption>Figure 112 Typical full-panoramic HDRI environment
CC0 license (Poly Haven - public asset library, 2023)</figcaption></figure>

<figure><img src="image_119.png" /><figcaption>Figure 113 Working principle of backward raytracing
(idealised, environment at infinite distance)</figcaption></figure>
limited to 1:255 (8-bit) and 1:65535 (16-bit) contrast ratios, the floating-point exponent format allows almost infinite contrast ratios. This can be used to include extended bright regions that effectively provide illumination in the *RadiCal* rendering method. In order to optimise the efficiency of the sampling, i.e. the calculation time for the rendering, these illumination spots should still be extended over relatively large areas of the image. Consequently, the final renderings appear as photos taken under diffuse lighting conditions.

The principle of the backward raytracing algorithm is depicted in Figure 113. A virtual screen is located between the eye point and the target object. The screen represents the final rendered images. Rays originating at the eye point are cast through the screen towards the object, successively scanning each pixel of the image. If the rays hit the target objects, the specific scattering functions (see section 7.8) are called to determine if the ray is absorbed, reflected or refracted. The process is repeated until the ray is absorbed or left without further intersection. The final direction of all rays that are not absorbed is used to determine a corresponding location on the environmental image (HDRI). The colour and brightness of the relevant location on the (HDRI) finally determine the contribution of this ray to the colour of the pixel in the screen.

Figure 114 and the related description provide additional information regarding the actual implementation of the algorithm. All parts of the algorithm were developed and implemented by the author. However, they cannot be presented in more detail here, as this would exceed the scope of this chapter.

(1) A sample ray originating at the eye point and going through the screen pixel (x,y) is generated. The polarisation state is generally assumed as unpolarised. The wavelength is randomly sampled on the global radiation spectrum but limited to the visible range.

(2) The collision detection algorithm is called (section 0) to detect intersections with target objects.

(3) If a surface is hit, the light-surface-interaction class is called to determine the result of the scattering event. If the ray is refracted or transmitted, the collision detection is repeated. If the ray is absorbed, the tracing of the current ray is terminated, and a new ray is generated.

(4) The wavelength of the sampling ray is converted into an

<figure><img src="image_120.png" /><figcaption>Figure 114 Process flow diagram of the RadiCal backward raytracing/rendering algorithm</figcaption></figure>
RGB vector. Note that the mapping of wavelengths into the RGB colour space is not unique. Therefore, the conversion provides just one plausible solution.

(5) Based on the final direction of the ray, a position on the environment HDRI is determined. The corresponding RGB colour vector is then determined by calculating a weighted average of the neighbouring four pixels (anti-aliasing).

(6) The components of the environment colour vector are weighted using the components of the $\lambda[r,g,b]$ vector. The resulting RGB vector is cumulated on the screen pixel.

(7) The convergence is monitored by applying a simple colour metric for the standard errors of the RGB channels of the screen pixel. If the determined colour deviations fall below a specified threshold, the sampling for this pixel is terminated. The standard error of the colours is again determined by applying the CLT (section 5.1.3). The approach allows for minimizing the sampling steps required for each pixel.

The entire process is carried out for all pixels of the screen. The final image is further processed using commonly applied colour and brightness calibration as well as a gamma correction. These post-processing steps are not part of the description above. Depending on the characteristics of the target object, the contrast of the selected HDRI environment and the targeted colour noise threshold, typically 100 to 1000 samples per pixel are required. Hence, the total number of scattering events that must be processed for a single image typically exceeds one billion operations. The occurring scattering events, therefore, cover an extensive range of potential wavelengths, angles and polarisation states. Any errors in the algorithms or numerical instabilities of the algebraic methods, such as algebraic singularities, would either cause a termination of the raytracing or faulty results in the final rendering. Generating rendered images by performing backward raytracing with the *Radi-Cal* algorithm is therefore considered a powerful validation step. In addition to validating the implemented algorithms, the renderings also allow

<figure><img src="image_121.png" /><figcaption>Figure 115 Glass/metal demo rendering using RadiCal raytracer
*Refractions, inter-reflections and total-internal reflection highlighted*</figcaption></figure>

<figure><img src="image_122.png" /><figcaption>Figure 116 RadiCal-rendering of validation window model</figcaption></figure>
checking the correct and appropriate assignment of material properties.

Exemplary renderings are shown in Figure 115 to Figure 117. The dashed circles in Figure 115 indicate regions with pronounced refraction, interreflections and total internal reflection. Figure 116 shows a rendering of the validation window model that was used for the full-system validation (chapter 8). In Figure 117, only one material's complex-valued refractive index function is altered. The material is used for the panes of the left window, the shading slats of the right window and the octopus sculpture. It can be seen how the reflectivity, refractivity, transparency and colour change significantly. Note that no explicit RGB colour information is contained in any model. The colours arise naturally from the wavelength dependencies of the refractive index and Fresnel functions.

<figure><img src="image_123.png" /><figcaption>Figure 117 Demo rendering "octopus" using RadiCal raytracing algorithm

The refractive index of the material used for the sculpture, left window pane and shades is varied: top-left: fused silica glass; top-right: float-glass; bottom-left: gold-brass; bottom-right: silver (Octopus 3D model by Lukáš Marek; CC BY-ND 4.0)</figcaption></figure>
# 8. Full system empirical validation

## 8.1. PyroScanner measurement method

### 8.1.1. Introduction

Since the *RadiCal* method aims to accurately predict energy flows induced by solar radiation for "real-world" applications, the ultimate benchmark is empirical data gathered in natural conditions. In order to perform a conclusive validation, empirical data should have the following features: variation of the shading system, variation of the shading state as well as variation of direct and diffuse irradiance components. The measured data must comprise high-quality data for internally and externally measured irradiance and resolved information regarding the sky diffuse and direct irradiance. As no data fulfilling these requirements had been available, a dedicated measurement equipment and process was developed in the course of the PhD. The equipment was assembled with the kind support of the AEE INTEC measurement department (AEE INTEC, 2023). The author implemented the required software and firmware for the various components and performed all measurements. The equipment was later used for an alternative measurement method, covered in a separate publication (Rüdisser et al., 2021).

Performing measurements under natural conditions is a valuable method to determine performance data; however, it is linked to challenging and elaborate tasks. The measurements depend on many external parameters that cannot be controlled but must be determined as accurately as possible. For the present case, it can be stated that the irradiance behind a window shaded by Venetian blinds depends on the following external, not controllable parameters:

* solar elevation angle (θ)

* solar azimuth angle (φ)

* diffuse irradiance (DHI)

* direct normal irradiance (DNI)

* global irradiance (GHI)

* albedo (albedo)

The three irradiance parameters are not independent of each other when an irradiance model is used to derive a parameter from the others. The albedo can be approximated as a constant value, although it actually exhibits a seasonal dependence as well as a dependence on solar angles and irradiances (as discussed in sections 8.1.11 and 6.2.4). The model for the irradiance represents an idealised case, as the actual irradiance profile is significantly more complex. Deviations from the modelled irradiance profile can mainly be attributed to inhomogeneous cloud patterns, varying reflections from objects on the ground, absolute air temperature and humidity levels or gradients, as well as to the concentration of scattering particles in the atmosphere.

Beyond these environmental factors, the irradiance behind the shaded window also depends on these controllable configuration parameters:

* shading type (sht)

* shading material (mat)

* window type or state (open/closed) (wst)

* slat angle (ang)

* window orientation (wor)
The usually uncontrollable parameter window orientation is here controllable since the testing facility at AEE INTEC allows rotation around its vertical axis. Still, the window orientation is considered an environmental parameter, whereas the remaining controllable parameters define the model that has to be processed using the *RadiCal* method. Each combination of controllable parameters [sht, mat, wst, ang] defines a specific test case, represented by a distinct 3D model. Following the empirical measurements, *SIOPs* derived for the specific test cases are used to perform virtual measurements based on the recorded environmental parameters. In order to analyse the accuracy of the model versus the measurement, the quantities of internal irradiance *E*<sub>int</sub> and transmittance *τ* are determined based on the model and the measurements. The transmittance relates the internal irradiance to the external irradiance. Its definitions for the model and measurements are provided in the subsequent sections. While the two quantities are not independent, they allow different interpretations regarding the model's accuracy.

Thus, for each test case, defined by a unique combination of controllable parameters, datasets represented by the two pairs form the basis to assess the quality of the *RadiCal* modelling:

$$ \begin{gathered} (E_{\text{int,MEAS,sht,mat,wst,ang}}(\theta, \varphi, \text{DHI}, \text{DNI}, \text{albedo}, \text{wor}), \quad E_{\text{int,MODEL,sht,mat,wst,ang}}(...)) \\ (\tau_{\text{MEAS,sht,mat,wst,ang}}(\theta, \varphi, \text{DHI}, \text{DNI}, \text{albedo}, \text{wor}), \quad \tau_{\text{MODEL,sht,mat,wst,ang}}(...)) \end{gathered} $$

As mentioned, the parameters albedo and window orientation are considered constants for each test case. Compared to the relevant definitions in the standards, the here-applied quantities of internally measured irradiance *E*<sub>int</sub> and transmittance *τ* are based on a more complex definition, as any real-world effects (shading and reflection by the reveal and window frame) and the significant spatial inhomogeneity of the irradiance caused by Venetian blinds are taken into account. In order to address this adequately, a novel measurement device and evaluation method have been established.

### 8.1.2. Measurement location

The measurement device is installed at the roof-top testing facility of AEE INTEC in Gleisdorf (lat: 47.1097°N, lon:15.7097°E; see Figure 118). The local altitude is 361 m MSL. The bottom side of the windows is located 10 m above the ground level (371 m MSL). The testing facility can be rotated into different compass directions. In south facing position, the distance between the testing window and the weather station is 9 m. The roof of the neighbouring building dominates the view to the south below the horizon (see Figure 120 top); however, the view on and above the horizon is mostly unobstructed, apart from one tree in close distance and a few trees at a further distance. Likewise, the west view (Figure 120 bottom left) offers a largely unobstructed view of the horizon.

<figure><img src="image_124.png" /><figcaption>Figure 118 Aerial view of measurement location
red circ.: weather station, yellow circ.: test window/PyroScanner
(source: Digitaler Atlas Steiermark)</figcaption></figure>
<figure><img src="image_125.png" /><figcaption>Figure 120 top: south panorama view; bottom-left: west panorama view; bottom-right: lab roof</figcaption></figure>

The lab's roof is surrounded by a glass railing reaching a height that is approximately 1 m lower than the bottom level of the window. The floor construction on the laboratory roof consists of steel grating panels mounted approximately 30 cm above a black roofing membrane (see Figure 120 bottom right).

<figure><img src="image_126.png" /><figcaption>Figure 119 Weather station with sun-tracker SOLYS2</figcaption></figure>
### 8.1.3. Measurement setup

The basic measurement setup consists of two window-based pyranometer units and an additional, nearby weather station (see Figure 119). The two window-related sensors comprise a stationary externally mounted pyranometer and a moveable internal pyranometer. The external sensor is of type CM11 (Kipp & Zonen, 2023a) and is mounted at a lateral distance of approx. 40 cm to the window. The device is oriented upright, i.e. parallel to the plane of the window element. The external sensor is used to determine the incident irradiance on a plane parallel to the window. Additionally, the high-performance sun tracker unit SOLYS2 (Kipp & Zonen, 2023b), part of the AEE INTEC weather station, monitors the incident diffuse and direct irradiance with high accuracy. All internal and external irradiance measurements are recorded with a 1-second resolution.

<figure><img src="image_127.png" /><figcaption>Figure 121 left: PyroScanner indoor pyranometer unit on moveable platform; right: overview of PyroScanner 3-axis router platform</figcaption></figure>

The unique feature of the *PyroScanner* is that the internal pyranometer is mounted on a platform that can move freely behind the window. This is realized by a self-developed, lightweight 3-axis leadscrew router platform, as used in CNC manufacturing. The necessary components were assembled in the AEE INTEC laboratory with the assistance of the measurement department. Many required components were created using a 3D printer. Table 15 provides an overview of the key hardware components used.

<table><caption>Table 15 Key hardware components of the PyroScanner</caption><thead><tr><td>type</td><td>model</td><td>manufacturer</td><td>Purpose/remarks</td></tr></thead><tbody><tr><td>pyranometer</td><td>SP-510</td><td>APOGEE instruments</td><td>internal, moveable pyranometer</td></tr><tr><td>pyranometer</td><td>CM11</td><td>KIPP & ZONEN</td><td>external pyranometer</td></tr><tr><td>A/D converter</td><td>ADS1115</td><td>Texas Instruments</td><td>precision analogue-to-digital converter (16-bit)</td></tr><tr><td>solar monitoring station</td><td>SOLYS2</td><td>KIPP & ZONEN</td><td>Sun Tracker as part of the AEE weather station for monitoring the incidence diffuse and direct irradiance</td></tr><tr><td>microcontroller unit</td><td>NodeMCU LUA<br>WIFI ESP8266</td><td>Espressif Systems</td><td>Microcontrollers with WiFi connectivity used for control and data acquisition in the measurement process</td></tr><tr><td>microcontroller unit</td><td>Arduino Uno</td><td>unspecified</td><td>microcontroller used for CNC control</td></tr><tr><td>CNC expansion board</td><td>CNC Shield V3.0</td><td>unspecified</td><td>stepper motor control board with drivers</td></tr></tbody></table>
To be able to perform the com-
plex measurements, a framework
of specialized software had to be
implemented. The implementa-
tion was performed in several dif-
ferent programming languages:
C++ for the microcontroller
firmware, PHP for the web ser-
vice used to control the measure-
ment, JavaScript and HTML for
the web client to monitor and
launch the measurements and fi-
nally, Pascal (Lazarus) in combi-
nation with SQL for the evalua-
tion and measurement database.

The *PyroScanner* features two

different operation modes: a point-mode allowing stationary measurements for extended time spans at changing positions and a tracking-mode in which the internal pyranometer is constantly moving, and measurements are taken on the fly. The latter operation mode is not covered in the thesis, as it was later added for a different project (Rüdisser *et al.*, 2021).

<figure><img src="image_128.png" /><figcaption>Figure 122 PyroScanner data and control flow diagram</figcaption></figure>

<figure><img src="image_129.png" /><figcaption>Figure 123 Web interface of PyroScanner to remotely control and monitor measurements</figcaption></figure>

The measurement sequences can be defined in a versatile way. The GRBL programming language, used in CNC manufacturing, was extended to additionally provide functions for controlling the pyranometer measurements. This allows the definition of specific measurement scripts to perform measurements at a series of defined locations automatically. The following parameters can be controlled individually:

- location of the pyranometer (x,y,z)

- speed and acceleration of the movement

- interval between measurements (in milliseconds)

- duration of the measurement (in milliseconds)
### 8.1.4. Calibration

A co-location calibration was performed to calibrate the voltage signal of the internal sensor against the already calibrated irradiance values provided by the external pyranometer. The co-location calibration ensures that both pyranometers receive the same irradiance, in terms of absolute value, but importantly also regarding the spectral and angular distribution. For this purpose, the internal pyranometer is mounted on a 3D-printed fixture to ensure that both pyranometers are aligned in the same plane and face in the same direction (see Figure 124).

The actual calibration turned out to be challenging, as the indoor pyranometer or its combination with the used A/D converter resulted in non-linear behaviour for low irradiance values. The values above 75 W/m², however, showed an excellent linear correlation with $R^2 = 0.9997$. In order to ensure that hysteresis effects or erroneous signal processing were not the origins of the observed, strongly non-linear behaviour, the calibration measurements were repeated for increasing (sunrise) and decreasing irradiance (sunset). Furthermore, additional measurements were taken using an artificial light source (tungsten lamp). All measurements confirmed that the calibration results were repeatable and unambiguous. Therefore, the measured data could be used to derive a calibration curve for the *PyroScanner* measurements.

<figure><img src="image_130.png" /><figcaption>Figure 124 colocation calibration of internal pyranometer</figcaption></figure>

<figure><img src="image_131.png" /><figcaption>Figure 125 Calibration curve used for evaluation irradiance measured at the external pyranometer (CM11) as a function of output signal of internal pyranometer (SP-510); dots: measured values; line: calibration curve; circles: cubic spline knots determined in the optimisation process.</figcaption></figure>

The calibration was split into a linear and a non-linear segment to account for the specific behaviour detected. The non-linear segment was extended to 125 W/m² to provide a continuous transition into the linear segment. Again, the knots of a cubic spline function were optimised to fit the underlying data by applying the simulated annealing algorithm (see section 5.2).

The approach turned out to be appropriate since the established spline function is able to fit the measured calibration data with high precision. The coefficient of determination $R^2$ for the fit of the
3203 calibration points reached a value of $R^2 = 0,99992$. Hence, the sophisticated calibration is able to compensate for the severe non-linearity of the sensor. However, this was only possible by using a high-precision analogue-to-digital converter that allows to capture the voltages with high precision. The steep slope of the signal in the range below 50 W/m² still increases the statistical errors of the measurements in this range.

The time-resolved functions of the irradiance data measured by both pyranometers during the calibration process were subsequently used to analyse a potential time shift between the measurement signals. Correlation analysis revealed that the signal of the *CM11* device trailed the *SP-510* device signal at a mean time lag of 5.5 s. This is in good accordance with the exponential time responses stated by the data sheets of the two devices. The *SP-510* specifications indicate a (not further specified) "time response" value of 0.5 s, which presumably reflects the time constant at which the sensor reaches $1 - \frac{1}{e} (\approx 63\%)$ of the final value. The *CM11* data sheet is more concise, as it specifies 15 s for the sensor to reach 95%. This can be converted into an exponential time constant of 5 s, thus validating the results of the correlation analysis. Consequently, the determined time lag was applied to compensate for the time lag between the two sensors in all measurements.

### 8.1.5. Prevention of internally scattered light

In order to focus on the window and shading, a black catcher plane (see section 7.10.2) is located directly behind the defined measurement area in the model. It ensures that any light passing the measurement plane is absorbed. A similar configuration is established for the empirical configuration by placing a black curtain behind the *PyroScanner* device. However, the curtain had to be located at an approximate distance of 1 m behind the window plane to allow free movement of the cable chains.

The remaining surfaces of the floor, walls and ceiling are covered with white boards. The light reflected from these bright surfaces is reflected back to the window or shading and further into the sensor, especially on measurements with a negative slat angle and high solar elevation (see Figure 126 right). Dark covers should be applied to these surfaces for future measurements.

<figure><img src="image_132.png" /><figcaption>Figure 126 Potential sources for internally scattered light
left: metallic lead screw and rail; right: illumination of the floor with negative slat angle</figcaption></figure>

In addition, the metallic lead screw and rail of the *PyroScanner* reflect irradiance back to the window and shading and further into the sensor. Due to the limited area, this effect is considered to be less
relevant for diffuse irradiance. However, it might occasionally affect some measurements if direct radiation is reflected at specific angles.

### 8.1.6. Determination and adjustment of slat angle

The slat angle is an essential parameter that must be determined accurately. In order to do this efficiently, a gauging tool was designed and realized using a 3D printer. The tool can be clipped tightly to the slats to indicate the plane that is formed by the slat's two edges. Photos taken at a further distance in the normal direction to the tool were then used to determine the angle of the slats by applying photogrammetry. The angles were measured against vertical elements of the window frame but are defined as angles relative to the global horizontal plane. Hence, the correct vertical alignment of the window had to be verified first by using a plummet.

The angle determination was applied repeatedly to examine its accuracy and ensure that the position of the slat could be moved into the same position with sufficient accuracy. The total error of the method is estimated to be approximately 3 degrees. The process

was carried out for several different angles to derive a function that provides the slat angle based on the motor's driving time. This function was individually determined for each slat mounted for the tests.

In order to adjust the slats to a specific angle, the shading device was first fully opened, then closed (i.e. extended to full length) and finally rotated for the specific driving time provided by the determined function. The automation system of the testing facility allows controlling the motor driving time with a temporal resolution of Milliseconds.

<figure><img src="image_133.png" /><figcaption>Figure 127 Photogrammetric determination of slat angle using a 3D-printed tool</figcaption></figure>

### 8.1.7. Evaluation method for the empirical measurements

Evaluating the transmittance for homogenous irradiance profiles, e.g. behind unshaded windows, is trivial as the measured transmittance can simply be determined as the ratio of the internally and the externally measured irradiance values.

$$ \tau_{cog,sys} = \frac{E_{int}}{E_{ext}} \qquad (136) $$

$\tau_{cog,sys}$ ... transmittance determined at centre of glazing for entire system

$E_{int}$ ... internally measured irradiance [W/m²]

$E_{ext}$ ... externally measured irradiance [W/m²]

The subscript cog indicates that the measurement has been performed at the centre of glass (further specified in sections 8.1.8.1 and 8.1.8.2), whereas the subscript sys indicates that the measurements consider the full system. The sys subscript distinguishes the quantity from the theoretical quantity as this is based on a calculation that assumes infinitely extended objects. By contrast, the measurements performed here include all effects related to the specimen's finite size (e.g. shading and reflections of the reveal). In order to determine the transmittance of objects causing an inhomogeneous distribution of the irradiance profile, a more sophisticated evaluation is required. To determine the
transmittance of the entire window, many measurement locations can be evenly distributed in the entire internal reveal area, and a statistical evaluation method can be applied (Rüdisser et al., 2021). In order to validate the accuracy of the RadiCal method, it is more efficient to primarily focus on the light scattered by the shading device and the glazing. For this purpose, the just introduced centre-of-glass transmittance is generalized to be applicable for shading devices that cause a vertically periodic irradiance profile. The measurement comprises N sequential measurements taken at vertically distributed equidistant locations in a short period of time. The spacing of the measurements is chosen in a way that M measurements represent exactly one periodic length l, whereas N ≥ M (see Figure 128). The empirically determined transmittance τ<sub>cog,sys,meas</sub> can then be written as:

$$ \tau_{cog,sys,meas} = \frac{1}{N - M + 1} \sum_{i=1}^{N - M + 1} \overline{\overline{\overline{E_{int,i}}}} \overline{\overline{E_{ext,i}}} \quad (137) $$

N ... total number of samples
M ... number of samples covering one periodic length

While the transmittance is now defined as a mean value, the added superscript lines indicate that the used irradiance values are, in turn, based on mean values.

$$ E_{int,meas} = \overline{E_{int,i}} = \frac{1}{M} \sum_{j=i}^{i+M-1} E_{int,j} \qquad \overline{E_{ext,i}} = \frac{1}{M} \sum_{j=i}^{i+M-1} E_{ext,j} \quad (138) $$

$\overline{E_{int,i}}$ ... internal irradiance mean value over one periodic length
$E_{int,j}$ ... internal irradiance measurement reading at location j
$\overline{E_{ext,i}}$ ... external irradiance mean value concurrent with $E_{int,i}$
$E_{ext,j}$ ... external irradiance measurement reading concurrent with $E_{int,j}$

The applied averaging method creates a set of N - M + 1 virtual sub-measurements (see Figure 128). Each of these virtual measurements covers a distance representing the periodic length of the shading device. While these measurements cannot be considered independent because they are based, to varying degrees, on the same measurements, they still enhance the statistical quality of the results, as they cover

<figure><img src="image_134.png" /><figcaption>Figure 128 Averaging method applied to consider periodic length example with N=13 measurements and a periodic length of M=8 samples</figcaption></figure>

slightly deviating spatial and temporal regions. Hence, all parameters relevant for the evaluation of the sub-measurements, i.e. the irradiance values measured by the weather station, the location of the sun and parameters derived from these values, are determined for each sub-measurements specifically and stored along with it.
## 8.1.8. Measurement specifications and model evaluation

### 8.1.8.1. Centre-of-glass approach

A centre-of-glass approach was pursued in the validation measurements, applying to both the empirical and model measurements. This means that the measurements are located in the centre of the unobstructed, visible part of the glazing. The reason for choosing this approach is twofold. First, the measured data at this location allows, to some extent, comparison with results that are determined with the currently applied methods. Commonly used values are calculated based on the assumption of infinitely extended elements, neglecting boundary effects caused by surrounding elements. While the present measurements do consider these effects in both the model and the empirical measurements, their impact is least significant at the centre of the unobstructed glazing. Second, the virtual model measurements require the detection of the irradiance on a sufficiently large detector area. For efficiency reasons, it is not practical to perform model measurements on an area representing the actual size of the pyranometer detection area, as this would require significantly longer computation times. While the different dimensions of the detection areas have a limited impact on most measurements, a significant deviation occurs when the shadow boundary of the direct solar radiation passes the detection area. Therefore, these times, i.e. solar angles, were excluded from the evaluation (see shadow boundary filter in section 8.1.12). In order to optimise the times when this shadow boundary filter is active, the measurement region is located in the centre of the glazing. If the measurement areas is horizontally located near the window frame, the shadow boundary of the opposite side would hit the area at higher incidence angles, making its determination harder.

### 8.1.8.2. Measurement location and spatial resolution

For the reasons just stated, the measurement region is located in the horizontal centre of the visible glazing. Regarding the vertical position, the measurement region is positioned near the geometric centre of the glazing; however, slightly offset below the centre of the visible area by 85 mm (depicted as a red rectangle in Figure 129). The vertical offset takes the mean elevation of the sun into account, as most measurements were performed with relatively high solar elevation angles during the summer period with a south orientation. The unshaded area is, therefore, vertically offset to the bottom side.

The dimensions of the model measurement region are 160 mm x 145 mm. It is essential that the height of the rectangular measurement region in the model equals exactly one vertical slat distance of the blinds or a multiple thereof. For the validation measurements, a height equalling two slat dis-

tances is chosen. The detected irradiance behind the blinds generally exhibits a pronounced periodic spatial variation. Hence, the model measurements have to detect the mean irradiance value for (integer) N slats in the vertical direction, as only this value correctly represents the average irradiance representing the extended system.

<figure><img src="image_135.png" /><figcaption>Figure 129 Location of model- and PyroScanner-measurement regions</figcaption></figure>
The empirical measurement locations for the internal pyranometer, defined in the *PryoScanner* measurement script, are horizontally aligned with the centre of the glazing. Vertically they are centred around the model measurement region. It is not necessary to align the vertical extent of the measurement path to the slat distance, as the evaluation process considers the periodicity of the irradiance. Therefore, the *PryoScanner* was programmed to perform 80 measurements with a vertical spacing of 2 mm equaling a vertical track length of 160 mm (see the black track in Figure 129).

### 8.1.8.3. Temporal resolution and duration

Regarding temporal resolution, the programmed algorithm allows the definition of an idle lead interval before the measurement and an active interval in which the sensor data is recorded. The idle time was set to 400 ms to allow a mechanic and electric stabilization after the movement. The active measurement interval was set to 600 ms. Several hundred independent measurements are captured from the AD converter during this period. The mean value of all measurement readings and their standard deviation is stored in the database for each measurement spot.

The chosen interval parameters are a compromise, as they have to satisfy two conflicting requirements. On the one hand, the exponential response time of the pyranometer demands long measurement periods to allow for the stabilisation of the signal. On the other hand, the comparison against model measurements can be performed best if the conditions (solar angle, direct and diffuse irradiance) do not vary significantly over one measurement cycle. Regarding the typically delayed response of the sensor, it must be noted that even if the pyranometer might take several seconds to adapt from a shaded location to an exposed location, shorter measurement intervals can still be applied. The device's exponential response symmetrically lags on the rising and falling edges of

the signal (see Figure 130). Hence, the integrated signal can still reflect the expected true value as the measurement errors made on the rising and the falling edge compensate. Consequently, the required condition is not that a single interval is sufficiently long but that the measurement period of subsequent samples is long enough to reach the „full“ signal.

On top of the defined measurement intervals, the control and drive time of the *PryoScanner* has to be considered to determine the actual temporal resolution of the measurement. The relevant period can directly be derived from the recorded measurements, as every measurement reading is stored with a timestamp with sub-second precision. The average time interval between two neighbouring measurements is 4.16 seconds. Therefore, the total time required for one measurement cycle comprising 80 individual measurements is 333 seconds.

<figure><img src="image_136.png" /><figcaption>Figure 130 Sampling with exponential response</figcaption></figure>
### 8.1.9. Model used for determination of SIOP

The 3D model used for the raytracing process was generated with the tools FreeCAD (FreeCAD Project association (FPA), 2022) and Blender (The Blender Foundation, 2023). The window frame and glazing were modelled based on original 2D-CAD data provided by the window manufacturer (KATZBECK Fenster GmbH, 2022). The shading devices were modelled using 2D-CAD data provided by the shading manufacturer (Wo&Wo GmbH, 2022). The dimensions of the window reveal, as well as mounting specific parameters, were measured at the site. A rendering of the model performed with the *RadiCal* algorithm and the materials used for the *SIOP* generation is depicted in Figure 116. The relevant material assignments are provided in Table 25. In particular the materials of the glazing panes and shading slats are significant for the internal irradiance measured at the centre of the glazing. The materials of other surfaces have a minor impact. Therefore, generic material definitions were used to model these surfaces. The glazing's internal and external surfaces were modelled based on the specific spectral measurement data applying the inversion approach covered in section 4.12. Specific information regarding the material's spectral data and resulting refractive index functions for the coated and uncoated glazing are covered in the examples contained in section 4.12.3. The material of the polymer-coated slats was

derived using the inversion approach for the subsurface reflectance function covered in section 4.13. Again, the material data of the slats correspond to the examples covered in section 4.13.3. All related information regarding the spectral input data and results of the inversions can be found in this section.

<figure><img src="image_137.png" /><figcaption>Figure 131 3D model (sliced)</figcaption></figure>

<table><caption>Table 16 Material assignments used for validation model</caption><thead><tr><th rowspan="2">Object / surface</th><th rowspan="2">description</th><th colspan="2">Material</th><th rowspan="2">roughness α</th><th rowspan="2">transparent</th></tr><tr><th>refractive index</th><th>subsurface reflectance</th></tr></thead><tbody><tr><td>glass pane 1 ext.<br>glass pane 2 ext.<br>glass pane 2 int.<br>glass pane 3 int.</td><td>float-glass</td><td>determined by inversion<br>based on spectral data<br>(see 4.12.3)</td><td>none</td><td>0</td><td>yes</td></tr><tr><td>glass pane 1 int.<br>glass pane 3 ext.</td><td>low-E coated float<br>float glass</td><td>determined by inversion<br>based on spectral data<br>(see 4.12.3)</td><td>none</td><td>0</td><td>yes</td></tr><tr><td>shading slat<br>white</td><td>polymer coating<br>with pigments<br>white</td><td>generic, constant<br>n̄ = 1.5 – 0i</td><td>determined by inversion<br>based on spectral data<br>(see 4.13.3)</td><td>0</td><td>no</td></tr><tr><td>shading slat<br>anthracite</td><td>polymer coating<br>with pigments<br>anthracite</td><td>generic, constant<br>n̄ = 1.5 – 0i</td><td>determined by inversion<br>based on spectral data<br>(see 4.13.3)</td><td>0</td><td>no</td></tr></tbody></table>
<table><caption>Table 16 Material assignments used for validation model (CONTINUED)</caption><thead><tr><th rowspan="2">Object / surface</th><th rowspan="2">description</th><th colspan="2">Material</th><th rowspan="2">roughness α</th><th rowspan="2">transparent</th></tr><tr><th>refractive index</th><th>subsurface reflectance</th></tr></thead><tbody><tr><td>shading slat silver</td><td>polymer coating with pigments silver</td><td>generic, constant<br>ñ = 1.5 - 0i</td><td>determined by inversion based on spectral data (see 4.13.3)</td><td>0</td><td>no</td></tr><tr><td>window sill</td><td>polymer coating with pigments light grey</td><td>generic, constant<br>ñ = 1.5 - 0i</td><td>generic, constant<br>ρ<sub>ss</sub> = 0.80</td><td>0.02</td><td>no</td></tr><tr><td>window frame</td><td>white paint (on timber)</td><td>generic, constant<br>ñ = 1.6 - 0i</td><td>generic, constant<br>ρ<sub>ss</sub> = 0.95</td><td>0.02</td><td>no</td></tr><tr><td>window spacer</td><td>black paint</td><td>generic, constant<br>ñ = 1.5 - 0i</td><td>generic, constant<br>ρ<sub>ss</sub> = 0.01</td><td>0.02</td><td>no</td></tr><tr><td>pelmet slat tracks</td><td>dark grey/black powder coated</td><td>generic, constant<br>ñ = 1.5 - 0i</td><td>generic, constant<br>ρ<sub>ss</sub> = 0.05</td><td>0.35</td><td>no</td></tr></tbody></table>

### 8.1.10. Evaluation method for model measurement

Following the principles of the *RadiCal* method, the model measurement can be performed by evaluating the corresponding *SIOP* for the measured irradiance profile. The evaluation of the *SIOP* is performed by applying the Perez diffuse sky model (see section 6.2.3).

By applying equation (125) and denoting the *SIOP* that measures the irradiance behind the window at the centre-of-glass position (see Figure 129) as *SIOP*<sub>Ecog</sub> the internal irradiance *E*<sub>int model</sub> behind the window can be determined by:

$$ 
\begin{aligned}
E_{\text{int, model}} = & SIOP_{\text{Ecog}}(\theta, \varphi) \cdot \left( \text{DNI} + \text{DHI} \cdot F_1 \cdot \frac{1}{b} \right) && \text{(direct and circumsolar irr.)} \\
& + \Omega_{\text{sky}} \cdot \bar{S}_{\text{sky}} \cdot \frac{\text{DHI}}{\pi} \cdot (1 - F_1) && \text{(diffuse sky irr.)} \\
& + \Omega_{\text{gnd}} \cdot \bar{S}_{\text{gnd}} \cdot \frac{\text{albedo} \cdot \text{GHI}}{\pi} && \text{(diffuse ground irr.)} \\
& + \Omega_{\text{horizon}} \cdot \bar{S}_{\text{horizon}} \cdot \text{DHI} \cdot F_2 \cdot 4.426 && \text{(diffuse horizon irr.)}
\end{aligned}
\quad (139) $$

The evaluation is an application of the method described in section 6.2.3. The irradiance values *DNI*, *DHI*, *GHI* are the average values for the measurement interval of the corresponding measurement. *F*<sub>1</sub>, *F*<sub>2</sub> and *b* are the Perez parameters. *Ω*<sub>sky</sub>, *Ω*<sub>gnd</sub>, *Ω*<sub>horizon</sub> are the relevant visible solid angles. For the vertical orientation of the windows, the solid angles are *Ω*<sub>sky</sub> = *Ω*<sub>gnd</sub> = *πsr* and *Ω*<sub>horizon</sub> = 0.3556 *sr*. The *SIOP* average values *S̄*<sub>sky</sub>, *S̄*<sub>gnd</sub> and *S̄*<sub>horizon</sub> are precalculated, whereas the direct radiation term (*SIOP*<sub>Ecog</sub>(*θ*, *φ*)) is evaluated for every measurements using the average values for solar angles *θ* and *φ* during the corresponding measurements.

In analogy to equation (136), the modelled internal irradiance is related to the external radiance to determine the model transmittance *τ*<sub>cog,sys-model</sub>, measured at the center of glazing (*cog*) and considering the entire window system (*sys*):

$$ \tau_{\text{cog,sys,mode1}} = \frac{E_{\text{int,mode1}}}{E_{\text{ext,mode1}}} \quad (140) $$

The external irradiance *E*<sub>ext,mode1</sub> is the direct and diffuse irradiance on the window/façade plane calculated by applying the Perez model for the diffuse component, equation (119).
### 8.1.11. Albedo values

The albedo value characterising the reflectance of the environment remains the only undetermined quantity of the evaluation process based on equation (139). The evaluation can therefore be used to determine the albedo values empirically. It has been found that the sensitivity of the internally measured irradiance ranges from 0.1% to almost 1% per 1% change of albedo, depending on the test case. Since it would contradict the objective of the validation to individually optimise the albedo value for each case, three different constant average values were assigned to the 23 cases. It is important to note that the validation measurements and evaluations demonstrated well that the constant albedo approach is linked to significant inaccuracies, as already discussed in section 6.2.4. In particular, the bright roof located in the immediate south neighbourhood of the testing facility and tilted towards it (see Figure 120 top) affected the results of the measurements by directing significant amounts of direct solar radiation towards the target. In contrast, the dark membrane on the lab roof (see Figure 120 bottom-right) exhibited a relatively low albedo value. Its albedo value was determined in another project to be 8%. In order to accurately model such an environment, an approach as proposed in section 6.2.4 has to be followed. Since the approach is currently not fully implemented and should be the subject of a future separate study, a crude model based on the same idea is applied: For all south-facing cases with a significantly “downward-looking” SIOP the albedo value of 8% is used. This reflects the fact that the tilted slats effectively filter radiation originating from the floor (lab roof). The downward-looking criterion is defined to be relevant for all shaded cases with at least 45 degrees tilt angle. For all other cases, the opposite roof's reflection must be considered; therefore, a value of 15% is applied for these cases.

Finally, an albedo value of 11% is applied for the single, unshaded west-facing case since meadows, trees and shrubbery dominate the environment in this direction.

The three-value approach considers the conflicting requirements of the evaluation. The constant albedo values were derived from the overall measurement results; however, they are not optimised for each of the 23 cases, which would lead to better validation results. All determined albedo values lie within a plausible range.

### 8.1.12. Comparison of model measurements and empirical measurements

#### 8.1.12.1. Definition of error measures

In order to assess the quality of the results achieved in the model measurements, the following parameters are derived and provided in tabular form for each test case:

**1. Mean error of transmittance**

$$ME_{\tau} = \frac{1}{N} \sum_{i=1}^{N} \Delta\tau_{i} \quad [1] \qquad (141)$$

The mean error indicates the bias (or offset) of the model reconstruction vs. the measured transmittance. This quantity is considered the most important measure. It is dimensionless and effectively relates the error of internal irradiance to external irradiance. Using relative errors is less appropriate, as these increase inversely proportionally to the transmittance values.

**2. Mean absolute error of transmittances**

$$\Delta\tau_{i} = \tau_{cog,sysmodel,i} - \tau_{cog,sysmeas,i} \qquad (142)$$
$$MAE_{\tau} = \frac{1}{N} \sum_{i=1}^{N} |\Delta\tau_{i}| \quad [1]$$

The mean absolute error of transmittance measures the average error of the model prediction vs. the measured values. Unlike the (signed) mean error, it also considers the statistical deviations resulting primarily from variations in the empirical values.

### 3. Standard deviation of error of transmittance

$$SE_{\tau} = \sqrt{\frac{1}{N} \sum_{i=1}^{N} (\Delta\tau_{i} - ME_{\tau})^{2}} \quad [1] \quad (143)$$

The quantity is a measure of the variation of the observed errors around the mean error of transmittance.

### 4. Mean absolute error of internal irradiance

$$MAE_{E} = \frac{1}{N} \sum_{i=1}^{N} |E_{int, model, i} - E_{int, meas, i}| \quad [W/m^{2}] \quad (144)$$

The mean error indicates the average absolute error of the modelled internal irradiance in comparison to the measured value. It considers statistical variation and a systematic offset.

### 5. Linear regression through the origin of internal irradiance values

A linear regression through the origin can be applied since both the calibrated and model measurements converge towards zero as the solar irradiance diminishes. In the evaluations, the slope of the fit line, as well as the coefficient of determination R² are provided.

### 6. Linear regression of internal irradiance values

Additionally, a standard linear regression is applied. In the evaluations, the slope and the intercept of the fit line, as well as the coefficient of determination R² are provided.

Note: It is essential to point out that the regression analysis on the irradiance measurements (5. and 6.) only allows relevant conclusions if sufficient data points are available and if their values are spread over a sufficiently large range. Only some of the test cases fulfil this requirement.

#### 8.1.12.2. Treatment of unshaded cases

In order to apply the same evaluation process as for the shaded cases (section 8.1.7), the internal and external irradiance values, as well as the weather station data, were averaged for 1-minute time periods. Each of these averaged values was subsequently evaluated as an individual measurement.

#### 8.1.12.3. Filtering

Before the comparison is performed, four different filters are applied to the data in order to improve the quality of the comparison. The function and purpose of each filter are described below.

#### 8.1.12.4. Irradiance consistency filter

The model measurement results are generated based on a SIOP evaluation using the Perez irradiance model (see section 6.2.3). For this reason, it is necessary to provide the irradiance components DNI, DHI and GHI. For the validation measurements, these values are provided by the weather station of the laboratory, based on the calibrated and highly accurate sun tracker device and the associated pyranometers. Apart from using the irradiance components for the SIOP evaluation, they are also used to determine the total on-plane irradiance for the incidence plane parallel to the glazing. The standard formula to determine the irradiance on a vertical plane is used for this purpose.
$$E_{on-plane} = \frac{DHI_{Perez}}{2} + \frac{\alpha GHI}{2} + \max(\vec{n} \cdot \vec{s}, 0) \cdot DNI \quad (145)$$

DHI<sub>Perez</sub> ... diffuse radiation according to equation (119)

α ... albedo

n̂ ... normal vector of glazing

ŝ ... vector pointing to direction of sun

As the external pyranometer of the *PyroScanner* directly measures this on-plane irradiance, the above equation was used to derive the initially unknown albedo value. By minimizing the mean errors of the calculation vs. the measurement, the average albedo value for the south-oriented un-shaded case was determined as $α = 0.15$. Two additional values were applied in the evaluations (see section 8.1.11).

Beyond the determination of the albedo, the comparison of the calculated vs. the directly measured on-plane irradiance is performed for every measurement. Significant deviations between the two irradiance values can arise. They can be attributed to two different causes:

1. Reflection spots from highly reflecting surfaces on the laboratory roof or the surrounding environment can hit only one of the two sensors or both with varying irradiance.

2. The Perez diffuse sky model only inadequately describes the current irradiance profile (e.g. due to pronounced broken cloud patterns).

As both cases do not allow a meaningful comparison of the Perez-based *SIOP* evaluation vs. the measurement, variations of the modelled to the measured on-plane irradiance are only accepted up to a specified threshold. In order to not exclude too many measurements from the evaluation, a threshold of 10% is defined. The irradiance consistency filter can, therefore, be stated as:

$$\frac{|E_{on-plane} - E_{ext}|}{0.5 \cdot E_{on-plane} + 0.5 \cdot E_{ext}} < 0.1 \quad (146)$$

$E_{on-plane}$ ... on plane irradiance, determined with Perez model

$E_{ext}$ ... irradiance measured at external pyranometer

It must be noted that the application of this filter is not a shortcoming of the model but a requirement that results from the chosen validation method based on real conditions. Since the *SIOPs* are evaluated using the Perez irradiance model, it has to be ensured that it adequately describes the current sky condition, as only these conditions allow meaningful validation. Alternative sky models or lab lighting conditions would be associated with different irradiance profiles that could as well be evaluated using the *SIOP* concept.

### 8.1.12.5. Minimum irradiance filter

Only measurements with global irradiance on the incidence plane greater than $10 W/m^2$ were considered in the evaluation. The reason for that is that even pyranometers of class A typically show errors in the range of a few $W/m^2$ that are related to environmental temperature change, thermal radiation and electronics. The relevant standard (ISO 9060, 2018) states an acceptable “total zero-offset of $\pm 10 W/m^2$” for class A devices.
### 8.1.12.6. Shadow boundary filter

As stated above, the slightly different definitions of measurement regions in the model and empirical measurement lead to systematic variations of the determined irradiance values in times when the shadow boundary of the direct radiation passes over them. It is important to point out that the shadow boundary filter does not have to be applied because one of the measurements provides faulty results but because the results of the model and empirical measurements do not correspond during these periods.

Based on geometric considerations, the shadow boundary filter was defined based on the incidence angle. In vector form, the filter can be written as:

$$0.104 < \vec{n} \cdot \vec{s} < 0.387 \qquad (147)$$
$$\begin{array}{l}
\vec{n} \text{ ... normal vector of glazing} \\
\vec{s} \text{ ... vector pointing to direction of sun}
\end{array}$$

This effectively means that measurements are excluded when the incidence angle of direct sunlight ranges between 67.2° and 84.0°. Hence, all measurements with an incidence angle in the range of 0° to 67.2° (potential direct exposure) and 84.0° to 180° (diffuse exposure only) pass the filter.

### 8.1.12.7. Outlier filter

Outlier removal is only applied very cautiously, as significant outlier removal could cause an increase in the correlation measures. Outliers can, e.g., be caused by specular reflections on the internal or external pyranometer, shadings by staff working on the testing facility or electromagnetic interference. A 3-sigma approach was applied to filter such outliers:

$$|\Delta E_{int}| < 3 \cdot \sigma_{|\Delta E_{int}|}$$
$$\Delta E_{int} = E_{int, model} - E_{int, meas} \quad \text{and} \quad \sigma_{|\Delta E_{int}|} = \sqrt{\frac{1}{N-1} (\Delta E_{int})^2} \qquad (148)$$

In general, less than 2% of the data is effectively filtered based on this criterion. Hence, the impact on the correlation and error analysis is considered minimal.
### 8.1.13. Test cases and specimen specification

Figure 132 shows the 3D models used for the slats. The models as well as the actual shading devices were provided by WO&WO Sonnenlichtdesign GmbH & Co KG (Wo&Wo GmbH, 2022). The model K80 has the traditional C-shape, with folded edges on either side. The slat width is 80.0 mm, the total height is 11.7 mm, and the folded edges' external radius is 4.3 mm. The model VZ90 is a Z-shaped profile with folded edges on either side. The slat width is 93.5 mm, the total height is 10.4 mm and the external folded edge radius is 3.0 mm. Both slats are modelled as 3D solids. The thickness of the coated aluminium profile of both slats is 0.44 mm.

<figure><img src="image_138.png" /><figcaption>Figure 132 3D models demonstrating the applied level of detail of slat types used</figcaption></figure>

The wireframe models of the slats as depicted in Figure 132 are sufficiently detailed as they include information on the orientation of the surface normals. This information is used in the smooth-shading feature (see sections 7.3 and 0) to reconstruct the slat's reflectance smoothly.

Three different coatings were tested for the K80 slat type and one coating for the VZ90 slat type. Figure 133 provides an overview of the four slat specimens that were used in the validation measurements.

<figure><img src="image_139.png" /><figcaption>Figure 133 Overview of the four slat specimens used</figcaption></figure>

Several different slat positions were analysed in the tests. In order to gather more data for the validations, the tests were performed with and without the influence of the glazing. To perform measurements without glazing, the window sash was opened at an angle of 90° degrees into the room. The cases that could be evaluated do not cover the potential parameter range in a regular way. Some cases had to be discarded due to the unavailability of the testing facility. In other cases, the measurement readings were erroneous or missing, either from the weather station or the internal measurements. Many of the tests were controlled and monitored remotely, and all tests had been carried out before the evaluation method was implemented. Hence, it was not possible to repeat unusable measurements in the course of the project. Still, 1061 single measurements for 21 separate shadowed cases were evaluated. Additionally, the shading was fully retracted, and measurements were performed with a stationary internal pyranometer located at the centre of the glazing. The measurements were averaged over 1-minute periods to apply the same evaluation as for the shaded cases.
These measurements were performed for two different orientations of the window (south and west facing). Based on this approach another 9438 measurement points were evaluated for the two unshaded cases. An overview of all validation test cases is provided in Table 17.

<table><caption>Table 17 Overview of test cases</caption><thead><tr><th>case no</th><th>case name</th><th>orient.</th><th>window state</th><th>slat type</th><th>slat angle</th><th>meas. count</th></tr></thead><tbody><tr><td>1</td><td>UNSHADED_SOUTH</td><td>south</td><td>win. closed</td><td>-</td><td>-</td><td>3347</td></tr><tr><td>2</td><td>UNSHADED_WEST</td><td>west</td><td>win. closed</td><td>-</td><td>-</td><td>6091</td></tr><tr><td>3</td><td>ANTHRA_wC_00</td><td>south</td><td>win. closed</td><td>anthracite, type K80</td><td>0</td><td>206</td></tr><tr><td>4</td><td>ANTHRA_wC_45</td><td>south</td><td>win. closed</td><td>anthracite, type K80</td><td>45</td><td>58</td></tr><tr><td>5</td><td>ANTHRA_wC_m45</td><td>south</td><td>win. closed</td><td>anthracite, type K80</td><td>-45</td><td>42</td></tr><tr><td>6</td><td>SILVER_wC_00</td><td>south</td><td>win. closed</td><td>silver, type K80</td><td>0</td><td>78</td></tr><tr><td>7</td><td>SILVER_wC_m45</td><td>south</td><td>win. closed</td><td>silver, type K80</td><td>-45</td><td>64</td></tr><tr><td>8</td><td>WHITE_wC_00</td><td>south</td><td>win. closed</td><td>white, type K80</td><td>0</td><td>136</td></tr><tr><td>9</td><td>WHITE_wC_53</td><td>south</td><td>win. closed</td><td>white, type K80</td><td>53</td><td>22</td></tr><tr><td>10</td><td>Z-SILVER_wC_00</td><td>south</td><td>win. closed</td><td>silver, type VZ90</td><td>0</td><td>6</td></tr><tr><td>11</td><td>Z-SILVER_wC_58</td><td>south</td><td>win. closed</td><td>silver, type VZ90</td><td>58</td><td>40</td></tr><tr><td>12</td><td>ANTHRA_wO_00</td><td>south</td><td>win. open</td><td>anthracite, type K80</td><td>0</td><td>123</td></tr><tr><td>13</td><td>ANTHRA_wO_45</td><td>south</td><td>win. open</td><td>anthracite, type K80</td><td>45</td><td>50</td></tr><tr><td>14</td><td>ANTHRA_wO_m45</td><td>south</td><td>win. open</td><td>anthracite, type K80</td><td>-45</td><td>44</td></tr><tr><td>15</td><td>SILVER_wO_45</td><td>south</td><td>win. open</td><td>silver, type K80</td><td>0</td><td>27</td></tr><tr><td>16</td><td>SILVER_wO_m21</td><td>south</td><td>win. open</td><td>silver, type K80</td><td>-21</td><td>23</td></tr><tr><td>17</td><td>SILVER_wO_m45</td><td>south</td><td>win. open</td><td>silver, type K80</td><td>-45</td><td>27</td></tr><tr><td>18</td><td>WHITE_wO_00</td><td>south</td><td>win. open</td><td>white, type K80</td><td>0</td><td>37</td></tr><tr><td>19</td><td>WHITE_wO_45</td><td>south</td><td>win. open</td><td>white, type K80</td><td>45</td><td>24</td></tr><tr><td>20</td><td>WHITE_wO_m45</td><td>south</td><td>win. open</td><td>white, type K80</td><td>-45</td><td>24</td></tr><tr><td>21</td><td>Z-SILVER_wO_00</td><td>south</td><td>win. open</td><td>silver, type VZ90</td><td>0</td><td>15</td></tr><tr><td>22</td><td>Z-SILVER_wO_45</td><td>south</td><td>win. open</td><td>silver, type VZ90</td><td>45</td><td>2</td></tr><tr><td>23</td><td>Z-SILVER_wO_58</td><td>south</td><td>win. open</td><td>silver, type VZ90</td><td>58</td><td>13</td></tr></tbody></table>
## 8.2. Results

### 8.2.1. Introduction

In the following sections, only the two unshaded cases and the four cases with horizontal slats (0°) of the total 23 test cases are presented in the main text. These cases were selected because they have the highest number of evaluated measurements. Hence, the results provide the highest statistical relevance. However, the results of all test cases are provided in the same way in Appendix A.

All evaluation results, as specified in section 8.1.12, are presented in tabular form for each test case. In addition to the evaluation results, relevant specifications of the measurement conditions are provided in a statistical form. The resulting irradiances and transmittances depend on the actual irradiance profile and therefore cover a specific range. Since these parameters' mean or average values are of little relevance, the quantities minimum, first quantile, third quantile and maximum are provided only. Finally, two charts offer a graphical overview of the main results for each test case: a correlation chart shows the modelled irradiance values as a function of the measured values; a second chart depicts the determined transmittance values as a function of the solar azimuth angle. In both charts, colours are used to indicate the sky condition at the time of the measurement. The diffuse ratio that was used for this purpose relates the diffuse irradiance on the horizontal plane (DHI) to the total irradiance (GHI).

Additionally, the measured and modelled irradiance values are presented as time series for the two unshaded cases (see Figure 136 and Figure 139). Presenting the shaded cases as time series is less relevant, as the duration of a single measurement is longer and the temporal resolution is significantly more coarse.

### 8.2.2. Case #1 – UNSHADED_SOUTH

<figure><img src="image_140.png" /><figcaption>Figure 134 Case #1 – int. irradiance correlation
modelled values as a function of measured values</figcaption></figure>

<figure><img src="image_141.png" /><figcaption>Figure 135 Case #1 – modelled transmittances
as a function of solar azimuth angle</figcaption></figure>
<table><caption>Table 18 Validation case #1 – results and specifications</caption><thead><tr><th colspan="5">Case specifications</th></tr></thead><tbody><tr><td><strong>case #1</strong><br>case name:<br>UNSHADED_SOUTH</td><td>case description:</td><td>win.state: window closed<br>slats: -<br>slat-angle: -°</td><td></td><td></td></tr><tr><td colspan="5"><strong>Transmittance τ<sub>COG,SYS</sub></strong></td></tr><tr><td colspan="3">range of determined values (min | q.1 | q.3 | max) [1]:</td><td>errors:</td><td>MAE<sub>τ</sub>: 0.024<br>ME<sub>τ</sub>: 0.013<br>SE<sub>τ</sub>: 0.040</td></tr><tr><td colspan="5"><strong>Internal irradiance E<sub>int-model</sub></strong></td></tr><tr><td colspan="3">range of determined values (min | q.1 | q.3 | max) [W/m<sup>2</sup>]:</td><td>errors:</td><td>MAE<sub>E</sub>: 1.8 W/m<sup>2</sup></td></tr><tr><td colspan="5"><strong>Internal irradiance correlation model vs. measurement:</strong></td></tr><tr><td colspan="3">number of points: 3347</td><td></td><td></td></tr><tr><td>lin. regression through origin</td><td colspan="2">y=0.9785·x</td><td>R<sup>2</sup>: 0.999</td><td></td></tr><tr><td>lin. regression</td><td colspan="2">y=0.9738·x + 0.4422</td><td>R<sup>2</sup>: 0.998</td><td></td></tr><tr><td colspan="5"><strong>Measurement specifications (min | q.1 | q.3 | max):</strong></td></tr><tr><td colspan="3">sun azimuth range [°]: 56.7 | 84.7 | 279.6 | 302.3</td><td></td><td></td></tr><tr><td colspan="3">sun elevation range [°]: 0.6 | 15.3 | 37.3 | 65.8</td><td></td><td></td></tr><tr><td colspan="3">incidence angle range [°]: 64.9 | 84.8 | 106.2 | 123.3</td><td></td><td></td></tr><tr><td colspan="3">diffuse ratio (DHI/GHI) range [1]: 0.08 | 0.21 | 0.91 | 1.00</td><td></td><td></td></tr><tr><td colspan="3">meas. period: 2021.07.02 | 2021.07.05 | 2021.07.10 | 2021.07.12</td><td></td><td></td></tr></tbody></table>

<figure><img src="image_142.png" /><figcaption>Figure 136 Case #1 – modelled and measured internal irradiance values as time series exemplary data for five days (unshaded; south)</figcaption></figure>

The grey regions indicate time periods excluded from the evaluation due to the shadow boundary filter criterion (see section 8.1.12).
### 8.2.3. Case #2 – UNSHADED_WEST

<figure><img src="image_143.png" /><figcaption>Figure 137 Case #2 – int. irradiance correlation modelled values as a function of measured values</figcaption></figure>

<figure><img src="image_144.png" /><figcaption>Figure 138 Case #2 – modelled transmittances as a function of solar azimuth angle</figcaption></figure>

<table><caption>Table 19 Validation case #2 – results and specifications</caption><thead><tr><th colspan="5">Case specifications</th></tr></thead><tbody><tr><td>case #2<br>case name:<br>UNSHADED_WEST</td><td>case description:</td><td>win.state: window closed<br>slats: 0<br>slat-angle: -°</td><td></td><td></td></tr><tr><td colspan="5">Transmittance τ<sub>cog,sys</sub></td></tr><tr><td>range of determined values (min | q.1 | q.3 | max) [1]:</td><td></td><td>errors: MAEτ: 0.030<br>MEτ: 0.003<br>SEτ: 0.040</td><td></td><td></td></tr><tr><td>range of determined values (min | q.1 | q.3 | max) [W/m²]:</td><td></td><td>errors: MAE<sub>E</sub>: 6.1 W/m²</td><td></td><td></td></tr><tr><td colspan="5">Internal irradiance E<sub>int model</sub></td></tr><tr><td>number of points: 6091<br>lin. regression through origin y=1.0105·x<br>lin. regression y=1.0119·x - 0.3049</td><td></td><td>R²: 0.996<br>R²: 0.992</td><td></td><td></td></tr><tr><td colspan="5">Internal irradiance correlation model vs. measurement:</td></tr><tr><td>sun azimuth range [°]: 59.4 | 104.8 | 258.8 | 303.1<br>sun elevation range [°]: -1.0 | 22.4 | 52.1 | 64.6<br>incidence angle range [°]: 16.8 | 38.8 | 138.4 | 162.8<br>diffuse ratio (DHI/GHI) range [1]: 0.12 | 0.36 | 1.00 | 1.00<br>meas. period: 2021.07.12 | 2021.07.17 | 2021.08.01 | 2021.08.05</td><td></td><td></td><td></td><td></td></tr></tbody></table>
<figure><img src="image_145.png" /></figure>

Figure 139 Case #2 – modelled and measured internal irradiance values as time series exemplary data for five days (unshaded; west)
### 8.2.4. Case #3 – ANTHRA_wC_00

<figure><img src="image_146.png" /><figcaption>Figure 140 Case #3 – int. irradiance correlation modelled values as a function of measured values</figcaption></figure>

<figure><img src="image_147.png" /><figcaption>Figure 141 Case #3 – modelled transmittances as a function of solar azimuth angle</figcaption></figure>

<table><caption>Table 20 Validation case #3 – results and specifications</caption><thead><tr><th colspan="5">Case specifications</th></tr></thead><tbody><tr><td>case #3</td><td>case description:</td><td>win.state:</td><td>window closed</td><td></td></tr><tr><td>case name:</td><td></td><td>slats:</td><td>anthracite, type K80</td><td></td></tr><tr><td>ANTHRA_wC_00</td><td></td><td>slat-angle:</td><td>0°</td><td></td></tr><tr><td colspan="5">Transmittance τ<sub>cog,sys</sub></td></tr><tr><td>range of determined values (min | q.1 | q.3 | max) [1]:</td><td></td><td>errors:</td><td>MAEτ: 0.010</td><td></td></tr><tr><td>0.029 | 0.041 | 0.120 | 0.207</td><td></td><td>MEτ: 0.003</td><td></td><td></td></tr><tr><td></td><td></td><td>SEτ: 0,016</td><td></td><td></td></tr><tr><td colspan="5">Internal irradiance E<sub>int model</sub></td></tr><tr><td>range of determined values (min | q.1 | q.3 | max) [W/m²]:</td><td></td><td>errors:</td><td>MAE<sub>E</sub>: 2.1 W/m²</td><td></td></tr><tr><td>2.8 | 17.3 | 25.7 | 39.3</td><td></td><td></td><td></td><td></td></tr><tr><td colspan="5">Internal irradiance correlation model vs. measurement:</td></tr><tr><td>number of points:</td><td>206</td><td></td><td></td><td></td></tr><tr><td>lin. regression through origin</td><td>y=0.9741·x</td><td>R²:</td><td>0.985</td><td></td></tr><tr><td>lin. regression</td><td>y=0.9184·x + 1.3131</td><td>R²:</td><td>0.870</td><td></td></tr><tr><td colspan="5">Measurement specifications (min | q.1 | q.3 | max):</td></tr><tr><td>sun azimuth range [°]:</td><td>85.3 | 153.4 | 223.5 | 289.8</td><td></td><td></td><td></td></tr><tr><td>sun elevation range [°]:</td><td>3.6 | 45.8 | 58.3 | 61.5</td><td></td><td></td><td></td></tr><tr><td>incidence angle range [°]:</td><td>53.4 | 60.5 | 65.6 | 109.8</td><td></td><td></td><td></td></tr><tr><td>diffuse ratio (DHI/GHI) range [1]:</td><td>0.11 | 0.17 | 0.96 | 1.00</td><td></td><td></td><td></td></tr><tr><td>meas. period:</td><td>2020.07.27 | 2020.08.01 | 2020.08.23 | 2020.08.25</td><td></td><td></td><td></td></tr></tbody></table>
### 8.2.5. Case #6 – SILVER_wC_00

<figure><img src="image_148.png" /><figcaption>Figure 142 Case #6 – int. irradiance correlation modelled values as a function of measured values</figcaption></figure>

<figure><img src="image_149.png" /><figcaption>Figure 143 Case #6 – modelled transmittances as a function of solar azimuth angle</figcaption></figure>

<table><caption>Table 21 Validation case #6 – results and specifications</caption><thead><tr><th colspan="4">Case specifications</th></tr></thead><tbody><tr><td><strong>case #6</strong><br>case name:<br>SILVER_wC_00</td><td>case description:</td><td>win.state: window closed<br>slats: silver, type K80<br>slat-angle: 0°</td><td></td></tr><tr><td><strong>Transmittance τ<sub>cog,sys</sub></strong><br>range of determined values (min | q.1 | q.3 | max) [1]:<br>0.067 | 0.084 | 0.230 | 0.237</td><td>errors:</td><td>MAEτ: 0.011<br>MEτ: -0.006<br>SEτ: 0,014</td><td></td></tr><tr><td><strong>Internal irradiance E<sub>int model</sub></strong><br>range of determined values (min | q.1 | q.3 | max) [W/m²]:<br>4.3 | 15.1 | 35.6 | 65.4</td><td>errors:</td><td>MAE<sub>E</sub>: 1.9 W/m²</td><td></td></tr><tr><td colspan="4"><strong>Internal irradiance correlation model vs. measurement:</strong></td></tr><tr><td>number of points: 78</td><td></td><td></td><td></td></tr><tr><td>lin. regression through origin y=1.0433·x</td><td></td><td>R²: 0.996</td><td></td></tr><tr><td>lin. regression y=1.0123·x + 1.0775</td><td></td><td>R²: 0.981</td><td></td></tr><tr><td colspan="4"><strong>Measurement specifications (min | q.1 | q.3 | max):</strong></td></tr><tr><td>sun azimuth range [°]:</td><td>86.5 | 181.6 | 272.9 | 298.6</td><td></td><td></td></tr><tr><td>sun elevation range [°]:</td><td>4.9 | 28.5 | 63.4 | 65.7</td><td></td><td></td></tr><tr><td>incidence angle range [°]:</td><td>62.4 | 65.7 | 93.1 | 118.5</td><td></td><td></td></tr><tr><td>diffuse ratio (DHI/GHI) range [1]:</td><td>0.09 | 0.14 | 0.70 | 1.00</td><td></td><td></td></tr><tr><td>meas. period:</td><td>2020.07.03 | 2020.07.04 | 2020.07.06 | 2020.07.27</td><td></td><td></td></tr></tbody></table>