# 4.7. Generalized Fresnel coefficients

## 4.7.1. Background and theory

The derivations and equations presented below are provided in different forms in various publications. A comprehensive derivation that starts with the Maxwell equations in vector form can be found in the book *Thin-film Optical filters* (Macleod, 2001). For the practical implementation, the publications by Craig (1987) and Keller (2001) were helpful.

The commonly known Fresnel's equations are used as a starting point for reference. They can be derived by solving Maxwell's equations for non-absorbing media. They are commonly stated as:

$$r_s = \frac{n_1 \cos\theta_1 - n_2 \cos\theta_2}{n_1 \cos\theta_1 + n_2 \cos\theta_2} \qquad t_s = \frac{2 n_1 \cos\theta_1}{n_1 \cos\theta_1 + n_2 \cos\theta_2} \qquad (41)$$

$$r_p = \frac{n_2 \cos\theta_1 - n_1 \cos\theta_2}{n_2 \cos\theta_1 + n_1 \cos\theta_2} \qquad t_p = \frac{2 n_1 \cos\theta_1}{n_2 \cos\theta_1 + n_1 \cos\theta_2} \\ n_1 \text{ ... refractive medium 1 (incidence)} \qquad \theta_1 \text{ ... angle of incidence [rad]} \\ n_2 \text{ ... refractive medium 2} \qquad \theta_2 \text{ ... angle of refraction [rad]}$$

The ratio of the incidence angle $\theta_1$ and transmittance angle $\theta_2$ is determined by Snell's law, see eq. (18). The Fresnel coefficients $r_{s,p}$ and $t_{s,p}$ determine the amplitudes of the reflected and transmitted radiation for each polarisation state. The coefficients for the reflected and transmitted power, referred to as reflectance $R_{s,p}$ and transmittance $T_{s,p}$ are given by:

$$R_{s,p} = |r_{s,p}|^2 \qquad T_{s,p} = \frac{n_2 \cos\theta_2}{n_1 \cos\theta_1} |t_{s,p}|^2 \qquad (42)$$

This takes into account that the power is proportional to the square of the amplitude and that the angle as well as the impedance of the medium will change for the transmitted beam.

For normal incidence ($\theta_1 = \theta_2 = 0$) the amplitude coefficients $(r, t)$ and power coefficients $(R, T)$ simplify to the elegant form:

$$r = \frac{n_2 - n_1}{n_1 + n_2} \qquad t = \frac{2 n_1}{n_1 + n_2} \qquad (43)$$

$$R = \left( \frac{n_2 - n_1}{n_1 + n_2} \right)^2 \qquad T = \frac{n_2}{n_1} \left( \frac{2 n_1}{n_1 + n_2} \right)^2$$

By using complex-valued refraction indices and introducing the concept of effective refractive indices, more general expressions valid for absorbing media and oblique (=non-normal) incidence angles will be derived now. The finally obtained relations show a striking similarity to equations (43).

The derivation is carried out by solving Maxwell's equation for interfaces of materials and demanding continuity of the electric and magnetic field at the boundary. While this can be done by applying vector calculus in three dimensions, a more straightforward and more efficient approach is possible when the energy flow normal to the interface is considered. In this case, the relevant electric field vector $\vec{E}$ and the magnetic field vector $\vec{H}$ are oriented parallel (tangential) to the boundary. It is, therefore, sufficient to consider the amplitudes only, denoted as $\vec{E}$ and $\vec{H}$. In order to include phase information, the amplitudes are represented by complex-valued numbers. The tilde in superscript is further used to indicate all quantities that are complex numbers. $s$ and $p$ subscripts are used in the following to indicate the direction of polarisation for which the relation is valid. If the letters for
both directions are indicated in an equation (separated by a comma), it will mean that the equation is separately valid for either polarisation state (using only *s* or only *p* polarised quantities throughout the equation).

Following the approach of considering the energy flows perpendicular to the boundary, it is helpful to introduce an *effective index of refraction* $\tilde{\eta}$, which links the two tangential electromagnetic field components $\tilde{E}$ and $\tilde{H}$:

$$ \tilde{\eta}_{s,p} = \frac{\tilde{H}_{s,p}}{\tilde{E}_{s,p}} \qquad (44) $$

In the applied approach considering the flow of energy perpendicular to the boundary, it can be shown that for the s-polarised and p-polarised components, the effective indices of refraction can be derived of the complex-valued index of refraction $\tilde{n}$ and the propagation angle $\tilde{\theta}$ by:

$$ \tilde{\eta}_s = \tilde{n} \cos(\tilde{\theta}) \qquad \tilde{\eta}_p = \frac{\tilde{n}}{\cos(\tilde{\theta})} \qquad (45) $$

It is important to note that the general propagation angle $\tilde{\theta}$ is now a complex quantity. Consequently, the angle of refraction $\tilde{\theta}_2$ for the transmitted beam has to be derived from the angle of incidence $\tilde{\theta}_1$ by solving Snell's law in the complex-valued form:

$$ \sin \tilde{\theta}_1 \tilde{n}_1 = \sin \tilde{\theta}_2 \tilde{n}_2 \qquad (46) $$

Unlike its real-valued counterpart, the *complex-valued effective index of refraction* $\tilde{\eta}$ contains angular information, as well as information on the attenuation (through the complex-valued index of refraction $\tilde{n}$). The amplitude transmission coefficients for reflection ($\tilde{r}$) and transmission ($\tilde{t}$) and the resulting reflectance ($R$) and transmittance ($T$), at the boundary of two absorbing or non-absorbing media 1 and 2, can for any incidence angle be stated as:

$$ \tilde{r}_{1 \to 2}^{s,p} = \frac{\tilde{\eta}_1^{s,p} - \tilde{\eta}_2^{s,p}}{\tilde{\eta}_1^{s,p} + \tilde{\eta}_2^{s,p}} \qquad \tilde{t}_{1 \to 2}^{s,p} = \frac{2 \tilde{\eta}_1^{s,p}}{\tilde{\eta}_1^{s,p} + \tilde{\eta}_2^{s,p}} \qquad (47) $$

$$ R_{1 \to 2}^{s,p} = |\tilde{r}_{1 \to 2}^{s,p}|^2 \qquad T_{1 \to 2}^{s,p} = \frac{\tilde{\eta}_{s,p} \tilde{E}_1^{s,p}}{\tilde{\eta}_{s,p} \tilde{E}_1^{s,p} + \tilde{\eta}_{s,p} \tilde{H}_1^{s,p}} \qquad (48) $$

It can be seen how this very elegant and simple form of the generalized form of Fresnel equations results only through the introduction of the complex-valued effective index of refraction. The similarity to the earlier, non-complex-valued form is remarkable, see equation (43).

However, taking the generalization a step further, a different definition for the amplitude coefficients is chosen to additionally consider the effect of potentially present thin-films in between the two media. It can be shown (Macleod, 2001) that they can be calculated by the tangential electric and magnetic field amplitudes in the incident medium:

$$ \tilde{r}_{1 \to 2}^{s,p} = \frac{\tilde{\eta}_1^{s,p} \tilde{E}_1^{s,p} - \tilde{H}_1^{s,p}}{\tilde{\eta}_1^{s,p} \tilde{E}_1^{s,p} + \tilde{H}_1^{s,p}} \qquad \tilde{t}_{1 \to 2}^{s,p} = \frac{2 \tilde{\eta}_1^{s,p}}{\tilde{\eta}_1^{s,p} \tilde{E}_1^{s,p} + \tilde{H}_1^{s,p}} \qquad (48) $$

If the incident field amplitudes $\tilde{E}_1^{s,p}$ and $\tilde{H}_1^{s,p}$ are determined by:

$$ \begin{pmatrix} \tilde{E}_1^{s,p} \\ \tilde{H}_1^{s,p} \end{pmatrix} = \bar{C}^{s,p} \begin{pmatrix} 1 \\ \tilde{\eta}_2^{s,p} \end{pmatrix} \qquad (49) $$

Where $\bar{C}^{s,p}$ is the *characteristic matrix* of a thin-film assembly consisting of N thin-films.
In the absence of thin-films, the matrix $\bar{\mathcal{C}}^{s,p}$ can be omitted, i.e. replaced with the unity matrix. For this case, it is evident that equations (48) will turn into equations (47) again. Hence, the definition for the reflection and transmission coefficients provided by (48) is the most general, as it allows consideration of 0 to N thin-films at the boundary of two materials. For this reason, this definition has been chosen for the actual implementation.

In order to get the characteristic matrix of the thin-film the matrix multiplication is repeatedly applied using the characteristic matrix for each layer $c_j$:

<figure><img src="image_38.png" /><figcaption>Figure 33 Thin-film assembly</figcaption></figure>

$$ \bar{\mathcal{C}}^{s,p} = \bar{\mathcal{C}}_{cN}^{s,p} \cdots \bar{\mathcal{C}}_{c2}^{s,p} \cdots \bar{\mathcal{C}}_{c1}^{s,p} \qquad (50) $$

It is crucial, of course, to keep the correct order of the layers, as matrix multiplication is not commutative. The coatings are counted starting at the 2nd medium (see Figure 33). The characteristic matrix for every layer $j$ is given by:

$$ \bar{\mathcal{C}}_{j}^{s,p} = \begin{bmatrix} \cos \delta_{j} & \frac{i}{\tilde{\eta}_{j}^{s,p}} \sin \delta_{j} \\ i \tilde{\eta}_{j}^{s,p} \sin \delta_{j} & \cos \delta_{j} \end{bmatrix} \qquad (51) $$

The characteristic matrix transforms the complex amplitudes of the electric and magnetic field according to the layer's electrodynamic properties and thickness $d_j$. The phase factor $\tilde{\delta}_j$ results of the ratio of the effective optical thickness $\tilde{\Delta}_j$ and the wavelength $\lambda$:

$$ \tilde{\delta}_j = 2\pi \frac{\tilde{\Delta}_j}{\lambda} \qquad \tilde{\Delta}_j = d_j \tilde{n}_j \cos \tilde{\theta}_j \qquad (52) $$

The specific propagation angles $\tilde{\theta}_j$ and the corresponding effective refractive indices $\tilde{\eta}_j^{s,p}$ for each layer $j$ can again be found by applying Snell's law, eq. (46).

Finally, the reflectance coefficient $R$ and transmittance coefficient $T$ to calculate the power of the reflected and the transmitted beam for each polarisation state can again be found by:

$$ R = |\tilde{r}|^2 \qquad T = Re \left[ \frac{\tilde{\eta}_2}{\tilde{\eta}_1} \right] |\tilde{t}|^2 \qquad (53) $$

The additional factor required to calculate the transmittance coefficient considers that the impedance, as well as the propagation direction, is generally different in the two media. In the further text, these coefficients will be referred to as generalized Fresnel coefficients. In the following application, it is essential to additionally consider the phase change $\epsilon$ of the reflected and transmitted beam. They are simply represented by the arguments of the complex coefficients $\tilde{r}$ and $\tilde{t}$:

$$ \epsilon_r = \arg(\tilde{r}) \qquad \epsilon_t = \arg(\tilde{t}) \qquad (54) $$

The corresponding coefficient for absorptance $A$ can be found by applying the law of energy conservation:

$$ A = 1 - R - T \qquad (55) $$

Note that, again, all relations in (53), (54) and (55) are valid for each polarisation state separately; for the purpose of clarity, the indices s and p were omitted here.
4.7.2. Implementation

The implementation of above formulae is performed in the module TFresnel of the library rc_fresnel :

type T22cmat = record m11,m12,m21,m22 : complex;End;
type TeightLayerThicknessList = array [1..8] of Double;
type TeightLayerRIList = array [1..8] of Complex;
type TTRpolRecord = record T,R : TSPDouble;tamp,ramp : TSPComplex;rcross,tcross : Complex;End;
type
TFresnel = class(TObject) { TFresnel }
private
lambda : Double; // wavelength
nm,ns : complex; // ref.index medium (med1) and substrate (med2)
nLay : TeightLayerRIList; // thin-film refract. (1..substrate side, N..medium side)
layer-order
dLay : TeightLayerThicknessList; // thin-flim thickness (1..substrate side, N..medium side)
layCnt : Integer; // number of thin films considered
function eta(const sPol : Boolean; const n : complex; const theta : complex) : complex;
function deltaLayer(const n,d,theta : complex) : complex;
function layerMat(const sPol : Boolean; const n,d,theta : complex) : T22cmat;
function matMult(const a,b : T22cmat) : T22cmat;
function unitMat : T22cmat;
function snellThetaT(n1,n2 : complex; thetaI : complex) : complex;
public
constructor Create;
procedure setup(const layCnt_ : Integer); // initialize material interface
procedure updateLambdaNmNs(const lambda_,nmm,nmk,nsn,nsk : Double); // set ref.indicies & WL
procedure setLayThickness(const no : Integer; const d_ : Double); // thin-film thickness
procedure updateLayNo(const no : Integer; const n_,k_ : Double); // thin-film ref.ind.
procedure getRandTMuMat(const TRpol: TTRpolRecord; var RMM, TMM: TMuMat);
procedure getRT(const thetaI : Double;var TRpol : TTRpolRecord);
procedure getRTandThetaS(const thetaI: Double; var TRpol: TTRpolRecord;
var thetaS: Complex);
function reflectionMuMat(const TRpol : TTRpolRecord) : TMuMat;
function transmissionMuMat(const TRpol : TTRpolRecord) : TMuMat;
depositor Destroy; override;
end;

The module contains, on the one hand, private functions which perform the actual calculation in line
with the relations given above. The operator overloading features of the Pascal language allows a
convenient and transparent implementation of the relevant equations, as the complex-number alge-
bra can be converted into programming code straightforwardly. On the other hand, the public func-
tions of the module allow the initialisation of the materials and the direct calculation of the Müller
matrices for reflection and transmission (see section 4.9). Alternatively, only the parameters required
to derive the Müller matrices can be generated. These essential parameters are contained in the type
definition TTRpolRecord and reflect the results of equations (48) and (53). In the actual implementa-
tion of the raytracing (see chapter 7), the function getRTandThetaS is called to perform the Fresnel
calculations, while the function getRandTMuMat is called to establish the Müller matrices of the calcu-
lated parameters.
4.7.3. Testing and validation

The full validation is implicitly included in the Müller matrix section (4.9.3), as the Müller matrix
elements rely on the Fresnel calculations. During the implementation, the Fresnel module was addi-
tionally tested based on commonly known material behaviours.

Exemplary, Figure 34 illustrates the angular reflectivity for three different materials, resolved for three different wavelengths and the main polarisation states. As can be seen, the dielectric (i.e. non-conducting) glass exhibits limited wavelength dependence and a pronounced Brewster angle at 57° where only s-polarised light is reflected off the surface. Gold shows high reflectivity in the visual and infrared spectral range; however, less than half of the incident radiation is reflected in the UV range. Aluminium generally exhibits a high reflectivity for all wavelengths considered.

<figure><img src="image_39.png" /><figcaption>*Figure 34 Angular reflection of material/air interfaces for three wavelengths and two polarisation states.*</figcaption></figure>
# 4.8. Stokes vectors

## 4.8.1. Background and theory

Basically, two different, complementary representations are used to describe the state of polarisation of light, i.e. electromagnetic waves. One is called the *Jones vector representation*, while the other is referred to as the *Stokes vector representation*. Since the Jones representation is based on the electric field vector, it can consider the coherent superposition of waves but is limited to representing fully polarised light. In contrast, the Stokes notation cannot be used to model interference effects but is inherently able to cover unpolarised or partially polarised radiation. As mentioned in the introduction, the *RadiCal* approach considers effects relying on coherence only for interactions at material boundaries (including thin-film effects) but not for propagation within a medium. Still, the state and angle of polarisation of the traced light beams are crucial information for the light-surface interactions and consequently have to be tracked during the raytracing process. The Stokes representation is the ideal model for this purpose and has, thus, been chosen for the implementation of the raytracer. The Stokes parameters are the four components $S_i$ of the Stokes vector $\vec{S}$. Alternatively, they are often denominated as $I, Q, U, V$. They can be defined as (Kliger et al., 1990; Bass et al., 1995; Hecht, 2015):

$$ \vec{S} = \begin{pmatrix} S_0 \\ S_1 \\ S_2 \\ S_3 \end{pmatrix} = \begin{pmatrix} \langle E_x^2 \rangle + \langle E_y^2 \rangle \\ \langle E_x^2 \rangle - \langle E_y^2 \rangle \\ \langle 2E_x E_y \cos\delta \rangle \\ \langle 2E_x E_y \sin\delta \rangle \end{pmatrix} = \begin{pmatrix} I_{0\circ} + I_{90\circ} \\ I_{0\circ} - I_{90\circ} \\ I_{45\circ} - I_{135\circ} \\ I_{CR} - I_{CL} \end{pmatrix} \quad (56) $$

$I_{0\circ}, I_{45\circ}, I_{90\circ}, I_{135\circ}$ ... (measured) linearly polarized irradiance values [W/m²]

$I_{CR}, I_{CL}$ ... (measured) right, left hand circularly polarized irradiance values [W/m²]

$E_x, E_y$ ... orthogonal amplitude of electric field vector [V/m]

$(\cdot, \cdot)$ ... expectation value / time average

$\delta$ ... phase angle inbetween $E_x$ and $E_y$ [rad]

The definition based on the amplitudes of the electric field vector ($\langle E^2 \rangle$) links the Stokes parameters to the electrodynamic fundamentals. The second definition is based on irradiance values ($I$), measured behind ideally linearly and circularly polarising filters. This definition reflects the empirical approach.

While the definition of the Stokes vector might seem impractical at first glance, it indeed provides a very efficient method to describe the polarising or depolarising effects to which light beams can be subjected, as is shown below. Based on its definition, the primary polarisation states are represented by the normalised Stokes vectors provided in Table 6.

<figure><img src="image_40.png" /><figcaption>Table 6 Stokes vectors for basic polarisation states</figcaption></figure>
Additional (scalar) parameters used to describe the state of polarisation can easily be calculated from the Stokes parameters by (Bass et al., 1995; Garcia-Caurel et al., 2013):

$$ 
\begin{aligned}
DOP &= \frac{\sqrt{S_1^2 + S_2^2 + S_3^2}}{S_0} \\
\eta_{max} &= \frac{1}{2} \arctan \left( \frac{S_2}{S_1} \right)
\end{aligned} 
\qquad (57) $$

DOP ... degree of polarisation

DOLP ... degree of linear polarisation

$\eta_{max}$ ... angle of max.polarisation

Since the Stokes parameters are usually provided in vector form, it is convenient to perform transformations of them by applying matrix multiplication (see the following chapter).

Ellipsometry is a powerful method to investigate the electrodynamic properties, i.e. the complex refractive index, of material surfaces or thin-flims. For this field, a third definition of the Stokes vector is provided by introducing the ellipsometric angles $\Psi$ [rad] and $\Delta$ [rad] as:

$$ \Psi = \arctan\left(\frac{E_x}{E_y}\right) \qquad \Delta = \phi_x - \phi_y \qquad (58) $$

$E_x, E_y$ ... orthogonal amplitude of electric field vector [V/m]

$\phi_x, \phi_y$ ... initial phases of $\vec{E}_x, \vec{E}_y$ [rad]

$\delta$ ... phase angle inbetween $E_x$ and $E_y$ [rad]

$\Psi$ represents the amplitude ratio of the orthogonal components, whereas $\Delta$ is the phase difference between them. While the phase difference as already been introduced using the lowercase Greek letter, the capital letter will be used here to stay in line with ellipsometric conventions. Using the ellipsometric angles, equation (56) can be rewritten as:

$$ \vec{S} = \begin{pmatrix} S_0 \\ S_1 \\ S_2 \\ S_3 \end{pmatrix} = \begin{pmatrix} \langle E_x^2 \rangle + \langle E_y^2 \rangle \\ -S_0 \cos(2\Psi) \\ S_0 \sin(2\Psi) \cos(\Delta) \\ -S_0 \sin(2\Psi) \sin(\Delta) \end{pmatrix} \qquad (59) $$

The inversion of this definition allows to express the ellipsometric angles by the Stokes parameters. This reflects the empirical process, as the Stokes components are measured in ellipsometry:

$$ \Psi = \frac{1}{2} \arctan \left( \frac{\sqrt{S_2^2 + S_3^2}}{S_1} \right) \qquad \Delta = \arctan \left( \frac{-S_3}{S_2} \right) \qquad (60) $$

The normalised Stokes vector $\vec{S}'$, that describes the polarisation state independent of the total irradiance, is defined as:

$$ \vec{S}' = \frac{1}{S_0} \vec{S} = \begin{pmatrix} 1 \\ S_1/S_0 \\ S_2/S_0 \\ S_3/S_0 \end{pmatrix} \qquad (61) $$
### 4.8.2. Implementation

The implementation of the Stokes vector definition is trivial. A new type called `TStokesVect`, representing an array of 4 scalar values (of type `double`), is defined in the `rc_globals` module. Additionally, the definitions for the basic polarisations states are included in the same module:

```c
Type TStokesVect = array[0..3] of Double;

Const unpolStokes : TStokesVect = (1,0,0,0);
Const horizStokes : TStokesVect = (1,1,0,0);
Const vertStokes : TStokesVect = (1,-1,0,0);
Const p45Stokes : TStokesVect = (1,0,1,0);
Const m45Stokes : TStokesVect = (1,0,-1,0);
Const rightCircStokes : TStokesVect = (1,0,0,1);
Const leftCircStokes : TStokesVect = (1,0,0,-1);
```

For the elementary manipulation of Stokes vectors and the calculation of characteristic parameters according to equations (57), the following methods have been implemented in the globally applied module `rc vectors`.

```pascal
function SVect(const s0,s1,s2,s3:Double): TStokesVect;
function SVadd(const s1,s2:TStokesVect): TStokesVect;
function SVNormalize(const s:TStokesVect): TStokesVect;
function SVMult(const f:Double; const s:TStokesVect): TStokesVect;
function SVtotPolDegree(const s:TStokesVect): Double;
function SVtotPolIntensity(const s:TStokesVect): Double;
function SVlinPolDegree(const s:TStokesVect): Double;
function SVangleMaxPol(const s:TStokesVect): Double;
function SVrotPol(const s: TStokesVect; const negAngle: Boolean; const cosT: Double)
  : TStokesVect;
function SVhorizPolIntensity(const s:TStokesVect): Double;
function SVvertPolIntensity(const s:TStokesVect): Double;
function SVp45PolIntensity(const s:TStokesVect): Double;
function SVM45PolIntensity(const s:TStokesVect): Double;
function SVellipsometryDELTAofStokes(const s:TStokesVect): Double;
function SVellipsometryPSIofStokes(const s:TStokesVect): Double;
```

### 4.8.3. Testing and validation

Simple tests of the methods were performed during the implementation process. Extensive validation is included in the following section.
# 4.9. Müller matrix operations

## 4.9.1. Background and theory

Since the Stokes representation uses a 4D vector to describe the polarisation state, a 4x4 matrix can be used to model scattering events that alter the polarisation state of the beam. Generally, these matrices are referred to as *Müller matrices* (after Hans Müller, who introduced them in 1943). Müller matrices are primarily used for evaluating empirical data gathered with optical methods (e.g. by ellipsometry) or for designing optical systems when polarisation effects have to be considered. In such cases, the polarisation properties of the optical elements used (e.g. linear polariser) are characterised by their specific Müller matrices.

For the *RadiCal* raytracer, currently, only two Müller matrices are relevant. However, these are able to cover a wide range of polarisations effects as the implemented, generalized approach for computing the Fresnel coefficients (see section 4.7.1) also provides information for the phase change occurring on the scattering event ($\epsilon_p, \epsilon_s$). The two implemented Müller matrices are used to calculate the polarisation states, as well as the intrinsically contained powers of the (specularly) reflected and transmitted beam at interfaces of materials. Additionally, coatings (e.g. low-E, antireflective) or naturally occurring thin-films (e.g. oxide layers) can be considered, as they can be modelled with the implemented approach to determine the generalized Fresnel coefficients and resulting phase shifts. In general, the reflection of an unpolarised beam on a surface will lead to a transmitted and reflected beam. The polarisation state, as well as the intensities of the incident, the reflected and transmitted beam are described by their Stokes vec-

tors $\vec{S}_i, \vec{S}_r$ and $\vec{S}_t$.

The Müller matrices to calculate the induced components ($\vec{S}_r, \vec{S}_t$) of the incident beam ($\vec{S}_i$), can be derived by demanding continuity of the electric field, as well as linear relations among the incident, reflected and transmitted beam governed by the generalized Fresnel coefficients, see Eq. (49). The derivation can, for example, be found in a general approach formulated by Chipman (see Bass et al., 1995, chapter “linear diattenuator and retarder”) and Stenflo (1994a) or more specifically, and more recently, by Garcia (2011).

Using the *Müller matrix for reflection*

$\bar{M}_r$, the Stokes vector of the reflected beam $\vec{S}_r$, based on the generalized Fresnel coefficients (see eq. (53) and (54)) can be found by:

<figure><img src="image_41.png" /><figcaption>Figure 35 Reflection and transmission on material interface considering polarization</figcaption></figure>

$$ \vec{S}_r = \bar{M}_r \cdot \vec{S}_i = \begin{pmatrix} \frac{1}{2}(R_s + R_p) & \frac{1}{2}(R_s - R_p) & 0 & 0 \\ \frac{1}{2}(R_s - R_p) & \frac{1}{2}(R_s + R_p) & 0 & 0 \\ 0 & 0 & -\sqrt{R_p R_s} \cos(\epsilon_p - \epsilon_s) & -\sqrt{R_p R_s} \sin(\epsilon_p - \epsilon_s) \\ 0 & 0 & \sqrt{R_p R_s} \sin(\epsilon_p - \epsilon_s) & -\sqrt{R_p R_s} \cos(\epsilon_p - \epsilon_s) \end{pmatrix} \begin{pmatrix} S_0 \\ S_1 \\ S_2 \\ S_3 \end{pmatrix} \quad (62) $$
Likewise, by using the *Müller matrix for transmission* $\bar{M}_t$, the Stokes vector of the transmitted beam $\vec{S}_t$ is determined by:

$$ \vec{S}_t = \bar{M}_t \cdot \vec{S}_i = \begin{pmatrix} \frac{1}{2}(T_s + T_p) & \frac{1}{2}(T_s - T_p) & 0 & 0 \\ \frac{1}{2}(T_s - T_p) & \frac{1}{2}(T_s + T_p) & 0 & 0 \\ 0 & 0 & \sqrt{T_p T_s} \cos(\epsilon_p - \epsilon_s) & \sqrt{T_p T_s} \sin(\epsilon_p - \epsilon_s) \\ 0 & 0 & -\sqrt{T_p T_s} \sin(\epsilon_p - \epsilon_s) & \sqrt{T_p T_s} \cos(\epsilon_p - \epsilon_s) \end{pmatrix} \begin{pmatrix} S_0 \\ S_1 \\ S_2 \\ S_3 \end{pmatrix} \quad (63) $$

The Müller matrices $\bar{M}_r$ and $\bar{M}_t$ are a function of wavelength and incidence angle, as the generalized Fresnel coefficients used to calculate them are so. Consequently, they have to be computed individually for each scattering event.

While in the above formulation of the Müller matrix, the explicit dependence of the phase angles and amplitudes is easily accessible, an alternative formulation of the Matrix based on the complex amplitudes $\tilde{r}$ and $\tilde{t}$ allows a computationally more efficient evaluation.

$$ \bar{M}_r = \begin{pmatrix} \frac{1}{2}(|\tilde{r}_s|^2 + |\tilde{r}_p|^2) & \frac{1}{2}(|\tilde{r}_s|^2 - |\tilde{r}_p|^2) & 0 & 0 \\ \frac{1}{2}(|\tilde{r}_s|^2 - |\tilde{r}_p|^2) & \frac{1}{2}(|\tilde{r}_s|^2 + |\tilde{r}_p|^2) & 0 & 0 \\ 0 & 0 & \mathrm{Re}[\tilde{r}_s \tilde{r}_p^*] & \mathrm{Im}[\tilde{r}_s \tilde{r}_p^*] \\ 0 & 0 & -\mathrm{Im}[\tilde{r}_s \tilde{r}_p^*] & \mathrm{Re}[\tilde{r}_s \tilde{r}_p^*] \end{pmatrix} \quad (64) $$

$$ \bar{M}_t = \begin{pmatrix} \frac{1}{2}(|\tilde{t}'_s|^2 + |\tilde{t}'_p|^2) & \frac{1}{2}(|\tilde{t}'_s|^2 - |\tilde{t}'_p|^2) & 0 & 0 \\ \frac{1}{2}(|\tilde{t}'_s|^2 - |\tilde{t}'_p|^2) & \frac{1}{2}(|\tilde{t}'_s|^2 + |\tilde{t}'_p|^2) & 0 & 0 \\ 0 & 0 & \mathrm{Re}[\tilde{t}'_s \tilde{t}'_p^*] & \mathrm{Im}[\tilde{t}'_s \tilde{t}'_p^*] \\ 0 & 0 & -\mathrm{Im}[\tilde{t}'_s \tilde{t}'_p^*] & \mathrm{Re}[\tilde{t}'_s \tilde{t}'_p^*] \end{pmatrix} \quad (65) $$

In order to consider the different angles of propagation and different impedance in the two media, a modified amplitude coefficient for transmission $\tilde{t}'$ is used here. It is defined as

$$ \tilde{t}'_{s,p} = \mathrm{Re} \left[ \frac{\tilde{\eta}_{s,p}^2}{\tilde{\eta}_{s,p}^2} \right] \tilde{t}_{s,p} \quad (66) $$

The phase information is now contained in the complex-valued amplitudes. Hence, these alternative formulations of the Müller matrices do not rely on computationally expensive trigonometric functions (sine and cosine). Further, the explicit determination of the phase-angles (see (54)) can also be omitted now. This additionally increases the algorithm's performance, as the explicit evaluation of the two phase-angles requires the computationally expensive use of the arctangent function.

By definition, the generalized Fresnel coefficients are only valid for radiation polarised strictly parallel ($R_p, T_p$) and perpendicular ($R_s, T_s$) to the plane defined by the surface normal and the propagation direction. Hence, the Stokes vector for the incident radiation, defined for an arbitrary orthogonal coordinate system $\vec{p}', \vec{s}',$ has to be provided in the aligned coordinate system $\vec{p}, \vec{s}$ determined by the surface's normal vector $\vec{n}$ (see Figure 35). To perform this, the Stokes vector has to be transformed accordingly. After the corresponding (clockwise) rotation angle $\varphi$ is determined, the aligned Stokes vector $\vec{S}$ can simply be calculated by:

$$ \vec{S} = R(\varphi) \cdot \vec{S}' = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & \cos(2\varphi) & \sin(2\varphi) & 0 \\ 0 & -\sin(2\varphi) & \cos(2\varphi) & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} S_0' \\ S_1' \\ S_2' \\ S_3' \end{pmatrix} \quad (67) $$
The rotation matrix $R(\varphi)$ is not referred to as Müller matrix here, as the actual polarisation state and the total intensity of the Stokes vector are invariant to this transformation. Solely the reference coordinate system of the Stokes vector is rotated.

### 4.9.2. Implementation

In the implementation, the Müller matrix class `TMuMatModule` in the library `rc_mumat` has been chosen to encapsulate all data and classes needed to perform the Fresnel and Stokes operations relevant for a specific material interface.

```typescript
Type TMuMat = array[1..4,1..4] of Double; // 4x4 Müller matrix
Type TTRpolRecord = record T,R : TSPDouble;tamp,ramp : TSPComplex;rcross,tcross : Complex;End;
Type
  TMuMatModule = class(SharpObject) { TMuMatModule }
  private
    Fr : TFresnel; // instance of TFresnel class for calc.
    matMedium : TRIbaseclass; // material data for mat1 (incidence)
    matSubstrate : TRIbaseclass; // material data for mat2
    mathinFilmCnt : Integer; // number of thin-films (max 8) (0=none)
    mathinFilm : TEightLayerRIBaseClassList; // material data for thinfilms
    mathinFilmD : TEightLayerThicknessList; // thicknesses of thinfilms
  public
    constructor Create;
    procedure setup(const calcBiCu: Boolean; const MapSettings: TMuMatMapSettings;
                    const withTrans: boolean; var matSubstrate, matMedium: TRIbaseclass;
                    const mathinFilmCnt: Integer; var mathinFilm: TEightLayerRIBaseClassList;
                    var mathinFilmD: TEightLayerThicknessList);
    function getExplicitMuMatR(const angle, w1: double): TMuMat;
    function getExplicitMuMatT(const angle, w1: double): TMuMat;
    procedure stokesCalcRT(const angle, w1: double; const s: TS Stokes Vect;
                          var sRefl, sTrans: TS Stokes Vect; var n1Divn2: double; var alphaAbs1,
                          alphaAbs2: Double); // returns the stokes vector for a reflected and transmitted
  beams
    procedure getMuMatRTangle(const angle, w1: double; var RMM, TMM: TMuMat;
                                 var thetaS: Complex; var n1Divn2: Double; var alphaAbs1, alphaAbs2: Double );
    destructor Destroy; override;
end;
```

As can be seen in the header, one instance of the `TFresnel` class (containing the methods to perform the Fresnel calculations), as well as several instances of `TRIbaseclass` (containing the material data), are included and used within the `TMuMatModule` class.

While the methods `getExplicitMuMatR` and `getExplicitMuMatT` are used for validation purposes, the raytracer solely calls the key method `getMuMatRTangle`, which performs all necessary operations to finally provide the Stokes vectors for the reflected and transmitted beam. The method only performs the matrix multiplications (equations (62) and (63)) and, in turn, calls the method `stokesCalcRT` to get the required values for the Müller matrices in the form of a `TTRpolRecord`. As mentioned above, the necessary parameters for the evaluation of the Müller matrices are provided in complex-valued form. Further, to avoid redundant calculation and minimize the calculation time, the functions `get-MuMatRTangle` and `stokesCalcRT` pass additional parameters of the Fresnel calculations relevant for the raytracing process (see inline comments in the source):

```pascal
procedure TMuMatMap.getMuMatRTangle(
    const angle,w1: double; // input: incident angle, wavelength
    var RMM,TMM: TMuMat; // output: Müller matrices for refl. and trans. beam
    var thetaS: Complex; // output: outgoing/refracted angle
    var n1Divn2: double; // output: fraction of n1/n2
    var alphaAbs1, alphaAbs2: Double); // output: (beer) attenuation coeff., both media
var
    c: Integer;
n,k: double;
nm,km: double;
TRpol: TTRpolRecord;
ang: double;

Begin
    // get and n,k for material 2
```
```csharp
if matSubstrate = nil then Begin n := NairConst; k := 0; End // use default air const
else matSubstrate.getnk(wl*1E6, n, k);
// get and n, k for material 1 (incidence)
if matMedium = nil then Begin nm := NairConst; km := 0; End // use default air const
else matMedium.getnk(wl*1E6, nm, km);
// return additional parameters for further raytracing calculation
nldivn2 := nm/n; // medium to substrate for VRefract function
alphaAbs1 := (4*pi*km)/wl; // get attenuation coeff. acc. to beers law for mat1
alphaAbs2 := (4*pi*k)/wl; // get attenuation coeff. acc. to beers law for mat2

Fr.updateLambdaNmNs(wl, nm, km, n, k);

For c:=1 to mathinFilmCnt do
    Begin
        MathinFilm[c].getnk(wl*1E6, n, k);
        Fr.updateLayNo(c, n, k);
    End;

Fr.getRTandThetaS(angle, TRpol, ThetaS);
Fr.getRandTMuMat(TRpol, RMM, TMM);
End;
```

### 4.9.3. Testing and validation

While the standard tests during implementation have been performed with due care, extended validation tests were performed to ensure the validity and functionality of the implemented code. The calculation of the outgoing Stokes vectors by multiplication with the corresponding Müller matrices is trivial. However, the determination of the Müller matrices depends fundamentally on the material module `TRIbaseclass` and particularly on the `TFresnel` module, which provides the necessary coefficients. Hence, the following examples prove the validity of all three modules.

As mentioned above, the detailed modelling of optics considering polarisation is currently mainly applied for advanced instrument design, astronomy and optical measurement systems (e.g. ellipsometry). Therefore the following examples are mainly of these fields. The complexity of the following validation cases gradually increases.

#### 4.9.3.1. Validation Case 1 - Simple Silver Mirror
[constant wavelength, constant angle, no thin-films)]

<figure><img src="image_42.png" /><figcaption>Figure 36 "Example 1" of Tetsu Anan working paper (Anan, 2018)
with reference to Stenflo (Stenflo, 1994b)</figcaption></figure>

This example has been found in a working paper by Tetsu Anan (2018) in which he evaluates the validity of calculations performed with the Zemax OpticsStudio software, which is widely used for advanced optical instrumentation design and also offers polarimetric raytracing features. In the paper, the Müller matrix of a reflection on a 45° tilted mirror is calculated. The precise wavelength and refractive index are provided in the example (see Figure 36).

In order to demonstrate how compact the validation can be done based on the implemented libraries, the entire source code is provided here:

```csharp
mmat := TRCMaterial.Create('Silver mirror'); // create a new material interface
```
```
with mmat do
Begin
  Medium:=nil; // medium1: air interface
  // medium2: Silver, use constant model to provide apply the specified values directly
  Substrate:=TRIconstNK.Create('Silver const.', '', '',
                             'n-const: 0.150160'+#13#10+'k-const: 3.4727');
  ThinFilmCnt:=0; // no coating considered
  SetupMMM(false); // initialize Mueller matrix module
  RMuellerMat:=mmat.MMM.getMuMatRprecise(45*deg2rad,567E-9); // calculate MM for 45°/567nm
  Result:=dumpMuMat(normalizeMM(RMuellerMat)); // dump the normalized matrix
end;
```

The resulting Müller matrix for reflection is:

<table><tr><td>1,000000</td><td>0,015913</td><td>0</td><td>0</td></tr><tr><td>0,015913</td><td>1,000000</td><td>0</td><td>0</td></tr><tr><td>0</td><td>0</td><td>-0,923511</td><td>0,383242</td></tr><tr><td>0</td><td>0</td><td>-0,383242</td><td>-0,923511</td></tr></table>

Considering the number of digits provided, the result perfectly matches the Zemax calculation results and Anan's verification thereof (see Figure 36).

### 4.9.3.2. Validation Case 2 – Aluminium mirror with an oxide layer [constant wavelength, variable angles, thin-film]

In a publication by van Harten *et al.* (2009) the polarisation properties of a “real mirror” are modelled using the Müller matrix and thin-film approach as applied in *RadiCal*. In the paper, the effect of the (usually present) amorphous aluminium oxide layer and its thickness is investigated.

Again, the complex refractive index values for the bulk material (Aluminium: $\tilde{n} = 1.2 - 6.85i$) and the thin-film (Aluminium oxide: $\tilde{n} = 1.6 - 0i$) is provided as a constant value for the wavelength used (600 nm).

The thickness of the oxide layer of the final mirror was determined as $d = 4.12nm$.

These parameters were used for the implementation of the model with *RadiCal*. The entire source code used for the validation is provided for this example as well:

```
mmat:=TRCMaterial.Create('Al/AlO "real" mirror'); // create a new material interface
with mmat do
Begin
  // medium1: air/vacuum interface
  Medium:=nil;
  // medium2: Aluminium, use constant model to provide apply the specified values directly
  Substrate:=TRIconstNK.Create('Al const.', '', '', 'n-const: 1.20'+#13#10+'k-const: 6.85');
  ThinFilmCnt:=1; // 1 coating considered
  // thin-flim: AlxOx, use constant model to provide apply the specified values directly
  ThinFilm [1]:=TRIconstNK.Create('AlO porous', '', '', 'n-const: 1.6'+#13#10+'k-const: 0.0');
  ThinFilmD[1]:=4.12*1E-9;

  SetupMMM(false); // initialize Mueller matrix module
  for c:=0 to 89 do // iterate incident angles from 0° to 89°
    Begin
      // calculate MM for c° and 600nm
      RMuellerMat:=mmat.MMM.getMuMatRprecise(c*deg2rad,600E-9);
      //dump the normalized matrix inline
      Result:=Result+dumpMuMatInline(normalizeMM(RMuellerMat))+#13#10;
    end;
end;
```

Since van Harten *et al.* have provided the Müller matrix elements in a commonly used table of charts representation, it is possible compare the results for the six non-zero matrix elements. For this purpose, the *RadiCal* calculation results are superimposed on the original charts. The red curves and grid lines represent the results of the *RadiCal* calculation; the black and grey elements are the original data of the publication.
<figure><img src="image_43.png" /><figcaption>Figure 37 Angular dependence of Müller matrix elements computed with RadiCal (red line), measured data (dots), van Harten et al. (van Harten et al., 2009): calculation including oxide-layer (black line) and excluding oxide-layer (grey line)</figcaption></figure>

As can be seen, the *RadiCal* calculations show a perfect match with the calculations of van Harten *et al.* Further, it is evident (see elements *m12* and *m21* in Figure 37) that the ability to include thin-films increases the potential and accuracy of the model considerably (compare grey curve vs. red/black curves and dots in Figure 37).

#### 4.9.3.3. Validation Case 3 – A space mirror (*Sciamatchy* on Envisat)
[wavelength range, two angles, two layer thin-film]

In their publication ‘Mirror contamination in space I: mirror modelling’ J.Krijger *et al.* (2014) have extended the approach of van Harten *et al.* (2009) by assuming an additional contamination layer in the form of an oil thin-film. The new model should be able to better predict the mirror’s reflectance in the UV-range, where the simple model by van Harten showed the largest deviations. Based on their two-film model (oil and oxide), they provide reflectance data for the mirror for two angles of incidence and in a wavelength range from 300 nm to 1500 nm. The modelled data is compared against measurements performed by ESA for the calibration of the aluminium scan mirror of the *Sciamatchy* spectrometer of the *Envisat* satellite.

Since they do not provide the entire information regarding the refractive index functions for oil, the amorphous aluminium oxide and the bulk aluminium in the relevant wavelength range, the *RadiCal* implementation is used to redo the work performed in the paper, i.e. based on the given refractive index functions, the optimal thin-film thicknesses for the two-layer model is determined. This approach also allows testing the TRIClasses used to import refractive index information from external sources. The datasets of Rakic (1995) for Aluminum, Malitson *et al.* (1958) for Al₂O₃ and Moutzouris *et al.* (2014) for benzene (based on corresponding RefractiveIndex.INFO (Polyanskiy, 2022) database files) were applied.

Only the implementation of the material definition, based on the imported datasets, is shown below:

```cpp
// medium1: air/vacuum interface
Medium:=nil;

// medium2: Aluminium, Rakic dataset
Substrate:= riLib.ReadRI('Al','\main\Al\Rakic.yml',err);
ThinFilmCnt:=2; // 2 coating considered

// thin-flim 1: Al2O3, Malitson et al.
ThinFilm [1]:=riLib.ReadRI('Al2O3','\main\Al2O3\Malitson-o.yml',err);

// thin-flim 2: benzene, Moutzouris et al.
ThinFilm [2]:= riLib.ReadRI('Benzene','organic\C6H6 - benzene\Moutzouris.yml',err);
```

Of course, whether the chosen refractive index functions match the real thin-film properties is debatable, as thin-films of only a few nanometers usually show deviations from the behaviour of the
bulk material. However, this argument is also relevant to the assumptions regarding refractive index and layering in the publication by Krijger et al..

After optimisation, the surprisingly best fitting result, measured by the average RMS of all data points, is a thickness of *d* = 5.3 nm for the oil film and the elimination (*d* = 0 nm) of the oxide layer. The dispersion curves of both thin-films are of similar shape; however, benzene shows a steeper rise of the refractive index towards the UV wavelengths, which seems to lead to better results in this region. Physically, the 5.3 nm thin-film of benzene represents the single layer of contaminated aluminium oxide but is not an actual homogeneous oil film.

<figure><img src="image_44.png" /><figcaption>Figure 38 Sciamachy scan mirror reflectance as a function of wavelength for two incidence angles (29° / 61°), for s- and p-polarisation. red dotted curves: RadiCal results, black dashed/solid curves: 2-film model (Krijger et al.), black dotted/dash-dotted curves: bulk aluminium model (Krijger et al.), markers: measured data (ESA)</figcaption></figure>

Again, the original chart is used to superimpose the data computed with *RadiCal*. As is shown in Figure 38, the results of the presented model match the results of Krijger *et al.* for wavelengths greater than 500 nm. Due to the different choices of film thickness and the different refractive index datasets, deviations can be seen in the UV range (300 nm measurement). While the S-61° measurement shows a good match, the other model results at 300 nm show a slight deviation to the upside, while the Krijger model generally shows deviations of the same scope to the downside.

The validation case proves that the *RadiCal* method is implemented correctly and able to perform advanced optical calculations, as applied in space technology.
4.9.3.4. Validation Case 4 – Determination of thin-film layer thickness
[thin-film with variable thickness, ellipsometry]

In the ‘Ellipsometry’ chapter of the
‘Encyclopedia of Spectroscopy and
Spectrometry’ Jellison (1999) shows
how ellipsometric measurements can
conveniently be used to determine the
layers of thin films. He provides a com-
monly applied example for the meas-
urement of silicon dioxide (SiO₂) layer
thicknesses on silicon wafers using a
Helium-Neon laser (λ = 633 nm). Since
Jellison provides the complex refractive
indices for SiO₂ (ñ = 1.46 – 0i) and
silicon (ñ = 3.86 – 0.018i) and the in-
cidence angle (θi = 60°), the imple-
mented libraries can be used to deter-
mine the ellipsometric angles and com-
pare the results with the original re-
sults. Equations (60) are used to com-
pute the ellipsometric angles Ψ and Δ.

Again, the original chart is presented, and the *RadiCal* calculation results are superimposed on it (see Figure 39). The results perfectly match the original calculations of Jellison. As can be seen, the characteristic variation of the ellipsometric angles Ψ and Δ allows a very precise determination of the thin-flim layer thickness from 0 up to a thickness of 283 nm. For higher thicknesses, the results will be cyclically mapped on the same graph with a modulo of that value.

<figure><img src="image_45.png" /><figcaption>Figure 39 Ellipsometric angles (Ψ, Δ) as function of SiO₂ film thickness;
black curve: original paper (Jellison, 1999);
red curve: RadiCal calculations</figcaption></figure>
### 4.9.3.5. Validation Case 5 – Ellipsometric measurements with variable angles [thin-film with variable thickness and incidence angles, ellipsometry]

In their publication 'The influences of roughness on film thickness measurements by Mueller matrix ellipsometry' Ramsey and Ludema (1994) have conducted research on the effect or layer thickness and roughness of various materials using a dual rotating-compensator Mueller matrix ellipsometer.

<figure><img src="image_46.png" /><figcaption>Figure 40 Müller matrix ellipsometer
original Figure of Ramsey and Ludema (1994)</figcaption></figure>

The measured data for the Müller matrix elements of three different angles of incidence and the angular variation of the elliptic angles for steel samples with five different coatings is replicated by performing “virtual ellipsometric measurements” using the *RadiCal* libraries. The computed data is then compared against the measurement data. The data of the five polished steel specimens 1A-1E coated with variable thickness of magnesium fluoride ($MgF_2$) was used for this comparison. The wavelength of the laser (632.8 nm), the incidence angles as well as the complex index of refraction values of the materials (steel: $\tilde{n} = 2.501 - 3.408i$, $MgF_2$: $\tilde{n} = 1.38 - 0i$) were defined in line with the parameters stated in the publication.

First, the full Müller matrix is calculated and compared for the steel sample coated with 315 nm $MgF_2$. The publication provides measured matrix elements for three different incidence angles. The original measured data, computed data and absolute deviation are presented in Table 7.

<table><caption>Table 7 Müller matrices of specimen 1E – measured data vs. RadiCal calculation results</caption><thead><tr><th colspan="5">measurement</th><th colspan="4">calculation</th><th colspan="2">abs. deviation</th></tr></thead><tbody><tr><td colspan="11">TABLE V. Typical Mueller matrices—specimen 1E [roughness=polished (0.01 µm R<sub>a</sub>) with 315 nm film (λ=632.8 nm)].</td></tr><tr><td colspan="2">Angle of incidence=45°</td><td></td><td></td><td></td><td colspan="2">1,0000</td><td>0,1568</td><td>0,0</td><td>0,0</td><td>0,0000</td><td>0,0033</td></tr><tr><td>1,0000</td><td>-0.1535</td><td>-0.0069</td><td>-0.0009</td><td></td><td>0,1568</td><td>1,0000</td><td>0,0</td><td>0,0</td><td>0,0</td><td>0,0031</td><td>0,0012</td></tr><tr><td>-0.1537</td><td>0.9988</td><td>0.0714</td><td>-0.0157</td><td></td><td>0,0</td><td>0,0</td><td>-0.8470</td><td>0,5080</td><td></td><td></td><td>0,0028</td><td>0,0138</td></tr><tr><td>-0.0063</td><td>0.0692</td><td>-0.8498</td><td>0.4942</td><td>×0.3523</td><td>0,0</td><td>0,0</td><td>-0.5080</td><td>-0.8470</td><td></td><td></td><td>0,0116</td><td>-0.0095</td></tr><tr><td>0,0014</td><td>0.0187</td><td>-0.4964</td><td>-0.8565</td><td></td><td>0,0</td><td>0,0</td><td>-0.5080</td><td>-0.8470</td><td></td><td></td><td>0,0116</td><td>-0.0095</td></tr><tr><td colspan="2">Angle of incidence=60°</td><td></td><td></td><td></td><td colspan="2">1,0000</td><td>0,3270</td><td>0,0</td><td>0,0</td><td>0,0000</td><td>0,0006</td></tr><tr><td>1,0000</td><td>-0.3264</td><td>-0.0167</td><td>0.0003</td><td></td><td>0,3270</td><td>1,0000</td><td>0,0</td><td>0,0</td><td>0,0002</td><td>0,0020</td></tr><tr><td>-0.3268</td><td>0.9980</td><td>0.0620</td><td>-0.0228</td><td>× 0.4072</td><td>0,0</td><td>0,0</td><td>-0.6331</td><td>0,7016</td><td></td><td></td><td>0,0012</td><td>0,0080</td></tr><tr><td>0,0117</td><td>0.0557</td><td>-0.6319</td><td>0.6936</td><td></td><td>0,0</td><td>0,0</td><td>-0.7016</td><td>-0.6331</td><td></td><td></td><td>0,0096</td><td>0,0120</td></tr><tr><td>0,0020</td><td>0.0263</td><td>-0.6920</td><td>-0.6451</td><td></td><td>0,0</td><td>0,0</td><td>-0.7016</td><td>-0.6331</td><td></td><td></td><td>0,0096</td><td>0,0120</td></tr><tr><td colspan="2">Angle of incidence=75°</td><td></td><td></td><td></td><td colspan="2">1,0000</td><td>0,6107</td><td>0,0</td><td>0,0</td><td>0,0000</td><td>0,0001</td></tr><tr><td>1,0000</td><td>-0.6106</td><td>-0.0285</td><td>0.0017</td><td></td><td>0,6107</td><td>1,0000</td><td>0,0</td><td>0,0</td><td>0,0003</td><td>0,0017</td></tr><tr><td>-0.6110</td><td>0.9983</td><td>0.0537</td><td>-0.0238</td><td>× 0.4291</td><td>0,0</td><td>0,0</td><td>-0.2638</td><td>0,7466</td><td></td><td></td><td>0,0038</td><td>0,0004</td></tr><tr><td>-0.0182</td><td>0.0401</td><td>-0.2600</td><td>0.7462</td><td></td><td>0,0</td><td>0,0</td><td>-0.7466</td><td>-0.2638</td><td></td><td></td><td>0,0007</td><td>0,0014</td></tr><tr><td>0,0018</td><td>0.0294</td><td>-0.7459</td><td>-0.2624</td><td></td><td>0,0</td><td>0,0</td><td>-0.7466</td><td>-0.2638</td><td></td><td></td><td>0,0007</td><td>0,0014</td></tr></tbody></table>
It can be seen that the calculated data shows
a good agreement with the empirical data.
Since the polarisation axes are defined differ-
ently, the m12/m21 elements show the op-
posite sign. Further, the empirical data also
shows non-zero elements in the upper-right
and lower-left quadrants, since the measure-
ment as well as the material exhibit non-
ideal properties. As stated in the original
publication, the measurement results ap-
proach the ideal values for higher incidence
angles. This is reflected in the mean absolute
deviation of the matrix elements for each of
the three incidence angles, which are 0.0057
for 45°, 0.0042 for 60° and 0.0011 for 75°.

Further, the ellipsometric angles Ψ (ampli-
tude ratio) and Δ (phase difference) were de-
termined from the Müller matrix for five dif-
ferent coatings in an angularly resolved way.
Again, the results of the virtual ellipsometric
measurements are compared against the
measured data, and again, the calculated re-

sults are superimposed on the original chart by Ramsey and Ludema (1994).

The comparison of the ellipsometric angle functions reveals a remarkably high level of agreement with the measured data (see Figure 41) as the ellipsometric angles are usually very sensitive to small changes of the influencing parameters (i.e. layer thickness, refractive indices, incidence angle).

<figure><img src="image_47.png" /><figcaption>Figure 41 Ellipsometric angles as function of angle
MgF₂ on steel samples with five coating thicknesses
black curves: original paper (Ramsey and Ludema, 1994);
red curves: RadiCal calculations</figcaption></figure>
# 4.10. Roughness model and microfacet theory

## 4.10.1. Background and theory

Based on Fresnel's equations and Snell's law, the theory presented so far allows accurate modelling of reflection and refraction (transmission) for many opaque and transparent surfaces, such as glass, metal, marble and other polished surfaces. The models will provide an accurate description of the scattering of light at the surface of these materials and, therefore, will also allow a photorealistic appearance of the surface if used for raytracing. However, the majority of surfaces show significant diffuse reflection. This effect cannot be described (or modelled in raytracing) by only applying Fresnel's approach, as this would lead to a purely specular reflection. Mainly two processes are responsible for diffuse reflection: surface roughness and subsurface scattering. The latter process will be covered in the next section, while this section covers surface roughness.

Surfaces can be modelled directly with the Fresnel approach as long as their boundary is ideally flat or ideally smooth, i.e. locally flat. While the production or finishing processes for glass and metals result in a smooth surface, most naturally occurring surfaces will not fulfil this condition. An important - and presumably the only - exception to this is the still surface of water, ice, or other liquids. While Fresnel's equation and Snell's law are still valid on a microscopic scale for most naturally rough surfaces, they cannot directly be applied by using the macroscopic orientation of the surface only. A statistical approach is required to model the uneven surface. If this distribution of the varying surface orientations is considered appropriately, most of the rough surfaces can still be modelled using the Fresnel approach. The implementation of diffuse reflection via microfacet modelling further allows considering the polarisation state of the diffusely reflected or refracted radiation, as the complete Fresnel calculation is applied for every micro-surface.

However, it should be noted here that the Fresnel model can also break down for some regular or nanostructured surfaces with structural dimensions in the range of the wavelengths, as in the derivation of Fresnel's equations, homogeneous electromagnetic properties are assumed for the two media that form the interface.

### 4.10.1.1. Microfacet models

The adaptation required for rough surfaces is performed by replacing the constant normal vector of the macrosurface $\vec{n}$ with statistically distributed normal vectors $\vec{m}$ describing the orientation of microsurfaces or micro-facets (Figure 42). This approach is called the *microfacet model*. The model can, of course, not be used to reconstruct the actual

<figure><img src="image_48.png" /><figcaption>Figure 42 Microfacet and macrosurface normal vectors</figcaption></figure>

topology of the surface accurately. It only provides distributions relevant for modelling the surface's directional reflection and transmission profiles, based on optical geometry assumptions. The model can therefore be regarded as randomly oriented, specularly reflecting facets (mirrors), whereby the orientations of the facets follow a distribution determined by a roughness parameter.

The microfacet model was introduced by Cook and Torrance (1982) for rendering in computer graphics. Based on this work, many authors have proposed variations or computationally cheaper implementations. Significant improvements were, e.g. provided by Ward (1992), Schlick (1994), Walter *et al.* (2007) as well as by Ashikhmin and Shirley (2011). Most of the early approaches, which
are still widely used, do not satisfy energy conservation since the main purpose of the models is to allow computer graphics rendering with high computational efficiency.

### 4.10.1.2. Shadowing and masking

While the different versions show some alterations, usually, the microfacet models consist of three components: a distribution of microfacets surface normals $D(\vec{m})$, a shadowing function and a masking function. The distribution of surface normals $D(\vec{m})$ depends on a roughness parameter and describes the orientations of the microfacets. In order to determine the visible facets $D(\vec{m})$ has to be masked by a masking function, which depends on the viewing angle in backward raytracing. In the forward raytracing approach applied in the *RadiCal* method, this corresponds to the angle of the incident light $\vec{v}_i$. Applying this masking function $G_1(\vec{v}_i, \vec{m})$ to the distribution of surface normals $D(\vec{m})$ the distribution of visible surface normal $D_v(\vec{m})$ is obtained. This distribution has to be sampled for any incoming ray to identify the correct orientation of the microfacet. The Fresnel coefficients are then determined for the reflection on the microfacet. If the ray is reflected, a second function, called the shadowing function, is applied to detect if another microfacet blocks the outgoing ray. Since in computer graphics, usually only one collision (or “bounce”) of the ray on a surface is considered, the two functions are often combined to form a shadowing-masking function $G_2(\vec{v}_i, \vec{v}_o, \vec{m})$. This function depends on the incoming and outgoing direction and determines if the considered ray is unobstructed and, thus, should be used to sample the surface's appearance.

In a simple approach, the masking function $G_1$ can be used to construct the shadowing-masking function

$$G_2(\vec{v}_i, \vec{v}_o, \vec{m}) = G_1(\vec{v}_i, \vec{m}) \cdot G_1(\vec{v}_o, \vec{m}) \qquad (68)$$

In this approach, the model is referred to as uncorrelated surface, as the orientations of neighbouring microfacets are assumed to show no correlation. On a more detailed look, it becomes evident that depending on the assumed topology, correlation relations must be considered. Therefore, the shadowing function for the outgoing ray has to depend on the orientation of the hit microfacet. For such correlated surfaces, equation (68) will not hold, and $G_2(\vec{v}_i, \vec{v}_o, \vec{m})$ has to be provided in closed form. In the implementation of the *RadiCal* raytracer, forward raytracing is applied. Contrary to the standard rendering approach using back-tracing, the masking function determines the illuminated microfacet, whereas the shadowing function determines if the reflection path is unobstructed. In the *RadiCal* approach, every ray is followed from “cradle to grave” to ensure energy conservation. Hence, unlike in the standard rendering methods, it is not possible to simply omit obstructed light paths as usually applied with shadowing-making functions. Since usually only the implicit shadowing-masking function and masking functions are provided, they are used to evaluate the shadowing function $G_s$:

$$G_s(\vec{v}_i, \vec{v}_o, \vec{m}) = \frac{G_2(\vec{v}_i, \vec{v}_o, \vec{m})}{G_1(\vec{v}_i, \vec{m})} \qquad (69)$$

Expressed in terms of probability, it can be stated that the probability of an obstruction in the outgoing direction ($G_s$) is obtained by dividing the probability of obstruction on either path ($G_2$) by the probability of the obstruction for the incoming path ($G_1$).
### 4.10.1.3. Implementation of the Heitz Smith GGX model

While many different BSDF models are available, most do not meet the requirements of the implemented Monte-Carlo method. For rendering purposes, usually, importance sampling is applied, i.e. many rays are computed and weighted, or rejected, according to the masking and shadowing functions. For the use in the *RadiCal* raytracer, however, it is essential to use a model which explicitly and always provides a statistically distributed surface normal orientation for any incident direction. In the year 2014 Heitz and D'Eon (2014a, 2014b) introduced a very efficient method which allows direct sampling of the distribution of visible surface normals $D_v(\vec{m})$ for V-cavity or Smith surfaces. Instead of sampling the distribution of normals first and then weighting the normals according to the masking function, the weighted distribution can directly be sampled in a single step. Three years later Eric Heitz (2017) developed the method further and proposed a new algorithm, which not only simplifies the implementation but also significantly improves the algorithm's efficiency even further. The suggested method has been implemented, and the widely used *GGX distribution of normals* has been used in this first version of the *RadiCal* raytracer. It allows to modelling of correlated and uncorrelated surface topologies defined by one isotropic ($\alpha$) or two anisotropic ($\alpha_x, \alpha_y$) dimensionless roughness parameters.

The corresponding *masking function* $G_1(\vec{v}, \vec{m})$ for the Smith GGX surface is (Heitz, 2014)

$$ 
\begin{array}{ll}
& G_1(\vec{v}, \vec{m}) = \frac{\chi^+(\vec{v} \cdot \vec{m})}{1 + \Lambda(\vec{v})} \quad \text{with} \quad \Lambda(\vec{v}) = \frac{-1 + \sqrt{1 + \frac{1}{a^2}}}{2} \\
\text{and} & a = \frac{1}{\alpha \tan \theta} \qquad \text{for isotropic roughness} \qquad (70) \\
\text{or} & a = \frac{1}{\sqrt{\cos^2 \varphi \alpha_x^2 + \sin^2 \varphi \alpha_y^2 \tan \theta}} \qquad \text{for anisotropic roughness}
\end{array}
$$

$\chi^+$ is the positive Heaviside function which returns +1 for positive values and 0 for negative arguments. It is used to remove any back-facing normal of the distribution.

For the case of uncorrelated surface roughness, the related masking-shadowing function $G_2(\vec{v}_i, \vec{v}_o, \vec{m})$ can be written as

$$ G_2(\vec{v}_i, \vec{v}_o, \vec{m}) = G_1(\vec{v}_i, \vec{m}) \cdot G_1(\vec{v}_o, \vec{m}) = \frac{\chi^+(\vec{v}_i, \cdot \vec{m}) \chi^+(\vec{v}_o, \cdot \vec{m})}{1 + \Lambda(\vec{v}_i) 1 + \Lambda(\vec{v}_o)} \quad (71) $$

For the correlated case, Heitz (2014) suggests the following relation, which considers height correlation:

$$ G_2(\vec{v}_i, \vec{v}_o, \vec{m}) = \frac{\chi^+(\vec{v}_i, \cdot \vec{m}) \chi^+(\vec{v}_o, \cdot \vec{m})}{1 + \Lambda(\vec{v}_i) + \Lambda(\vec{v}_o)} \quad (72) $$

As mentioned above, the shadowing-function is explicitly in the *RadiCal* implementation. While the shadowing-function $G_s$ matches the masking function $G_1$ in the uncorrelated case, equation (69) is used to derive $G_s$ of $G_1$ and $G_2$.

The distribution of visible surface normals $D_v(\vec{m})$ has been implemented based on the very efficient sampling algorithm provided by Heitz and D'Eon (2014a) and subsequently improved by Heitz (2017). Its implementation is provided in the next section.
4.10.2. Implementation

The functions relevant to modelling the microfacets are contained in the core of the raytracer. In
order to determine the orientation of the relevant microfacet, i.e. to sample the distribution of visible
surface normal according to introduced algorithm (Heitz and D'Eon, 2014b; Heitz, 2017), the function
sampleGGXVNDF is used. The function argument is the direction vector of the incident ray, while the
function result returns the surface normal vector.

function sampleGGXVNDF(const Ve: TVect3d): TVect3d;
// TRANSLATION OF ORIGINAL CODE SEE HEITZ2018GGX.PDF
var Nh,Vh,T1,T2: TVect3d;
s,r,phi,p1,p2,nm: Double;
u1,u2: Double;
Begin
  // transform the view direction to the hemisphere configuration
  u1:=RANDOM;
  u2:=RANDOM;
  Vh:= Vnormalize(vect(alpha_x * Ve.x, alpha_y * Ve.y, Ve.z));

  // setup orthonormal basis
  if Vh.z<0.9999 then T1:= Vnormalize(vcrossProd(ZAxis, Vh))
    else T1:= Vnormalize(vcrossProd(YAxis, Vh));

  // parameterization of the projected area
  r := sqrt(U1);
  phi := 2.0 * pi * U2;
  p1 := r * cos phi;
  p2 := r * sin phi;
  s := 0.5 * (1.0 + Vh.z);
  p2 := (1.0 - s) * sqrt(1.0 - p1*p1) + s*p2;

  // calculate the normal in this stretched tangent space
  nm:=1.0 - p1*p1 - p2*p2;
  if nm>0 then nm:=sqrt(nm) else nm:=0;
  Nh:=(p1*T1) + (p2*T2) + (nm*Vh);

  if nh.x>1E9 then nh.x:=1E9;

  // transform the normal back to the ellipsoid configuration,
  // unstretch and normalize the normal
  if Nh.z<0 then Result:=VNormalize(vect(alpha_x * Nh.x, alpha_y * Nh.y, 0))
    else Result:=VNormalize(vect(alpha_x * Nh.x, alpha_y * Nh.y, Nh.z));
end;

The shadowing function is required to check if the outgoing ray is blocked. For the correlated surface
the function SmithGGX_ShadowingCorr is used for that purpose. It is an implementation of equation
(69) using (77) and (79). For the uncorrelated case, the masking function implemented in SmithGGX_G1
can be used directly. To increase computational performance, the tangent function call is replaced
with an equivalent vector expression.

Function SmithGGX_ShadowingCorr : Double; // get shadowing func. GS=G2/G1 for correlated case
var G2,G1 : Double;
dotNL,dotNV : Double;
tanthL,tanthV : Double;
lamL,lamV : Double;
aL,aV : double;
sigL,sigV : double;
iso : boolean;
begin
dotNL := lout.z; // m=(0,0,1) dot product with surface normal
dotNV := lin.z; // m=(0,0,1) dot product with surface normal
iso:=alpha_x=alpha_y; // isotropic roughness?

if abs(dotNL)<1E-7 then dotNL:=sign(dotNL)*1E-7;
if abs(dotNL)>0.999999 then dotNL:=sign(dotNL)*0.999999;
tanthL:=sqrt(1-dotNL*dotNL)/dotNL; // tan = (1-cosT²)^(1/2)/cosT
if iso then aL:=1/(alpha_x*tanthL) else aL:=1/(alpha_i(lout)*tanthL);
if aL>=0 then sigL:=1 else sigL:-1;
lamL:=(-1+sigL*sqrt(1+(1/(aL*aL)))/2);

if abs(dotNV)<1E-7 then dotNV:=sign(dotNV)*1E-7;
if abs(dotNV)>0.999999 then dotNV:=sign(dotNV)*0.999999;
tanthV:=sqrt(1-dotNV*dotNV)/dotNV; // tan = (1-cosT²)^(1/2)/cosT
if iso then aV:=1/(alpha_x*tanthV) else aV:=1/(alpha_i(lin)*tanthV);
```matlab
if aV>=0 then sigV:=1 else sigV:-1;
lamV:=(-1+sigV*sqrt(1+(1/(aV*aV)))/2;

Result:=(1+lamV)/(1+lamV+lamL);
end;

function TLSISRoughPol.SmithGGX_G1 : Double;
// for anisotropic cases calculate effective alpha according to Heitz
var dotNV,denom : Double;
    lam : Double;
    tanth : Double;
    a : double;
    sig : Double;
begin
dotNV:=lout.z;

tanth:=sqrt(1-dotNV*dotNV)/dotNV; // tan = (1-cosT²)^(1/2)/cosT
if alpha_x=alpha_y then a:=1/(alpha_x*tanth) else
    a:=1/(alpha_i(lout)*tanth);
if a>=0 then sig:=1 else sig:-1;
lam:=(-1+sig*sqrt(1+(1/(a*a)))/2;

Result:=1/(1+lam);
end;
```

### 4.10.3. Testing and validation

In order to test the implemented algorithm, the reflection and transmission profile of float-glass and silver surfaces with different roughness values were analysed. Unpolarised light with a wavelength of 500 nm was used for incidence. The implementation and the calculation results were also briefly discussed with Eric Heitz, the author of the relevant publications.

The following polar plot representations show the reflected and transmitted angularly resolved radiance values. To be able to present the results within a single chart, the square root of the radiance values is used as the radius parameter for the polar plots. For the polar plots, only the correlated model was applied.

Figure 43 depicts the scattering of light for three different angles of incidence (grey: 10°, violet: 45°, red: 80°) and three different roughness values (left: 0.1, centre: 0.3, right: 1.0). It can be seen that most of the energy is transmitted below the surface. The surface roughness leads to a broadening of the transmitted and refracted beam, increasing with surface roughness. Likewise, the reflection on the surface is broadened/diffused. Consistent with the Fresnel equations, the reflectance generally rises with a higher angle of incidence. However, the effect is mitigated by a high roughness of the surface.

<figure><img src="image_49.png" /><figcaption>Figure 43 BSDF functions resulting of implemented microfacet roughness model for float glass surface with different roughness

[left: 0.1, centre: 0.3, right: 1.0] and different incidence angles [grey: 10°, violet: 45°, red: 80°]

(polar plot radiuses are scaled with square root of intensity)</figcaption></figure>
Figure 44 provides a closer view of the centre case of Figure 43. It can be seen that most of the incoming light is transmitted and refracted below the surface. While Snell's law determines the angle of refraction, the transmitted beam is broadened due to the variation of the surface normal. It can also be seen that the occurrence of secondary surface bounces adds a little “collar” to the lobe of the transmitted beam. This indicates that the secondary bounces enable an additional and even broader diffusion of the transmitted light into regions that would, due to the law of refraction, not be accessible on the first scattering process. However, from a quantitative point of view, the diffusion by multiple scattering events is very limited.

Analysing the reflectance of the rough silver surfaces, a generally higher reflectance and, therefore, a weaker dependence of the reflectance on incidence angle and roughness is observed. The significant effect of the surface roughness is that the reflectance loses its directional component, as can clearly be seen in Figure 45. The roughest case (right) does not show any specular reflection. The reflected light is equally diffused in all directions.

<figure><img src="image_50.png" /><figcaption>Figure 44 BSDF functions resulting of implemented microfacet roughness model for float glass surface with roughness 0.3 and different incidence angles</figcaption></figure>

[red: 80°, violet: 45°, grey: 10°]

(polar plot radiuses are scaled with square root of intensity)

<figure><img src="image_51.png" /><figcaption>Figure 45 BSDF functions resulting of implemented microfacet roughness model for silver glass surface with different roughness [left: 0.1, centre: 0.3, right: 1.0]</figcaption></figure>

and different incidence angles [grey: 10°, violet: 45°, red: 80°]

(polar plot radiuses are scaled with square root of intensity)
<figure><img src="image_52.png" /><figcaption>**Figure 46** BSDF functions resulting of implemented microfacet roughness model for silver surface with roughness 0.3 and different incidence angles
[left: 80°, centre: 45°, right: 10°]
(polar plot radiuses are scaled with square root of intensity)</figcaption></figure>

In order to analyse the difference between
the correlated roughness model vs. the
uncorrelated approach, the reflection on a
silver surface with a roughness value of
1.0 (at an incidence angle of 45°) is
analysed using both models. As shown in
Figure 47, the consideration of
correlation, on the one hand, reduces the
intensity of the reflection towards the
surface-normal, while the radiance in
directions more parallel to the surfaces is
increased. The uncorrelated model results
in an ideal Lambertian reflection profile.
The polar plot representation shows a
perfect spherelike structure, as can be
expected for a Lambertian surface.
However, the root-square scaling
approach in the presented figure slightly
distorts the sphere.

<figure><img src="image_53.png" /><figcaption>Figure 47 BSDF function silver surface with roughness 1.0 and incidence angle 45° comparison of correlated (red) vs. uncorrelated roughness model (black)</figcaption></figure>
In order to quantify the simulation results the reflectance was integrated over the entire hemisphere to determine the spectral directional-hemispheric reflectance (SDHR). Table 8 provides the total reflectance values for all cases depicted above.

<table><caption>Table 8 Spectral directional-hemispheric reflectance values for test cases</caption><thead><tr><th rowspan="3">case</th><th colspan="3">float-glass surface</th><th colspan="3">silver surface</th></tr><tr><th colspan="3">roughness</th><th colspan="3">roughness</th></tr><tr><th>0,1</th><th>0,3</th><th>1</th><th>0,1</th><th>0,3</th><th>1</th></tr></thead><tbody><tr><td>Incidence angle 10° - uncorr.</td><td>4,4%</td><td>4,1%</td><td>2,2%</td><td>98,1%</td><td>97,9%</td><td>96,0%</td></tr><tr><td>Incidence angle 45° - uncorr.</td><td>5,5%</td><td>5,1%</td><td>2,7%</td><td>98,1%</td><td>97,8%</td><td>96,2%</td></tr><tr><td>Incidence angle 80° - uncorr.</td><td>28,0%</td><td>13,8%</td><td>5,4%</td><td>98,2%</td><td>97,8%</td><td>96,9%</td></tr><tr><td>Incidence angle 10° - corr.</td><td>4,4%</td><td>4,1%</td><td>2,2%</td><td>98,1%</td><td>97,9%</td><td>96,3%</td></tr><tr><td>Incidence angle 45° - corr.</td><td>5,5%</td><td>5,2%</td><td>2,8%</td><td>98,1%</td><td>97,8%</td><td>96,6%</td></tr><tr><td>Incidence angle 80° - corr.</td><td>28,2%</td><td>14,3%</td><td>5,9%</td><td>98,2%</td><td>97,9%</td><td>97,5%</td></tr></tbody></table>

Additionally, the SDHR values were determined as a function of the roughness parameter $\alpha$ for both surface materials and the six incidence cases each (see Figure 48).

<figure><img src="image_54.png" /><figcaption>Figure 48 Spectral dir.-hemi. reflectance for silver and float glass surface as a function of surface roughness $\alpha$</figcaption></figure>

It can be seen that the reflectance of the surface is generally reduced with increasing roughness. This effect is particularly pronounced for the reflection off the glass with the high incidence angle (80°). In this case, the reflectance reduces from approx. 40% for the flat surface to almost 5% for the rough surface ($\alpha$ = 1). The results further show that choosing either the correlated or uncorrelated surface model for the float-glass has little effect on the generally low reflectance.

On the highly reflective silver surface, the reflectance decrease with roughness is generally lower, albeit clearly visible. Choosing the correlated model leads to apparently higher reflectance values; however, the absolute reflectance increase is only 0.5% at most. The reason for this increase can be identified if the average number of bounces on the surface structure is evaluated (not listed here): for the rough surface ($\alpha$ = 1) and incidence angle 80° the average number of bounces for each ray is lowered from 1.66 to 1.36 if the correlated model is chosen. This indicates that the different distribution of microfacet orientations allows the ray of light to exit more easily, reducing potential absorption.
Another less quantitative but powerful test for the validity of the implementation is the generation of rendered images, as the eye is very sensitive to inaccurate representations or singular errors. For this purpose, the raytracing algorithm is used in backward mode (see section 7.11). As can be seen in Figure 49, the microfacet algorithm provides a realistic appearance for transparent and opaque materials. It can also be seen that a roughness value $\alpha$ of 0.3 already provides a pronounced, "sandblasted" appearance of the surfaces. This indicates that higher values for $\alpha$ are less relevant for realistic modelling of surfaces.

<figure><img src="image_55.png" /><figcaption>Figure 49 Various surface roughnesses rendered with RadiCal algorithm
top row: float glass; bottom row: silver;
roughness $\alpha$ values: left column: 0.0, centre column: 0.1, right column: 0.3</figcaption></figure>
# 4.11. Subsurface scattering

## 4.11.1. Background and theory

In many rendering approaches, a three-component model is used to approximate the optical behaviour of surfaces. It involves specular reflection, diffuse surface reflection and diffuse subsurface reflection (see Figure 50).

Specular reflection directly results from solving Fresnel's equations for flat surfaces. Diffuse reflection on the surface occurs due to its topology, i.e. roughness. In the *RadiCal* approach, microfacet models are used to model this effect. Therefore, the specular (A) and diffuse surface (B) reflection on the surface are implicit parts of the extended Fresnel model, and the relevant dependence on roughness and incidence angle arises intrinsically. This approach also allows considering the presence of thin-films on the surface.

<figure><img src="image_56.png" /><figcaption>*Figure 50 Reflection types currently covered by the RadiCal implementation*</figcaption></figure>

Subsurface reflection (case C), however, has to be covered by an additional model. Subsurface reflection results from light rays penetrating through the (rough or flat) surface. If the material below the surface is a conglomerate showing a spatial variation of refractive indices, then the incoming light will suffer multiple scattering events, and finally, a fraction will exit beyond the surface again. The effect is, therefore, also referred to as ‘subsurface scattering’ or ‘volume scattering’. Materials with significant subsurface scattering are sometimes called translucent to express that they are neither transparent nor opaque. Common materials with strong subsurface scattering are wax, milk, skin, leaves or minerals. These materials allow relatively deep penetration of light below the surface. However, in many materials, subsurface scattering occurs directly below the surface.

For the indented technical application of the raytracer, paints represent the probably most relevant type of material requiring the modelling of subsurface scattering. Paints are made of a binder containing small particles referred to as pigments. Binder materials, e.g. oils, resins or gums, generally are transparent and have a low refractive index to allow the light to penetrate below the surface. In contrast, pigments exhibit very high refractivity in some parts of the spectrum and effectively scatter light of the corresponding wavelengths while other parts of the spectral range are absorbed.

The Monte Carlo raytracing approach applied in the *RadiCal* method allows a detailed simulation of multiple scattering events within the volume of an object. This enables the implementation of specific models for a variety of layered or homogeneous translucent materials. In this first version, however, a simpler approach has been implemented.

Based on a given spectrally resolved subsurface reflectance function $\rho_{ss}(\lambda)$, it is assumed that a fraction of light refracted below the surface is diffusely reflected to the outside again. This first implementation assumes a perfect diffuse directional distribution according to a Lambert reflector (see section 4.3.3). Further, it is assumed that the multiple scattering events lead to a full depolarisation of the reflected light. The two assumptions are approximations that are often applied. The
implementation of the ideal diffuse and depolarised subsurface reflection, therefore, allows modelling surfaces in line with many existing approaches, including the current ISO standard approach for modelling shading devices.

Additionally, more refined subsurface scattering models can be integrated in the future. The potential extension of the scattering models can be divided into two groups. On the one hand, commonly used BRDF scattering models, like the widely used Oren-Nayar model (Oren and Nayar, 1994), can be added. However, it has to be considered that many effects that these models are based on are already covered by the microfacet approach. On the other hand, following the raytracing approach, new subsurface models could comprise explicit subsurface scattering models and consider effects occurring when the potentially rough surface is passed again from the backside. This approach would also allow tracking the reflected light's polarisation state.

### 4.11.1.1. Sampling the outgoing direction on a Lambert surface

The subsurface reflectance module needs to provide an outgoing direction for the reflected ray following the principles of Lambert reflectance. While the BRDF function of a Lambert surface, being constant throughout the hemisphere (see section 4.3.3), might indicate uniform sampling of the hemisphere, sampling from the perspective of radiance must be considered. Following the Lambert assumption, it is evident that the directions on the hemisphere must be sampled with a cosine weighting, as directions near the surface normal will contribute stronger.

In order to obtain this distribution for the direction of the diffusely reflected light, the inversion method (see section 5.1.2) is applied to determine the z-coordinate. The PDF that has to be proportional to $cos \theta$ is integrated and subsequently inverted to obtain $CDF^{-1}$. The x and y coordinates are determined by equally distributed points on the corresponding circle of the hemisphere. Two independent random numbers $u_1$ and $u_2$ are used in this process.

$$PDF(z) = cos \theta = z \rightarrow CDF(z) = z^2 \rightarrow z = CDF^{-1}(u) = \sqrt{u} \quad (73)$$

$$ 
\begin{aligned}
z &= \sqrt{u_1} \\
x &= \sqrt{1 - z^2} \sin(u_2 \cdot 2\pi) \\
y &= \sqrt{1 - z^2} \cos(u_2 \cdot 2\pi)
\end{aligned} 
\quad (74) $$

The actual implementation follows this method; however, some numerical optimisations are applied to increase the algorithm's efficiency (see next section).

### 4.11.2. Implementation

In analogy to the refractive index module, an abstract base class for the subsurface scattering function, which provides the subsurface reflection as a function of wavelength, is defined and inherited to child classes, which carry the specific models. Currently, three different models are available; additional models can be implemented straightforwardly.

```python
TsubSurfaceRhoBaseClass = class(TObject) { TsubSurfaceRhoBaseClass }
private
    ssRhoName : String;
    cacheV : TRhoCacheV;
public
    constructor Create(const ssRhoName_ : String);
    function getRho(const w1 : Double; virtual; abstract;
                      destructor Destroy; override;
End;

TssRhoConst = class(TsubSurfaceRhoBaseClass) { TssRhoConst }
private
    rhoConst : Double;
public
```
```typescript
procedure initConst(const rhoConst_ : Double);
function getRho(const w1 : Double) : Double; override;
decorator Destroy; override;
End;

TssRhoSpline = class(TsubSurfaceRhoBaseClass) { TssRhoSpline }
private
  NPoints : Integer; // actually high value not length
  w1Point : array [0..NpointsMaxRhospline] of Double;
  rhoPoint : array [0..NpointsMaxRhospline] of Double;
  rhoPointD2 : array [0..NpointsMaxRhospline] of Double;
public
  procedure initSpline(const rhoPoint_,w1Point_ : array of Double);
  function getRho(const w1 : Double) : Double; override;
  destructor Destroy; override;
End;

TssRhoTableLinIP = class(TsubSurfaceRhoBaseClass) { TssRhoTableLinIP }
private
  minWL,maxWL : Double;
  wrPoints : array of TwrPoint;
public
  procedure initConst(const rhoConst_ : Double);
  function getRho(const w1 : Double) : Double; override;
  procedure readFromFile(const fn : String);
  destructor Destroy; override;
End;
```

The class `TssRhoConst` allows the definition of a constant subsurface reflectance value for the entire wavelength range (grey or white material). The class `TssRhoSpline` can be used to define reflectance values based on spline functions with a variable number of knots. Currently, mainly this class is used as it corresponds to the reflectance behaviour of pigments, which usually show smooth and continuous variation rather than edges. Finally, to be able to import tabulated data, the `TssRhoTableLinIP` is available. It allows describing the reflectance curve using linear interpolation between given points. In order to generate the outgoing direction, a unit vector aligned to the local system (where the z-axis equals the macroscopic surface normal) is generated. The implementation follows equation (74); however, to avoid the time-consuming sine and cosine function calls the random points on the circle, providing the x- and y-coordinates are generated using a rejection sampling method. The approach is based on an algorithm commonly referred to as ‘Malley’s method’ (Malley, 1988). The implementation was found to be 44% faster than the original algorithm. In the algorithm, firstly, a random point within the unit square is chosen. If this point lies within the unit circle, its coordinates are scaled to match the radius of the unit sphere; otherwise, a new point is processed. This method provides points that are equally distributed on a circle.

```typescript
Function LocalLambert: TVect3d;
var v : Double;
    nrm : Double;
    eta : Double;

begin
  eta:=random;
  result.z:=sqrt(eta); // cos theta distributed z coordinate, CDF-1

  (** SLOW VERSION ang:=random*2*pi;
    v:=sqrt(1-sqr(result.z));result.x:=v*sin ang;result.y:=v*cos ang; **)

  // 44% FASTER VERSION, rejection sampling of circle, avoiding tri functions
  v:=sqrt(1-eta);
  repeat
    result.x:=2*random-1;
    result.y:=2*random-1;
    nrm:=sqr(result.x)+sqr(result.y);
    until nrm<=1;

    nrm:=sqrt(nrm);
    v:=v/nrm;
    result.x:=result.x*v;
    result.y:=result.y*v;
end;
```
### 4.11.3. Testing and validation

In order to test the LocalLambert function, some basic statistical tests are performed on the distribution of the directions provided. The function is called for one billion times generating samples to perform two statistical tests. Firstly, the mean vector of all samples $|\vec{v}_i|$ is calculated. Additionally, the result is segmented into one million subsamples of sample size 1000 to analyse the validity and convergence of the three coordinates of the mean vector.

The calculation results for the average coordinates and standard deviation are:

$$|\vec{v}_i| = \begin{pmatrix} 0.000015 \\ -0.000014 \\ 0.666668 \end{pmatrix} \qquad \begin{pmatrix} \sigma_x \\ \sigma_y \\ \sigma_z \end{pmatrix} = \begin{pmatrix} 0.015812 \\ 0.015820 \\ 0.007455 \end{pmatrix} \quad \text{for } n = 1\,000\,000, m = 1000$$

For the cosine-weighted hemisphere, the mean value is expected at $(0,0,2/3)$, which can be shown by integration. In addition to that, the central limit theorem (see section 5.1.3) is applied to check for an error or bias in the result. For this purpose, the statistical error that has to be expected for the given number of samples is determined. Applying equation (91), the standard error of means for the three coordinates are:

$$\sigma_{SEM} = \begin{pmatrix} 0.000016 \\ 0.000016 \\ 0.000007 \end{pmatrix}$$

Therefore, the deviations of $|\vec{v}_i|$ are well within the range that can be expected statistically. No bias or error in the results is detected.

For the second test, a histogram of the results is created to validate the cosine dependence of the directions. If valid, the density of the random directional samples has to increase with the cosine of the angle between the given direction and the surface normal vector. Since the generated directions are aligned to the local surface plane with $\vec{n} = (0,0,1)$ and the angle between the direction vector $\vec{v}$ and $\vec{n}$ can be written as $\vec{v} \cdot \vec{n} = (0,0,v_z) = \cos \alpha$, a linear relationship between the density of the random directions and their z coordinate is expected. In order to test this, a billion directional samples were compiled in a histogram with 100 bins using the z-coordinates for classification.

<figure><img src="image_57.png" /><figcaption>Figure 51 Frequency/sample density of randomly generated Lambert directions as a function of z-coordinates</figcaption></figure>

<figure><img src="image_58.png" /><figcaption>Figure 52 Frequency/sample density of randomly generated Lambert directions as a function of angle to surface normal</figcaption></figure>

As shown in Figure 51 and Figure 52, the almost perfect correlation of the distribution of sampled directions with the z-coordinates confirms the cosine dependence and, therefore, the validity of the implemented Lambert algorithm.
# 4.12. Inversion method to derive complex-valued refractive index functions of transparent materials

## 4.12.1. Background and theory

As stated above, it is essential to have knowledge of the complex-valued refractive index (cRI) functions of all materials used in physically-based raytracing. While it is possible to retrieve models or empirical data for some “bulk” materials, e.g. of the *refractiveIndex* database (Polyanskiy, 2022), this data is unfortunately currently not available for many other materials. For building science applications, it is of particular importance to have models for coated and uncoated glass panes, to be able to model glazings accurately. While the necessary refractive index information is not available, Lawrence Berkley National Laboratories kindly maintains the freely accessible and continuously updated glazing database *IGDB* (LBNL Windows and Daylighting Group, 2021), containing spectral measurements for almost all relevant types of glass. The database provides transmittance and reflectance spectra for the global radiation range. The datasets typically include measured spectral functions for the frontside reflectance $R_f$, backside reflectance $R_b$ and transmittance $T$. These measurements are generally provided for normal, i.e. near-normal, incidence only. The data is provided for uncoated as well as for coated glass panes.

A method has been developed to generate the required cRI functions based on these measured spectral data by inversion, i.e. appropriate parametric models for the cRI functions are varied until the modelled reflectance and transmittance spectra match empirical spectra provided by *IGDB*. The two real-valued components forming the complex-valued refractive index function, the refractive index $n(\lambda)$ and the extinction coefficient $\kappa(\lambda)$, have to be determined for the entire global radiation spectrum. Due to the many resulting degrees of freedom, the inversion can only be performed using a multi-parameter optimisation algorithm.

### 4.12.1.1. Thick slab model

While the MC-raytracing approach could be used to determine the total reflectance and transmittance of the uncoated or coated glass panes, the optimisation algorithm to fit the refractive indices relies on fast computation. Since the total reflectance and transmittance of the glass panes, being the result of interreflections, can be solved analytically, there is no reason to apply the – for this application slower – Monte-Carlo method. Instead, the thick

<figure><img src="image_59.png" /><figcaption>Figure 53 Applied "thick-slab" model</figcaption></figure>

slab model is applied to efficiently model a coated or uncoated glass pane's front- and backside reflectances and transmittance. As discussed in section 4.3.7, this model assumes a non-coherent superposition of many interreflections within the glass pane. Figure 53 provides a schematic overview of the model. Since the measurements are performed with air contact, the media 0 and 2 are modelled as air (approximated as vacuum). Medium 1 is the actual glass material. Interfaces 01 and 12 are modelled using the generalized Fresnel approach (see section 4.7). At interface 01, a thin-film stack
consisting of N layers can additionally be considered. The interfaces' reflectivity and transmissivity values are denoted with ordered indices to indicate the direction, even though most are identical for both directions. All coefficients are determined using equations (53). Note that for the case of a thin-film stack present at interface 01, any absorption within the thin-film will also be included in the model, as for this case $T_{01} + R_{01} < 1$ will apply.

Beyond the Fresnel coefficients, relevant for the scattering processes at the interfaces, the transmittance $T_1$ of the glass pane is required. It is obtained by applying Lambert-Beer's law (see section 4.3.6):

$$T_1 = \frac{I_{12}}{I_{01}} = \frac{I_{01} \cdot e^{-\alpha l}}{I_{01}} = e^{-\alpha \frac{d_1}{\cos\theta_1}} \quad (75)$$

$I_{01}, I_{12}$ ... irradiance at top and bottom of glass pane [W/m²]

l ... optical path length [m]

$d_1$ ... thickness of glass pane [m]

$\theta_1$ ... propagation angle in glass pane [rad]

$\alpha$ ... attenuation coefficient [1/m]

Whereas the attenuation coefficient $\alpha$ is calculated using equation (40), and the angle $\theta_1$ is calculated based on Snell's law (46).

The derivation of the total reflectance and transmittance, being the superposition of all interreflections, can easily be comprehended if the coefficients $R$ and $T$ are understood as probabilities. In this view, the probabilities for each path have to be added, while the required probabilities forming a path have to be multiplied. A similar derivation, but with slightly different components, is, e.g., provided by Stenzel (2016, p. 131).

The infinite summation of all reflection paths (see Figure 53) then leads to a geometric series, which can easily be evaluated into a simple fractional term:

$$ 
\begin{aligned}
T_{pane} &= T_{01}T_{1}T_{12} + T_{01}T_{1}R_{12}T_{1}R_{10}T_{1}T_{12} + T_{01}T_{1}R_{12}T_{1}R_{10}T_{1}R_{12}T_{1}R_{10}T_{1}T_{12} + \cdots = \\
&= T_{01}T_{1}T_{12} \cdot \left[ \sum_{n=0}^{\infty} (R_{10}T_{1}^{2}R_{12})^{n} \right] \quad (76) \\
T_{pane} &= \frac{T_{01}T_{1}T_{12}}{1 - R_{10}T_{1}^{2}R_{12}}
\end{aligned}
$$

In complete analogy, the total reflectance for the front side and, by flipping the direction and exchanging the interface indices correspondingly, the backside of the glass pane can be derived. The frontside reflectance $R_{f,pane}$ and the backside reflectance $R_{b,pane}$ are:

$$R_{f,pane} = R_{01} + \frac{T_{01}T_{1}^{2}R_{12}T_{10}}{1 - R_{10}T_{1}^{2}R_{12}} \quad (77)$$

$$R_{b,pane} = R_{21} + \frac{T_{21}T_{1}^{2}R_{10}T_{12}}{1 - R_{12}T_{1}^{2}R_{10}} \quad (78)$$

Note that the expressions could further be simplified by expressing transmittance values by reflectivities and applying identities for transmittance values in reverse directions. However, since the coefficients are readily provided by the algorithms and in order to average out numerical inaccuracies, the expressions were implemented as stated in equations (76) and (77). Note also that all quantities
are functions of wavelength and incidence angle; for reasons of clarity, the corresponding arguments were omitted in the above equations.

### 4.12.1.2. A parametric model for coated and uncoated glass panes

After some testing, the following models were identified as best suited to determine the cRI functions by inversion. As is always the case for such optimisation problems, it is critical to provide a physically appropriate model with as many variable parameters as required - but only as few as necessary, to avoid overtraining the model.

In the implemented method, mainly natural cubic spline functions (G., Press et al., 1991) with N knots are used to individually approximate the real-valued functions $n(\lambda)$ and $\kappa(\lambda)$ that constitute the complex-valued refractive index function. Since $\kappa(\lambda)$ often ranges over several orders of magnitude, a logarithmic mapping of the spline was tested but was found not to perform better than the linear approach.

The following constraints were applied in the optimisation process:

* the wavelength ("x-coordinate") of the first knot of the spline was constrained to range within 0 nm and 300 nm.

* the wavelength ("x-coordinate") of the last knot of the spline was constrained to range within 2500 nm and 3000 nm.

* all $\kappa$ values for glazings were constrained in the range of [1E-8, 10].

* all $\kappa$ values for coatings were constrained in the range of [1E-8, 100].

* all $n$ values for glazings were constrained in the range of [1E-5, 4].

* all $n$ values for coatings were constrained in the range of [1E-5, 6].

Note that the actual $n$ and $\kappa$ functions for the coatings do not have to reflect the true physical properties if taken individually and, therefore, can seem entirely out of range. A two-layer substitute model is applied to approximate the behaviour of the thin-film stack, potentially consisting of a double-digit number of different materials, e.g. for a low-E coating (Ding and Clavero, 2017). The actual optical behaviour of the thin-film assembly is determined as the result of the multiplications of the individual characteristic matrices of each layer (see equations (50) - (52)). This matrix multiplication relies on the cRI functions of the individual layers, as well as on their thicknesses $d_i$. Therefore, physically not meaningful values for $n$, $\kappa$ and $d$ achieved for the single layers of the two-layer model can still be suitable to resemble the correct physical behaviour of a multi-layer assembly. While spline functions are used to model the refractive index and extinction coefficient functions of all coating layers, the refractive index of the glass can be modelled based on a less complex model. The following dispersion formula (see TRIsimpleDispersionAndKspline of Table 5) has proven to be suitable to model the relevant uncoated glass pane types:

$$n = C_1 + C_2\lambda^2 + \frac{C_3}{\lambda^2} \qquad (79)$$

As can be seen, this function relies on three parameters only, while the natural cubic spline function has N degrees of freedom.

Since IGDB contains data for coated and uncoated glass panes, the uncoated panes' spectra are always used to determine the model parameters for the plain glass material. The given thickness of the glass pane is provided as a constant model parameter.
<figure><img src="image_60.png" /><figcaption>Figure 54 Applied parametric models for uncoated (left) and coated glass panes (right); free model parameters marked in red</figcaption></figure>

In order to determine the model for coated glass panes, the model parameters of the uncoated glass, comprising three constants for the refractive index (eq. (79)) and a series of N spline knots for the extinction coefficient, are used as constant values. Only the parameters relevant for the two coating layers, comprising the thickness and the two series of knots coordinates for each coating layer, are the free parameters in the optimisation process to replicate the measured spectra of the coated glass pane.

### 4.12.1.3. An ideal optimisation cost function

Ideally, the optimisation should be controlled by applying transmittance and reflectance spectra for multiple incidence angles as the target quantity of the optimisation. Based on the experience gained when working with the optimisations, it seems likely that a number of only two to four different incident angles should be sufficient for this purpose.

The objective of the optimisation process is defined by providing a cost function that rates the optimality of a solution. For the present case, the following definitions are proposed: first, the deviations ($\Delta T$, $\Delta R_f$ and $\Delta R_b$) of the model vs. the measurement for a specific incidence angle $\theta_n$ and the wavelength $\lambda_m$ are defined as:

$$ 
\begin{aligned}
\Delta T(\theta_n, \lambda_m) &= T(\theta_n, \lambda_m) - \hat{T}(\theta_n, \lambda_m) \\
\Delta R_f(\theta_n, \lambda_m) &= R_f(\theta_n, \lambda_m) - \hat{R}_f(\theta_n, \lambda_m) \\
\Delta R_b(\theta_n, \lambda_m) &= R_b(\theta_n, \lambda_m) - \hat{R}_b(\theta_n, \lambda_m)
\end{aligned}
\qquad (80) $$

Where $T$, $R_f$ and $R_b$ are the model results, according to (76) - (78), and the corresponding values marked with the hat accent represent measurement values provided by the IGDB spectrum files. The global radiation spectrum is sampled in two dimensions. First, by using $M$ distinct, energetically equidistant wavelengths $\lambda_m$ (see section 4.4.2), and second, for $N$ available incidence angles $\theta_n$. The loss functions to assess the quality of the fitting are defined as:

$$ 
\begin{aligned}
sqr_err &= \sqrt{\frac{\sum_{n=0}^{N-1} \sum_{m=1}^{M} [\Delta T(\theta_n, \lambda_m)^2 + \Delta R_f(\theta_n, \lambda_m)^2 + \Delta R_b(\theta_n, \lambda_m)^2]}{N \cdot M}} \\
mean_err &= \frac{\sum_{n=0}^{N-1} \sum_{m=1}^{M} [\Delta T(\theta_n, \lambda_m) + \Delta R_f(\theta_n, \lambda_m) + \Delta R_b(\theta_n, \lambda_m)]}{N \cdot M} \\
max_err &= max_{m,n}[|\Delta T(\theta_n, \lambda_m)|] + max_{m,n}[|\Delta R_f(\theta_n, \lambda_m)|] \\
&\quad + max_{m,n}[|\Delta R_b(\theta_n, \lambda_m)|]
\end{aligned}
\qquad (81) $$
The weighted loss functions finally determine the cost function for the optimisation process:

$$C = c_1 \cdot \text{sqr\_err} + c_2 \cdot \text{mean\_err} + c_3 \cdot \text{max\_err} \quad (82)$$

The coefficients $c_i$ can be used to balance the three components of the cost function. Unfortunately, this approach could not be applied and tested fully, as only spectra for one angle (near-normal incidence) were available. However, based on the presented definitions, an extended cost function is introduced below.

### 4.12.1.4. The pragmatic optimisation cost function

It turned out that applying spectral data only for one angle (near-normal incidence) was insufficient to perform a meaningful determination of the model parameters for the refractive index functions. As the model offers a relatively large number of free parameters, overfitting occurs if insufficient fitting data is provided. Therefore, additional data or constraints are needed. Since no measurements for varying angles are available, two additional constraints are imposed:

1. the normal-incidence reflectance and transmittance spectra should also be valid for a specified low, but non-zero, incidence angle $\theta_1$.

2. the ratio of the total transmittance for a specified higher oblique angle $\theta_2$ to the total transmittance at normal incidence should match a specified value.

Both constraints are physically reasonable. The first condition helps to avoid spectrally sensitive thin-film solutions that frequently occur if the constraint is not applied. Such solutions will provide the correct reflectance behaviour for the normal incidence but will show strong spectral and angular dependences for oblique angles, e.g. as the colourful “rainbow” interference known from thin oil or soap films. Such effects are, in fact, avoided in the production of glazings, as this would impair the visual quality of the glazing. Since the angular characteristics of glazing generally show very little change at incidence angles less than 10 to 20 degrees (Rubin et al., 1999; Roos et al., 2001), an angle $\theta_1$ within this range can be chosen.

To further restrict the optimisation results to physically plausible solutions, a target value for the decrease of the transmittance with higher angles is specified in the second constraint. While the entire spectrally resolved reflectance and transmittance data is provided in the first constraint, here, only the total (=spectrally integrated) value of the transmittance is considered. Since the sampling wavelengths $\lambda_m$ are provided for energetically equidistant intervals, the total transmittance can simply be calculated by averaging the values. Hence, the additional loss function can be written as:

$$ang\_dec = \left| \frac{\sum_{m=1}^{M} T(\theta_2, \lambda_n)}{\sum_{m=1}^{M} T(\theta_0, \lambda_n)} - f_T \right| \quad (83)$$

Where $f_T$ represents the target value for the decay of the total transmittance. The other loss functions remain the same as stated above, however two angles are considered there now: $\theta_0$ representing near-normal incidence, and $\theta_1$ as low angle for which the normal-incidence spectra should still be valid.

$$C = c_1 \cdot \text{sqr\_err} + c_2 \cdot \text{mean\_err} + c_3 \cdot \text{max\_err} + c_4 \cdot \text{ang\_dec} \quad (84)$$

In the practical implementation, the relevant angles and weighting-factors were defined as: $\theta_0 = 10^{-3} \circ$ (nonzero for numerical stability), $\theta_1 = 15 \circ$, $\theta_2 = 70 \circ$ and $c_1 = c_2 = c_4 = 1$, while $c_3 = 0.5$.
## 4.12.2. Implementation

The implementation is contained in the module `rc_fittrrspectrum3`, whereby the version number indicates the evolution of the module. The class `TfitTRRSpectrum` contains all functions necessary for the optimisation process and additional tools used to import, export and report the optimisation results.

The progress of the potentially multi-hour-long optimisation can be monitored in text form, as well as in the form of an updating chart.

```python
Type TfitTRRSpectrum = class(TOBJECT) { TfitTRRSpectrum }
private
    [...] local variables [...]
procedure drawChart(const layMax: Integer; const logSc: boolean);
procedure getSplineN_D2; // spline setup N
procedure getSplineK_D2(const layer0also: Boolean); // spline setup K
function getN(const nlay: IntegerR; const w1: Double): Double; // n(lambda)
function getK(const nlay: IntegerR; const w1: Double): Double; // k(lambda)
// functions for annealing process
function getCost(const calcErr: boolean): Double;
procedure alter;
procedure copyToOld;
procedure saveBest;
procedure restoreFromOld;
procedure restoreBest;
public
    fresAIR_F, fresAIR_B, fresF_AIR, fresB_AIR: TFresnel; // Fresnels of front-/backside
    spect : TIGDBSpectrum; // IGDB spectral data
    constructor Create(var memo: TMemo; var oChart: Tchart); // constructor
    procedure readSpect(const fn: String); // IGDB file import
    procedure Anneal(const testFile: String; const RedoResults: Boolean); //optim.
    procedure ReadNKPoints(const fn: String); // write opt. results
    procedure WriteNKPoints(const fn: String); // read opt. results
    procedure DumpCurve;
    procedure TestRI;
    function ParamDump: String;
    destructor Destroy; override; // destructor
end;
```

The implementation is an application of the simulated annealing algorithm (see section 5.2), where equation (84) is used as a cost function. The relations for the thick slab model derived above, see (76) - (78), are used to compute the analytic transmittance and reflectance values used for the fitting to the empirical data during the optimisation process.

In the implementation, the `ipfisn` and `ipfspn` of the `NumLib` library (Free Pascal Team, 2021) are used to solve and evaluate the linear equation systems for the applied natural cubic spline fitting.

## 4.12.3. Testing and validation

In order to demonstrate the potential of the method, the complex-valued refractive index (cRI) functions of three different materials are determined from their transmission and reflection spectra. The objective of the test is two-fold:

1.) If applied in the raytracing (or thick-slab model), the determined cRI function should, with good accuracy, result in transmission and reflection spectra matching the measured data they were generated on. The match is quantified with three different parameters that are also applied during the optimisation process:

a. Mean error of samples: as the wavelength samples are distributed according to the global radiation spectral density, the total mean error also coincides with the integrated total error considering the entire spectrum.

b. Root-mean-square error of samples: this value is the best indicator for the average spectral deviation of the three spectra.
c. Max. error of samples: the maximum deviation for any sample in any of the three spectra (transmission, reflection-frontside, reflection-backside) is determined.

2.) The determined cRI functions should be physically correct, i.e. represent the “true” cRI functions.

The second point cannot fully be validated here, as a thorough validation would require the availa-
bility of measured transmission and reflection spectra, as well as of measured cRI functions (e.g. by
applying spectroscopic ellipsometry). While this data is not available for the materials considered
here, the plausibility of the data can still be assessed by comparison with available similar cRI
functions.

**4.12.3.1. Case 1 - determination of cRI functions for uncoated window pane**

In this first case, the cRI of an uncoated pane used for insulated glazing units (IGU) is determined.
The result is later used in the full-system validation (chapter 8). The pane is of type *Planibel Clearlite*
by the manufacturer *AGC Glass Europe*. The nominal thickness is 4 mm, while the actual thickness
is 3,85 mm. The related spectra were exported from the *IGDB* database (LBNL Windows and
Daylighting Group, 2021).

For the inversion process to determine the cRI of the uncoated glass, the dispersion relation stated
in equation (79) was used to model the refractive index function (n), while a spline function with 10
knots was used to model the extinction coefficient function (κ). The knots are moved in both dimen-
sions during the optimisation process. Only the wavelength coordinates of the first and last knot are
held constant at the boundaries of the global radiation spectrum. The optimisation process comprised
35 million iteration steps and lasted a few hours. The finally determined cRI is depicted in Figure
55. The dots on the blue line indicate the knots of the spline function.

<figure><img src="image_61.png" /><figcaption>*Figure 55 Determined complex-valued refractive index function of uncoated glass pane (Planibel Clearlite)*</figcaption></figure>

<figure><img src="image_62.png" /><figcaption>*Figure 56 Empirical complex-valued refractive index function of clear-glass determined by Rubin*
*(Rubin, 1985)*</figcaption></figure>

Figure 56 shows the cRI as determined by Rubin (1985) for unspecified “clear glass”. As Rubin notes in his publication, the absorptance and reflectance of the window glass are affected by the FeO concentration and further additives. This means that significant variations can be expected among different clear glasses. However, as can be seen by comparison of Figure 55 and Figure 56, the match between the values for a modern glass, determined with the novel method presented here and the values determined by Rubin in 1985, is remarkably high.
The match of the transmission and reflection spectra is depicted in Figure 57. The precision of the reconstruction is very high, as the only visible deviation occurs in the transmission spectrum at approx. 2.2 µm. However, the deviation is still small, and the energetic impact of the deviation is even less significant, as this part of the global radiation spectrum exhibits a low spectral flux density.

The maximum sample error is 1.1% for the transmittance and less than 0.08% for the reflectance, and the root mean squared error is 0.4% for all transmittance and reflectance samples. The mean sample error is 0.09% for the transmittance and + 0.001% for the reflectance samples. Since mean sample errors are spectrally weighted, the errors are consistent with the total error for the entire global radiation spectrum.

<figure><img src="image_63.png" /><figcaption>Figure 57 Uncoated window pane - measured T/R-spectra vs. model results using the derived cRI (at near-normal incidence)</figcaption></figure>

<figure><img src="image_64.png" /><figcaption>Figure 58 Uncoated window pane - angular dep. of transmittance and reflectance based on derived cRI (integrated over the entire global radiation spectrum)</figcaption></figure>

Figure 58 finally shows the angular dependence of the transmittance and reflectance of the glass pane for global radiation. As discussed in the state-of-the-art analysis (section 3.2), even the more advanced, commonly used modelling approaches apply simplified models in order to determine the angular dependence of transmittance and reflectance. The depicted angular dependencies are relevant for global radiation determined by spectrally weighted integration.

### 4.12.3.2. Case 2 - determination of cRI functions for low-E coated window pane

In this example, a low emissivity (low-E) coating (of type iplus Top 1.1) is applied to the previously considered uncoated glass pane. As discussed in the theory section, the coating, consisting of multiple layers, is approximated using a two-layer substitute model. The already determined cRI for the 3.85 mm thick glass pane is held constant in the optimisation process. The optimisation involves the $n$ and $\kappa$ functions of both layers, where each function is modelled using spline functions with 10 knots. Beyond that, the two thicknesses $d$ of the layers are variables during the optimisation process. As discussed above, an additional parameter for the total angular weakening of transmittance is required due to the lack of angular spectral measurements. The value was taken from measurement data contained in a publication by Rubin *et al.* (1999). The transmission of the specimen 'double Ag 3mm clear' was found to decrease to 63.8% of its near-normal value at 70°. This value was used as a target for the optimisation process. The optimisation terminated after 12.3 million iterations. The layer thicknesses for the substation model were determined as 9.0 nm for the first (=glass-sided) layer and 37.6 nm for the second (=top) layer, respectively.
<figure><img src="image_65.png" /><figcaption>Figure 59 Determined cRI functions for the layers of the two-layer substitute model for the low-E coating (iplus Top 1.1)</figcaption></figure>

As mentioned in the theory section, the two
individual coating layers do not necessarily
need to have physically meaningful cRI func-
tions or thicknesses, as their effect is deter-
mined by applying the matrix method to the
entire stack. However, it can clearly be seen
that the method again captures relevant phys-
ical properties of the materials, as the cRI of
layer 1 (n1, k1 in Figure 59) shows a clear
qualitative and quantitative resemblance to
the measured values (Ciesielski *et al.*, 2017;
Polyanskiy, 2022) for a 20 nm thick silver
coating deposited on a SiO₂ substrate (see Fig-
ure 60). This is remarkable, as the optimisa-
tion process has only a few restrictions and is
offered a large number of degrees of freedom
available. While the actual layering of the low-
E coating is not available to the author, it can

be assumed that it comprises several different materials, including one to three silver layers. Due to production and durability aspects, the actual layering is usually complex and involves additional materials (Ding *et al.*, 2017).

The comparison of the measured vs. the modelled spectra (Figure 61) shows that the applied two-layer model is able to replicate the transmittance and reflectance behaviour of the more complex layering very well. This is also supported by the error measures: transmission: max. sample error: 2.8%, mean sample error: <0.01%; reflection-frontside: max. sample error: 0.85%; mean sample-error: <0.01%; reflection-backside: max. sample error: 3.2%; mean sample-error: <0.01%;).

<figure><img src="image_66.png" /><figcaption>Figure 60 Measured cRI function of 20 nm-thick silver film deposited on SiO2 substrate

source: (Ciesielski, et al., 2017; Polyanskiy, 2022)</figcaption></figure>
<figure><img src="image_67.png" /><figcaption>Figure 61 Low-E coated window pane - measured
T/Rf/Rb-spectra vs. model results
(at near-normal incidence)</figcaption></figure>

<figure><img src="image_68.png" /><figcaption>Figure 62 Low-E coated window pane - angular dep.
of transm., refl. and absorption based on derived cRI
(integrated over the entire global radiation spectrum)</figcaption></figure>

Figure 62 shows the angular dependencies for transmittance, frontside- and back-side reflectance as well as absorptance. Additionally, two references for the angular transmittance are depicted. The first is taken from Rubin *et al.* (1999) and represents measured data for a silver double-coated pane with a nominal thickness of 3 mm. Since the transmission ratio at 70° was an optimisation target, it shows a perfect match at this angle. However, a good match is also visible regarding the general shape of the curve. The second reference was taken from Roos *et al.* (2001). Figure 10 of the cited publication provides a model fit for a silver-coated Planitherm pane with unspecified thickness and few angular measurements. The angular dependence modelled by Roos *et al.* shows slightly higher values than are determined here.

4.12.3.3. Case 3 - determination of cRI functions for ETFE membrane

In order to include a non-glass example of practical relevance, the cRI function of a clear ETFE membrane, as used for skylights in buildings, is determined. Transmission and reflection spectra for the relevant global radiation range were extracted from a chart provided by a manufacturer (Buitink Technology, 2022) (“standard foil” curves). The thickness of the sample that was used to determine the spectra is specified as 200 µm.

While an appropriate dispersion relation was used to model the *n* function for the uncoated glass, the cRI of the ETFE membrane is now determined by optimising two spline functions. The number of available knots was set to 9 for the *n*-component and 19 for the *κ*-component, as more fluctuation in the absorption can be expected. Even though the optimisation is offered a very high number of degrees of freedom for the spline functions, remarkably smooth and plausible *n* and *κ* functions are determined in the optimisation process (see Figure 63). As can be seen in the *n* function, the optimisation process concentrated most knots below 1 µm to capture the rise towards shorter wavelengths accurately.
<figure><img src="image_69.png" /><figcaption>Figure 63 Determined complex-valued refractive index function of ETFE-membrane</figcaption></figure>

Again literature was reviewed to assess the validity of the inversion result. A reference containing measurement data has been found in a publication on *Optical properties of materials for concentrator photovoltaic systems* (French *et al.*, 2009). The publication includes results of spectroscopic ellipsometer and spectrophotometer measurements performed on a clear ETFE membrane.

<figure><img src="image_70.png" /><figcaption>Figure 64 Measured index of refraction and absorption coefficients of ETFE membranes (blue lines)
source: French *et al.* (2009)</figcaption></figure>

The comparison of the measured *n*-function (blue curve in Figure 64 left) vs. the function determined by inversion (Figure 63) shows that for wavelengths longer than 1,3 µm, both functions quickly converge to an almost constant value around 1.40. In the UV range, a sharp rise towards shorter wavelengths, typical for polymers, is visible. While the shape of the rise is very similar, a quantitative deviation is clearly visible in this region. However, the products are of different manufacturers, and limitations of the measurement devices in the UV region can distort the measurement results at the end of the measurement range.

The limited vertical resolution of the absorption coefficient chart (blue curve in Figure 64 right) allows only a crude assessment of the validity of the κ-function, containing the extinction coefficients. As can be seen in the chart, three main absorption regions can be identified (indicated with blue arrows). These correspond to three absorption regions determined in the inversion, see Figure 63. Note that Figure 64 uses a linear scale, while the extinction coefficient in Figure 63 is depicted using a logarithmic scale. By converting the estimated absorption values using Beer-Lamberts law, eq.
(26), it can be shown that the values match, at least in order of magnitude, being the relevant scope here.

Figure 65 demonstrates that the match between the measured and modelled transmission and absorption spectra (the optimisation process's target) is excellent: transmission: max. sample error: 1.8%, mean sample error: <0.01%; reflection: max. sample error: 0.4%; mean sample-error: <0.01%. The derived angular transmittance and reflectance for global radiation are depicted in Figure 66.

<figure><img src="image_71.png" /><figcaption>Figure 65 ETFE membrane - measured T/R-spectra
vs. model results using the derived cRI
(at near-normal incidence)</figcaption></figure>

<figure><img src="image_72.png" /><figcaption>Figure 66 ETFE-membrane - angular dependence of transmittance and reflectance based on derived cRI (integrated over the entire global radiation spectrum)</figcaption></figure>
## 4.13. Inversion method for a generic model of diffusely reflecting, opaque surfaces

### 4.13.1. Background and theory

The following method is proposed to derive generic surface reflectance models for painted or coated opaque surfaces. The approach can be considered an intermediate model downgrade and follows the principle of fallback approaches (see section 2.2). While simplified assumptions are made for the complex-valued refractive index, the diffuse reflection based on subsurface scattering still contains spectrally resolved information. Beyond that, even the assumption of a constant refractive index is still able to model scattering events in a spectrally resolved way.

Highly resolved spectral measurement data is often unavailable for practically relevant, coated or painted materials, e.g. lammelas of shading elements. Nonetheless, the manufacturer will usually still provide integrated reflectance values for specific ranges of the global radiation spectrum. The reflectance values usually represent either the total diffuse or the total hemispherical reflectance, i.e. they exclude or include specular reflection. The values are usually determined by using incidence light with a low, near-normal direction. The commonly applied standard ASTM E 903 (ASTM E 903, 1996) advises using incidence angles in the 6 to 12 degrees range. In order to determine the solar reflectance value, the reflectance values are weighted based on a global radiation standard (see section 4.4).

As discussed in section 4.3.3, the diffuse reflection of painted or coated material is commonly the result of subsurface scattering processes on pigments submerged in a binder. It is, therefore, a reasonable approximation to apply a generic model for these materials, assuming a constant or smooth refractive index for the binder with little or no absorption. The subsurface scattering on the pigments, providing the characteristic absorption and colour of the surface, is then only modelled by using the subsurface reflectance function $ρ_{ss}(λ)$, see section 4.11.

The actual shape of the spectral reflectance function $ρ_{ss}(λ)$ is determined in an optimisation process. Again, simulated annealing (section 5.2) is utilized for this purpose. In the optimisation process, the reflectance measurement is performed virtually and in accordance with the relevant standard (ASTM E 903, 1996). The spectral distribution of the incidence light is modelled according to the standard the measurement data refers to. The intensity and spectral distribution of the reflected rays are tracked, and the function $ρ_{ss}(λ)$ is altered until the quantities measured in the model match the provided empirical data. The optimisation is performed by moving the knots of the implemented cubic-spline subsurface reflectance function (see section 4.11). The use of the spline functions to approximate the reflectance spectra takes into account that solid-state absorption (unlike gas absorption) usually exhibits a continuous and smooth spectral dependence.

The empirical reflectance data can be provided for hemispherical reflectance or for diffuse reflectance. In the latter case, a filter is used to exclude the specular reflection (gloss) from the measurement. The same features are implemented in the method to ensure that the model measurement mimics the real-world measurement.

As the empirical data commonly only covers four different reflectance values (UV, visible, NIR and total) and the cubic spline has comparably many degrees of freedom, an additional optimisation objective is generated by setting a target for the colour of the surface. For this purpose, the spectrum of all diffusely reflected rays is converted into an RGB triple. This is performed by first converting the wavelengths into the XYZ colour space of CIE 1931 and then further into the RGB colour space
(Fairman *et al.*, 1997; Service, 2016; Wikipedia, 2022). Since the brightness of the colour is already determined by the visible reflectance value, only the normalised colour value is defined as the optimisation target. A third, less significant optimisation target measures the variability of the reflectance functions. It is simply defined as the average vertical distance between successive spline knots.

The three optimisation targets and an additional precision term are formulated as loss functions to assess the quality of the fitting.

They are defined as:

$$ 
\begin{aligned}
R_{\text{err\_rmsq}} &= \sqrt{\frac{(\Delta R_{\text{UV}})^2 + 2 \cdot (\Delta R_{\text{VIS}})^2 + (\Delta R_{\text{IR}})^2 + 2 \cdot (\Delta R_{\text{TOTAL}})^2}{6}} \\
col_{\text{err}} &= \sqrt{\frac{(\Delta R_{\text{red}})^2 + (\Delta R_{\text{green}})^2 + (\Delta R_{\text{blue}})^2}{3}} \\
vari\_penalty &= \frac{\sum_{i=1}^{N} (R_{\text{spline.i}} - R_{\text{spline.i-1}})}{N} \\
prec\_penalty &= 3 * curPrec
\end{aligned}
\quad (85) $$

$$ \Delta R_x = \Delta R_{x,\text{model}} - \Delta R_{x,\text{measured}} $$

$$ R_{\text{spline.i}} = y\text{-coordinante of i-th spline knot} $$

$$ N \text{ ... maximum index of spline knots} $$

The weighted loss functions determine the entire cost function for the optimisation process:

$$ C = c_1 \cdot R_{\text{err\_rmsq}} + c_2 \cdot col_{\text{err}} + c_3 \cdot vari + c_4 \cdot prec_{\text{penalty}} \quad (86) $$

$R_{\text{err\_rmsq}}$ represents the deviation for the reflectance in the given spectral ranges, $col_{\text{err}}$ is the deviation of the normalised colour value and $vari\_penalty$ provides a penalty for strongly fluctuating functions. Finally, $prec\_penalty$ adds a penalty related to the precision of the current solution. In order to increase the performance of the optimisation, the first optimisation cycles are performed with fewer samples. The potential statistical error of the results is considered by adding three times the standard error to the cost function. This very effectively prevents early solutions from being identified as optimal solutions. The coefficients $c_i$ are used to balance the four components of the cost function. In the optimisations performed, the applied cost function weighting was: $c_1 = c_2 = c_4 = 1$ and $c_3 = 1/1000$.

The method has proven to provide stable results matching the targeted reflectance quantities with good accuracy. However, given the limited restrictions for the optimisation process and the relatively high degree of freedom, the simulation results can certainly not be considered unique. In general, multiple different solutions will be able to meet the optimisation targets. Even if the generated subsurface reflection includes a certain degree of randomness and deviates from the actual spectral reflectance, it will still be a good approximation for the unknown function as it exhibits the same integrated reflectance values within the specific, provided spectral ranges and is based on a continuous reflectance function. In commonly applied alternative methods, a constant diffuse reflection coefficient $\rho_{sol}$, which is integrated over the entire global radiation range, is generally used.
### 4.13.2. Implementation

The implementation is contained in the module `rc_TsubsurfaceOpt`. The method `OPT_ASTM_E903_Reflectance` can be called to perform an optimisation of subsurface reflectance function represented by a cubic spline with N Knots.

```c
function TsubsurfOpt.OPT_ASTM_E903_Reflectance (var fmat: TRCMaterialInterface;
    const precTarget: Double;          // precision target of MC sampling
    const NSampleMax: Integer;        // emergency brake value for MC sampling
    const sampleSize: Integer;        // sub-sampling size for precision monitoring
    const theta0,phi0: Double;        // incidence angles
    const specToUse: TGRType;         // incidence spectrum (global radiation standard)
    const incPolState: TStokesVect;  // incidence polarisation states (usually unpol.)
    const diffTargets: Boolean;      // true: provided reflectance data is diffuse reflection, excluding gloss
                                      // false: provided reflectance data is hemispherical reflection, including gloss
    const glossAngleDeg: TStokesVect; // angular size of gloss filter (°)
    const NCP: Integer;              // max spline knot index (number of knots-1)
    const Rtargets: array of Double; // refl. targets for ranges: TOTAL, UV, VIS, IR
    const targetR, targetG, targetB: Double; // normalized target colour
    var logMemo: TMemo;              // output optimisation progress to memo
    const masterName: String): String; // filename for optimisation result
```

As can be seen in the implementation, currently, an optimisation method performing virtual measurements according to the ASTM E 903 standard is implemented. In order to align the model with the empirical data, the optimisation method offers the same functionality as the “real-world measurement” method. The incidence direction, spectrum, size of the gloss filter and type of reflectance value can be defined freely. Beyond that, the colour target, the number of free knots for the spline optimisation and the polarisation state of the incidence beam have to be provided.

### 4.13.3. Testing and validation

The feasibility of the inversion method is demonstrated based on three different surface coatings for Venetian blinds. The materials reflect the materials of the actual testing specimens used in the full system validation (chapter 8).

#### Determination of Total Solar Diffuse Reflectance (based on PI-LTM-313)

05.07.2019

Total Solar diffuse Reflectance according to ASTM E 903 weighted with direct near normal solar irradiance ASTM G173-03

<table><thead><tr><td>Sample Name</td><td>Total Solar Diffuse Reflectance</td><td>UV Range (300-395 nm)</td><td>Visible Range (400-720 nm)</td><td>Near Infrared Range (725-2500 nm)</td></tr></thead><tbody><tr><td>VSR 140 = M73256M</td><td>59.9%</td><td>52.3%</td><td>55.6%</td><td>64.2%</td></tr><tr><td>7016 = R70166</td><td>16.8%</td><td>6.7%</td><td>7.9%</td><td>25.2%</td></tr><tr><td>VSR 010 = N00516</td><td>74.2%</td><td>12.1%</td><td>80.9%</td><td>72.5%</td></tr></tbody></table>

"silver"

"anthracite"

"white"

All values lie within the reported expanded uncertainty of ±0.03 reflectance units. This uncertainty is based upon the standard uncertainty multiplied by a coverage factor k = 2, providing a coverage probability of approximately 95%.

This test conformed with all the requirements of PI-LTM-313 which is based upon ASTM E 903.

<caption>Figure 67 Reflectance data as provided by the manufacturer</caption>

Figure 67 shows the material data as provided by the manufacturer. The table provides the integrated reflectance values for the total, UV, visible and NIR ranges for three different coatings/colours. The manufacturer also specifies the testing method and spectrum used. In the product data sheet, the
material of the coating is referred to as PE (polyester). The specifications are inconsistent, as PE is usually used for polyethylene and not polyester. As most polymers exhibit refractive indices in the range of 1.4 to 1.6 in the relevant spectral range, a constant value of 1.5 is commonly assumed to model polymers generically. Alternatively, additional information on the gloss at specified angles could be used to approximate the material's refractive index. A precise measurement of the gloss value can be used to narrow down the refractive index range for the applied spectral range, as the gloss value is a measure for the specular reflection (see sections 4.3.2 and 4.3.3).

The generic value of 1.5 is assumed for the refractive index for the inversions performed here. The absorption of the material representing the binder of the coating is considered to be negligible; therefore, the complex-valued refractive index is defined as $\tilde{n} = 1.5 - 0i$.

In order to model the diffuse subsurface reflection caused by the pigments, the subsurface reflectance functions $\rho_{ss}(\lambda)$ are determined in the optimisation process described above. The results are depicted in Figure 68. As shown in Table 9, the integrated reflectance values for the different spectral ranges of optimised subsurface functions match the targeted values with reasonable accuracy. The deviations are in the range of a few tenths of a per cent. The reason for the deviation is that the optimisation based on smooth spline functions additionally has to consider the colour target. Without this objective, the optimisation would be able to provide almost perfectly matching reflectance values. Without additional color measurements available, the normalised colour target was set to red = green = blue = 1.0 for all three materials, representing a desaturated colour (white or grey).

<figure><img src="image_73.png" /><figcaption>Figure 68 Subsurface reflectance functions determined in optimisation process</figcaption></figure>

The subsurface reflectance function for the white coating shows strong fluctuations in the UV range. This is required to fulfil the steep drop from 80.9% reflectance in the visible to only 12.2% in the UV range (stated in the manufacturer's datasheet). In order to assess the quality and relevance of the optimised functions, it has to be considered that the visible and lower parts of the NIR spectral range are most relevant. The UV range contributes only 3%, and the range above 1460 nm accounts for only 10% of global radiation (see Table 2).
<table><caption>Table 9 Analysis of optimisation results for generic coated surface models</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="4">spectral ranges</th></tr><tr><th>total<br>300-2500 nm</th><th>UV<br>300-395 nm</th><th>VIS<br>400-720 nm</th><th>IR<br>725-2500 nm</th></tr></thead><tbody><tr><td rowspan="3">target<br>diffuse<br>reflectance</td><td>anthracite</td><td>16,8%</td><td>6,7%</td><td>7,9%</td><td>25,2%</td></tr><tr><td>silver</td><td>59,9%</td><td>52,3%</td><td>55,6%</td><td>64,2%</td></tr><tr><td>white</td><td>74,2%</td><td>12,1%</td><td>80,9%</td><td>72,5%</td></tr><tr><td rowspan="3">determined<br>diffuse<br>reflectance</td><td>anthracite</td><td>16,8%</td><td>7,1%</td><td>7,6%</td><td>25,5%</td></tr><tr><td>silver</td><td>60,3%</td><td>53,1%</td><td>55,7%</td><td>64,8%</td></tr><tr><td>white</td><td>74,3%</td><td>12,2%</td><td>81,0%</td><td>72,4%</td></tr><tr><td rowspan="3">deviations</td><td>anthracite</td><td>0,0%</td><td>0,4%</td><td>-0,3%</td><td>0,3%</td></tr><tr><td>silver</td><td>0,4%</td><td>0,8%</td><td>0,1%</td><td>0,6%</td></tr><tr><td>white</td><td>0,1%</td><td>0,1%</td><td>0,1%</td><td>-0,1%</td></tr></tbody></table>
# 5. Additionally applied methods

This chapter covers three topics that are essential components of the *RadiCal* method but are not directly related to light and optics. First, a short overview of the principles of Monte-Carlo methods is provided. Afterwards, the related optimisation method, simulated annealing, which is applied for several tasks in this work, is briefly introduced. Finally, the chapter concludes with a description of the geometric Fibonacci algorithm and its implementation for the *RadiCal* method.

## 5.1. Monte Carlo method

### 5.1.1. Overview

Monte-Carlo (MC) methods is a rather general term used for a variety of methods for different purposes. The methods share the approach that random numbers (samples) and the laws of probability theory are used to solve problems. The problems that are solved with MC methods in many fields, such as physics, chemistry, traffic, and economy, usually have in common that finding their analytic solution is either elaborate or impossible. They are usually hard to solve because they depend on a large number of variables that can be random or unknown to some extent. An (in)famous example is the solution of diffusion problems, as performed by Enrico Fermi and others while working on the concept of an atomic bomb at Los Alamos during world war II (Metropolis, 1987). This application represents the birth of the modern MC simulation and is also the origin of its name. "Monte-Carlo" was used as a code name that referred to the casinos of Monaco as a reference to the stochastic nature of the method (Metropolis, 1987). The MC approach is basically an application of the *law of large numbers*. This essential law of statistics states that when independent experiments subject to statistical variations are repeated a large number of times, the average of all outcomes converges to the expected value (Bernoulli, 1713; Poisson, 1837; Hsu and Robbins, 1947).

In the present work, only basic Monte Carlo methods are required, and only these are covered briefly here. Even though only basic MC methods are applied here, the present application is a good example of a problem that can, in fact, not be solved in an analytic way. Determining the transmitted power behind a shaded window would require solving a complex and large set of nested spatial and spectral integrals representing the relevant scattering processes and the incidence radiation. Providing the analytic solution to such an equation system is virtually impossible.

However, this task can be solved efficiently and straightforwardly by applying the MC method. Instead of considering the entire distributions that describe spectrally or angularly resolved flux densities, only single discrete samples matching the corresponding probabilities are used. These samples, which are random by nature but governed by distinct distributions, are fed as input parameters into a deterministic calculation. The results or output parameters are consequently subject to variation. However, by iterating this process for a large number of samples, the average of the results will converge to the desired expected value, and the measured variation can be used to assess the accuracy of the results. For the considered optical process, this means that the infinite number of continuous paths that energy flows follow are substituted by a large number of discrete random paths governed by probability functions.

Consequently, a key task of MC simulation is the generation of random samples that follow given probability distribution functions (PDF). Three simple methods relevant to this work are briefly presented below.
## 5.1.2. Simple sampling strategies

### 5.1.2.1. Sampling discrete values

Choosing discrete outcomes based on given probabilities represents the simplest sampling case. A good example is the determination of powers going into reflection, absorption and transmission at a scattering event. Since the relevant coefficients are determined by Fresnel's equations and energy conservation demands that the absorption, reflection and transmittance coefficients add up to unity ($A + R + T = 1$), the MC implementation of the scattering can be performed using merely a single uniformly distributed random number:

```
r:=random;
if r<A then AbsorptionCase else
if r<(A+R) then ReflectionCase else
else TransmissionCase;
```

Instead of following all paths, one distinct path is determined based on its probability and followed further on.

### 5.1.2.2. Sampling continuous values – inverse transform method

Provided the targeted PDF can be integrated in order to determine its cumulated distribution function (CDF) and provided this function can be inverted (CDF⁻¹), a random sample following the PDF can efficiently be generated by evaluating this CDF⁻¹ for a random number r in the range of 0 to 1, see, e.g., (Robert and Casella, 2004):

$$ 
\begin{array}{l}
r \in [0,1] \\
x = \text{CDF}^{-1}(r) \qquad \text{CDF}(x) = \int_{x_0}^{x} \text{PDF}(x) \, dx \qquad (87) \\
\rightarrow X \sim \text{PDF}(x)
\end{array}
$$

If either the analytic integration or inversion is impossible, a numerical solution can be applied instead. However, in this case, it has to be considered that only numerically efficient solutions are suitable for application in sampling methods. A numerical inversion is, e.g., applied in this work to generate wavelength samples following the global radiation spectrum (see section 4.4.2). An analytic form of the CDF⁻¹ function is, e.g., used to determine the z component of the diffuse reflection in the method localLambert (see section 4.11.2) or the ray starting point on a triangle (see section 7.10.2).

### 5.1.2.3. Sampling continuous values – rejection method

This method is applied if the CDF⁻¹ cannot be derived analytically or if solving it is computationally too expensive, but the targeted PDF f(x) can be evaluated for any x easily; see e.g. (Robert and Casella, 2004). In this case, an alternative PDF h(x) that allows drawing samples based on the inversion method is applied. This function h(x) should form an envelope on the targeted PDF function and can be scaled using a constant C so that:

$$ C \cdot h(x) \geq f(x) \qquad (88) $$

In order to achieve high performance, the function C·h(x) should cover f(x) as narrow as possible. To determine a random sample x following the PDF f(x) the following steps are performed:

1. draw sample x following h(x)

2. calculate the probability $P = \frac{f(x)}{C \cdot h(x)}$

3. draw uniform random number $r \in [0,1]$

a. if $r < P$: accept sample x

b. if $r >= P$: discard the sample
Following this algorithm, a distribution of samples X proportional to the target PDF can be generated.

More sophisticated sampling methods include, e.g. importance sampling and the metropolis algorithm (Robert and Casella, 2004). While these sampling methods were not required in the present work, a special variant of the metropolis algorithm was used for the required optimisation processes in the thesis (see section 5.2).

### 5.1.3. Sampling accuracy – the central limit theorem

Apart from generating random samples, the second essential task of MC applications is to evaluate the statistical accuracy of the results. As mentioned, the Monte Carlo method relies on stochastic input parameters fed into a deterministic process producing results that are subject to variation of initially unknown scale.

In order to assess the uncertainty of the result, the fundamental central limit theorem (CLT) can be applied. The essential theorem is presented in many forms. E.g. Kreyszig et al. (2011) state it as:

Let $X_1, ..., X_n, ...$ be independent random variables that have the same distribution function and therefore the same mean $\mu$ and the same variance $\sigma$.

Let $Y_n = X_1 + X_2 + ... + X_n$. Then the random variable

$$Z_N = \frac{Y_n - n\mu}{\sigma\sqrt{n}} \qquad (89)$$

is asymptotically normal with mean 0 and variance 1; that is, the distribution function $F_N(x)$ of $Z_N$ satisfies

$$\lim_{n \to \infty} F_N(x) = \Phi(x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{x} e^{-u^2/2}$$

Expressed in simple words, this means that all distributions formed as a sum of distributions converge to a normal distribution centred around the expected value of all distributions. Note that the initial distributions can be of any “non-pathological” type, to put it casually, but their average, always converges to a normal distribution. More specifically, random, independent samples with finite expected values and variance are required; see e.g. Montgomery and Runger (2014).

In its practical application, subsamples of size $m$ of the MC results $r_i$ are used to generate $n$ mean values $\mu_{ss,i}$ in parallel. The sample means are used to estimate the population mean $\mu_e$, that converges to the targeted expected value $\mu$ of the result (see, e.g., Montgomery et al. (2010)):

$$\mu_{ss,i} = \frac{\sum_{i=0}^{i0+m-1} r_i}{m}, \quad \mu_e = \frac{\sum_{i=1}^{n} \mu_{ss,i}}{n} \qquad (90)$$

While dividing the samples into subsamples has no impact on the result of $\mu_e$, it is relevant for monitoring the accuracy of the results. Since, unlike the $r_i$ values, the $\mu_{ss,i}$ will be normally distributed according to the CLT, the standard error of the sample means ($\sigma_{SEM}$) can be used to assess the accuracy of the estimate $\mu_e$. It can be derived as (see, e.g., Montgomery et al. (2010)):

$$\sigma_{SEM}^2 = \mathrm{var}(\mu_e) = \mathrm{var}\left(\frac{\sum_{i=1}^{n} \mu_{ss,i}}{n}\right) = \frac{1}{n^2} \mathrm{var}\left(\sum_{i=1}^{n} \mu_{ss,i}\right) = \frac{ns^2}{n^2} = \frac{s^2}{n} \qquad (91) \\ \to \sigma_{SEM} = \frac{s}{\sqrt{n}}$$
This is the simple but fundamental relation used to assess the accuracy of MC results. Since based on a normal distribution, the expression $\sigma_{SEM}$ can be used to state confidence intervals for the results.

Note that the relation holds only for the case of independent sampling, as applied in the present method. If the samples are subject to correlation, like in the case of Markov chains, more complex approaches are required to assess the accuracy, see e.g. Robert and Casella (2004). The fundamental square root dependence on $n$ in equation (91) gives us crucial information regarding the convergence speed of any MC simulation based on uncorrelated sampling. It implies that in order to achieve twice the accuracy, the number of iterations must be quadrupled; to reduce the error to

one tenth, the required iterations have to be increased a hundredfold, etc. Therefore, the convergence of MC is often perceived as fast in the beginning and slow in the end (see Figure 69).

<figure><img src="image_74.png" /><figcaption>Figure 69 Convergence of MC simulations</figcaption></figure>
## 5.2. Simulated annealing

Simulated annealing is a probabilistic optimisation method proposed in 1983 by Kirkpatrick et al. (1983). However, the main idea behind the algorithm dates back significantly longer to a key paper by Nicholas Metropolis in the year 1953 (Metropolis et al., 1953). The combinatorial optimisation method allows finding the minima (or maxima) of functions that depend on a large number of variables. Thermodynamical observations have inspired the method, and an actual annealing process can serve as an illustrative example to explain its principles. It is commonly known that perfectly symmetrical single crystals can be grown from a melt or a saturated solution. While we are used to this fact, it is

actually quite astounding that trillions of molecules or atoms can arrange themselves in a way that their macroscopic arrangement reflects an energetically ideal shape. How is this possible? The answer to this question is the principle behind the simulated annealing optimisation algorithm: Firstly, excess energy determined a temperature measure is provided to the system to allow an efficient reconfiguration of its constituents. These reconfigurations also include transitions to energetically less optimal configurations. Second, the temperature, i.e. the associated excess energy, is slowly reduced, and the probability for less ideal configurations decreases along with it. In order to determine what configurations are possible in such a process, first some fundamentals of thermodynamics, or statistical mechanics to be more precise, have to be considered. At thermal equilibrium the probability that a system with N potential states can be found in the state i having energy $E_i$, is described by the Boltzmann distribution:

<figure><img src="image_75.png" /><figcaption>Figure 70 simulated annealing progression of cost function blue: cost of current solution every 10.000th iteration red: cost of optimal solution found so far</figcaption></figure>

$$p_i = \frac{e^{-\frac{E_i}{kT}}}{\sum_{j=1}^{N} e^{-\frac{E_j}{kT}}} \qquad (92)$$

Metropolis proposes an efficient way to populate these states numerically in his publication (Metropolis et al., 1953). In order to generate configurations matching the Boltzmann distribution, the following strategy is applied:

0. start with state j and calculate its probability

1. create a new state i and calculate its probability

2. if $p_i \ge p_j$ ($E_i < E_j$): accept the new state always

3. if $p_i < p_j$: accept the new state with the probability given by $p_i/p_j$

4. proceed to step 1
If the temperature is gradually reduced during this process, the system moves effectively to its ideal state, as the temperature dependence of $p_i/p_j$ increasingly rules out less ideal states (see Figure 70). However, since at higher temperature levels, less optimal configurations are still accepted, it can be avoided that the optimiser gets trapped in local minima early. If the annealing process were given infinite time, it would reach the optimal solution, i.e. the absolute minimum of the cost function, as it would escape from any local minima again. In its practical application, with limited processing time, it cannot be guaranteed that the ideal state can be identified. However, simulated annealing is able to identify good solutions for such problems very efficiently.

In its practical implementation, the following four components need to be implemented to perform a simulated annealing optimisation (see, e.g., Ziegel, *et al.* (2007), Penna (2002)):

**1. Definition of a state**

A concise and unambiguous structure containing all state parameters of the system has to be provided.

**2. Cost function**

A function that is used to evaluate the optimality of a state (its “energy”). The target of the optimiser is to find the absolute minimum of this function. Often the function is defined in negative values. However, only the difference between two evaluations of the cost function is considered in the algorithm. Therefore, its absolute value is irrelevant.

**3. Alteration algorithm**

An algorithm that is able to generate random new states is necessary. This function has a significant impact on the efficiency of the optimisation. While it is generally possible to create new states independent of the current state, it proved that an algorithm that is able to generate incremental adaptions of the current states is more efficient for most problems. Generally, a combination of both is applied.

**4. Temperature and cooling schedule**

Following the thermodynamic fundamentals, the cooling schedule usually provides a logarithmic temperature decay. However, the actual cooling schedule depends on the cost function and the nature of the optimisation task. In most cases, “phase changes” in the system can be observed. These represent temperature regions where the cost function remains constant to drop sharply afterwards. Since this implies that the system transitions into a path of new configurations, making it harder to access previous states, it is advisable to slow down the cooling at such critical temperature ranges. This can, e.g., effectively be achieved by tracking the average costs with two exponential moving averages of different period lengths. If the deviation of the two measures indicates a sharper decrease, the cooling speed should be reduced, whereas it can be increased if the ratio of the averages indicates little change.

In order to demonstrate how simple the core of the annealing algorithm is, a simple example is provided below. All applications of the algorithm used in this work follow the presented scheme. Only the usually applied additional features to adjust the cooling schedule dynamically are omitted below for reasons of clarity.
```pascal
TDec:=power(0.9,DurationScale/TargetDuration); // decrement for logarithmic cooling
T:=InitialTemp; // set starting temperature
setInitialState; // set initial state
it:=0;

Repeat
    Inc(it); // increment iteration counter
    saveOldState; // backup/save the current state
    Alter; // generate new state
    cost:=getCost(false); // eval. cost-function for new state
    delta:=cost-oldCost; // delta of new vs. Old

    if (delta<=0) or (random<exp(-delta/T)) then // accept if better, or
        if worse with Boltzmann probability
        Begin // additionally capture the current state if best so far
            if (cost<bestCost)and(it>100000) then saveBest;
        end
    else
        restoreOldState; // discard new state and go back to the old state

T:=T*TDec; // lower temperature
until (someAbortCriterion);
```

As can be seen, the algorithm described above can be simplified to calculating the difference of both cost functions and using it to determine the probability for a new, energetically less optimal state by calculating the probability $exp(-delta/T)$ as this is equivalent to the ratio $p_i/p_j$.
## 5.3. A Fibonacci hemisphere

### 5.3.1. Background and theory

In the method presented here, the scattering or absorption behaviour of the analysed system is not provided in a discrete form (like the Klems matrix) but in a functional form. However, in order to derive this functional form (*SIOP*, see section 6.1) it is necessary to provide discrete measurement samples of the target quantity for different incidence directions on the hemisphere. The most suitable distribution of incidence directions for the intended technical application should cover the hemisphere with equal spacing. A regular angular pattern, as provided by the Klems segmentation, does not seem ideal for that purpose. The reason for that is threefold: first, the two-axis symmetry (vertical/horizontal) of the angular Klems pattern will often lead to an inefficient sampling of redundant information, as the analysed object (e.g. shaded window), will often also show symmetry in the vertical and/or horizontal axis. In the case of a rotationally symmetric behaviour (e.g. of unshaded glazings or textile screens), the full Klems matrix can be derived from only nine incidence directions. Second, the regular spacing of the two angles is more prone to sampling issues, e.g. aliasing, since the analysed applications, like a window shaded by tilted Venetian blinds, will also show a pronounced regular angular pattern. Sampling this angularly equally spaced pattern at a few equally-spaced detection angles can easily lead to aliasing effects if these patterns are correlated. Finally, the angular resolution of the Klems matrix decreases for higher incidence angles $\theta$ and is therefore less able to capture gradients of shadings or reflectivity increases that are dominant at these angles (see below).

<figure><img src="image_76.png" /><figcaption>**Figure 71** Fibonacci pattern on the succulent plant *Epithelantha*

Photo by: Morgoth92/Wikimedia (CC 2.0)</figcaption></figure>

Following these arguments, it seems more appropriate to distribute the sampling directions evenly on the hemisphere with little symmetry regarding its angular distribution. Providing an even distribution of directions (points) on a sphere is a theoretical topic with a long history and high relevance in various scientific disciplines.

Various approaches are available for that purpose, like the icosahedron algorithm and numerous methods based on numerical optimisation. Considering all requirements stated above, the Fibonacci sphere method seems to provide the ideal distribution of directions here:

a) it provides an almost perfect evenly spaced distribution on the sphere (hemisphere)

b) it shows little symmetry regarding the vertical or horizontal axis

c) it can be generated for any desired number of sampling directions N

d) it is a direct and fast method (as it is not based on optimisation)

The angular or spatial patterns that are based on Fibonacci distributions are frequently found in nature, for example, the arrangement of seeds (e.g. sunflowers) or leaves (e.g. aloe plants). The frequent occurrence in nature is generally a good indication regarding the efficiency and optimality of the distribution.
### 5.3.1.1. The Fibonacci sphere algorithm

In order to derive a Fibonacci sphere, firstly, a Fibonacci lattice on the unit square is considered. The coordinates for the corresponding vertices $p_i$ can be written as (Marques et al., 2013):

$$p_i = \begin{cases} \left(i \frac{F_m}{F_{m-1}} \% 1, \frac{i}{F_m}\right) & (93) \\ \text{for } i = 1.. F_m & \end{cases}$$

Where $F_m$ and $F_{m-1}$ are the elements of the Fibonacci number series, and % stands for the modulo operator.

Using a Lambert cylindrical equal-area projection, the pattern can be mapped on a sphere. Additionally, the ratio $\frac{F_m}{F_{m-1}}$ is substituted with the golden ratio $\Phi = 1/2(1+\sqrt{5})$, as the ratio of consecutive Fibonacci numbers quickly converges to this ratio for higher m values (Graham et al., 1994). Consequently, the vertices for a Fibonacci sphere with $N = F_m$ elements can be expressed in polar coordinates as:

$$p_i = (\varphi_i, \theta_i) = \left(i \frac{2\pi}{\Phi}, \arcos\left(1 - \frac{2i}{N}\right)\right) \quad (94)$$

<figure><img src="image_77.png" /><figcaption>Figure 72 Fibonacci sphere generation
$\varepsilon$ offset applied for optimised distribution</figcaption></figure>

### 5.3.1.2. A (slightly) optimised Fibonacci hemisphere

This canonical Fibonacci distribution can be optimised by inserting an offset $\varepsilon$ at the beginning and the end of its vertical distribution (94). Further, the distribution is shifted in a way that it is symmetrical in its z coordinates, as the original distribution has its top vertex slightly offset ($z = 2/N$) and its bottom vertex at the bottom pole ($z = -1$). Hence, the z coordinates of the N points, that are initially equally spaced within the interval $[-1, 1 - 2/N]$, are now mapped onto the interval $[-1 + \varepsilon, 1 - \varepsilon]$. The distribution of the directions in polar coordinates can then be written as:

$$p_i(\varepsilon) = (\varphi_i, \theta_i) = \left(i \frac{2\pi}{\Phi}, \arcos\left(1 - \frac{2(i + \varepsilon)}{N - 1 + 2\varepsilon}\right)\right) \quad \text{for } i = 0..N-1 \quad (95)$$

In Figure 72, it can be seen how the first vertex is offset from the pole along the helix formed by the vertices. The vertical distance of this offset is defined by the parameter $\varepsilon$. The parameter $\varepsilon$ can now be used to “fine-tune” the distribution of the directions. Deviating from the proposal made by Roberts (2020), a constant value $\varepsilon_0$ for all N values is sought here. This shall, on the one hand allow simple replicability of the distribution and on the other hand the value of $\varepsilon$ converges quickly for - the here relevant - higher N values anyway. The maximization of the average minimum distance $\bar{\alpha}$ of all vertices has been chosen as the criterion for the optimisation process:

$$\bar{\alpha} = \frac{1}{N} \sum_{i=0}^{N-1} \alpha_i \rightarrow \max \quad (96)$$
with $\alpha_i$ being the minimum angle between the vertex $i$ and its closest neighbour. Additionally, since - for the application intended here - only the hemisphere is considered, a new constant $N_{\text{hemi}}$ is introduced. $N_{\text{hemi}}$ represents the number of Fibonacci vertices that are provided for the (upper) hemisphere. However, $N_{\text{hemi}}$ does not need to be identical with the half of the $N$. Rather, for any desired $N_{\text{hemi}}$, $N$ is determined in a way that precisely $N_{\text{hemi}}$ vertices of the upper hemisphere are located above a certain threshold $\xi$. This threshold value is introduced to avoid grazing incidence angles, i.e. almost parallel directions. These grazing directions can be detrimental, as they can increase the calculation time or lead to numerical issues. The height threshold value, which applies to the cosine of $\theta$ (=z-coordinate), has been defined as $\xi = 0.001$.

Figure 73 Results of the optimisation process to determine

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

In the optimisation process, the optimal values for $\varepsilon$ according to equation (96) for hemispheres with an increasing number of directions $N_{\text{hemi}}$ (up to 10,000) were analysed. The results of the optimisation process are depicted in Figure 73. It can be seen that while there are different solutions for smaller $N_{\text{hemi}}$ distributions, the value of $\varepsilon$ converges for higher numbers. Therefore, a constant value of $\varepsilon_0 = 0.3613$ is defined and universally applied in the following.

Figure 74 and Figure 75 show the angular segmentation of the hemisphere according to the classical Klems approach and the implemented Fibonacci approach as polar plots and 3D representations. As can be seen, the Klems matrix shows high symmetry in multiple axes as well as angular correlations. Both features are not ideal for the *RadiCal* approach. By contrast, the Fibonacci segmentation shows no axial symmetries and little correlation. It can also be seen that due to the actual implementation, the first Fibonacci patch is slightly offset from the pole, which is intended.

<figure><img src="image_79.png" /><figcaption>Figure 74 2D polar representation of angular segmentations of the hemisphere
Klems vs. RadiCal implementation of Fibonacci algorithm
(left: Klems N=145, centre: Fibonacci N=145, right: Fibonacci N=1296)</figcaption></figure>
<figure><img src="image_80.png" /><figcaption>Figure 75 3D view of angular segmentations of the hemisphere
(left: Klems N=145, centre: Fibonacci hemisphere N=145, right: Fibonacci hemisphere N=1296)</figcaption></figure>

### 5.3.1.3. Analysis of distributions

In order to analyse the quality of the angular distributions, the corresponding solid angles are analyzed. For the Klems segmentation, these are explicitly provided in its definition, while for the Fibonacci hemisphere, they are computed based on a proximity criterion. The solid angles of the Klems distribution show a significant variation due to the regular segmentation in $\theta$. The solid angles of the angular patches increase with the incidence angle $\theta$. This was obviously a design choice of Klems, as the relevance of these directions will, for many applica-

<figure><img src="image_81.png" /><figcaption>Figure 76 Comparison of relative solid-angle sizes for three different distributions of directions on the</figcaption></figure>

tions, usually decrease proportionally with $cos(\theta)$. On the other hand, it can be argued that higher gradients of the represented functions (such as the transmission of the window) can be expected mainly for these higher incidence angles. For example, the area shaded by protruding parts or the reveal will increase proportionally to $tan(\theta)$. Hence, correctly capturing this characteristic requires higher sampling resolution for higher incidence angles. Considering both arguments, sampling with equally spaced directions, as provided by the Fibonacci hemisphere, seems to be the most appropriate approach. Figure 76 shows a comparison of the solid angles that correspond to the angular patches provided by three different distributions. It can be seen that the patches of the Klems distribution for the highest incidence angles (>75°) are approximately six times larger than the ones near normal incidence (0°). In contrast, the Fibonacci distributions provide almost perfect equally-spaced directions for most incidence angles. Only in regions near grazing incidence (close to 90°) a fluctuation of the patch sizes is apparent. The reason for this is the clipping effect by the surface plane and the threshold height criterion.
### 5.3.2. Implementation

In the implementation, all algorithms relevant for setting up and evaluating the Fibonacci hemisphere are contained in the class `TFibonacciHemiSphere` of the library `rc_fibonacci_sphere`.

```typescript
TFibonacciHemiSphere = class (TDataModule) { TFibonacciSphere }
```

**private**

```typescript
eps: Double;
nHemi: Integer;
nFull: Integer;
fbDir: array of Tvect3d;
fbTheta: array of Double;
fbPhi: array of Double;
fbSolidAngles: array of Double;
goldenRatio: Double;
```

**public**

```typescript
// main functions
procedure getFiboHemisphereAngles(const i, nHemiSphere: Integer; var ph,th: Double);
function getFiboHemisphereDir(const i, nHemiSphere: Integer): Tvect3d;
function getSolidAngle(const i, nHemiSphere: Integer): Double;
procedure setupHemi(const nHemiSphere: Integer);
function nMax: Integer;

// aux. functions
constructor Create;
function meanMinAngle: Double;
procedure calcSolidAngles(const res: Integer);
function ReadSAFromCache(const nHemiSphere: Integer): boolean;
procedure WriteSAToCache;
function getClosestI(const v: Tvect3d; const nHemiSphere: Integer): Integer;
function dumpFiboHemi: String;
function findOptEps: String;
deleterctor Destroy; override;
end;

const minHeightThreshold = 1E-3;
const eps0 = 0.3613;
const goldenRatio = 1.6180339887499; // (1 + sqrt(5))/2
```

The angular and vector-based data is precomputed to increase performance. The precalculation can either be performed directly by calling the procedure `setupHemi` or indirectly by calling `getFiboHemisphereAngles` or `getFiboHemisphereAngles` for the first time, as they will, in turn, call the setup procedure.

The actual and remarkably simple Fibonacci algorithm is contained in the setup procedure. It is presented in full here:

```typescript
procedure TFibonacciSphere.setupHemi(const nHemiSphere: Integer);
var cosTheta: Double;
    i: Integer;
    per: double;
begin
    nHemi:=nHemiSphere; // map to global var

    // determine correct N for Fibo sphere to get nHemisphere directions above the height threshold
    nFull:=2*nHemi-4;
    Repeat
        Inc(nFull);
        // precompute cosTheta of last entry on hemisphere
        cosTheta:=1 - 2*(nHemi-1)+eps)/(nFull-1+2*eps);
        until cosTheta<minHeightThreshold; // make sure last direction (i-1) is slightly greater than parallel

        SetLength fbDir,nHemi);
        SetLength fbTheta,nHemi);
        SetLength fbPhi,nHemi);

        per:=0;
        for i:=0 to nHemi-1 do // note: (nHemi-1)=(i max)
            Begin
```
```matlab
fbphi[i]:=2*pi*i/goldenRatio;
if fbphi[i]-per*2*pi>=2*pi then per:=per+1; // map phi to 0..2*pi
fbphi[i]:=fbphi[i]-per*2*pi;

fbtheta[i]:=arccos(1 - 2*(i+eps)/(nFull-1+2*eps));
fbDir[i]:=vect(cos fbphi[i]) * sin fbtheta[i]), sin fbphi[i])
    * sin fbtheta[i]), cos fbtheta[i]]);

end;
end;
```

In the first part of the setup procedure, the correct value for the number of directions on the entire Fibonacci sphere *N* is determined. *N* is determined in a way that ensures that precisely the desired *N*<sub>hemi</sub> are located above the height threshold. Afterwards, the angles and vectors that define the *N*<sub>hemi</sub> directions are precomputed based on equation (95) and stored in the corresponding arrays.

### 5.3.3. Testing and validation

During the implementation of the algorithm and, in particular, during the optimisation of the parameter $\epsilon_0$ several tests and validations of the algorithm were performed. The distribution of solid angles and angles among the directions were computed and analysed for many different $N_{hemi}$ values. Beyond that, the visual representations, as depicted in Figure 74 and Figure 75, are a reliable indicator for the validity of the implemented code.
# 6. Solar Incidence Operators
## Definition, generation and evaluation

### 6.1. SIOP – Solar Incidence Operator

#### 6.1.1. SIOP concept and motivation

In order to describe the scattering characteristic of a surface (or an object comprising multiple surfaces), usually a scalar function $f(\theta, \varphi)$ of the zenith and azimuthal angle defined on the unit sphere's surface is used. The scalar function then represents either the transmittance coefficients or reflectance coefficients in an angularly resolved way. However, to capture the scattering characteristics comprehensively, it is necessary to consider the direction of the incoming, as well as of the outgoing light direction. For this purpose, the concept of the bi-directional scattering-distribution-function BSDF, which is then a four-dimensional function $f(\theta_i, \varphi_i, \theta_o, \varphi_o)$, depending on the incoming as well as the outgoing direction, is applied. In practical application, the function is either defined by an idealized analytic expression or generated based on empirical data. For the latter, the measured data is provided as discrete samples, e.g. based on the Klems-matrix containing 145x145 elements (see section 3.2). In this approach, incoming and outgoing directions are only considered on the hemispheres, as this is sufficient for most (but not all) applications.

The SIOP approach introduced here is primarily intended for energetic modelling and not for optical modelling. Therefore, the SIOP is a two-dimensional scalar function considering only the incident light's angular distribution. Consistent with the common approaches, only the incident radiation distributed on the hemisphere ($\theta = 0...\pi/2$) is considered, as this will be sufficient for most practical applications.

The main objective of the SIOP definition is to overcome issues resulting from the discrete nature of sampled data. While transforming the data in the four-dimensional space of the Klems matrix into a smooth function is a hardly solvable task, the restriction of the SIOP definition to the two dimensions of the incidence light makes this much easier. Representing the scattering properties in a smooth, differentiable functional form provides significant benefits:

* integration: the SIOP can easily and precisely be integrated over any region resulting again in a smooth function (see section 6.2).

* differentiation: the smooth derivative of the function means that it can provide smooth time series at any resolution as optimal input for dynamic simulation

* interpolation: using suitable basis functions, optimal interpolation is possible

* fitting: the use of suitable basis functions allows optimal reconstruction of underlying functions that are provided in the form of “noisy” measurement samples

* data-compression: suitable basis functions allow a high level of data-compression

In the implementation, spherical harmonics have proven to represent a suitable set of basis functions for this purpose, as will be shown in the following sections. The developed method shows analogies to MP3 coding for audio data (Musmann, 2006) or the JPEG coding algorithm for images (Usevitch, 2001), as original sampled data is represented using a suitable set of basis functions.
The SIOP function can be used to determine any physical quantity q that represents a scalar value and depends on an incidence profile I(θ, φ):

$$q = \iint_{hemi} (SIOP_q \cdot I) d\Omega = \int_{\theta=0}^{\pi/2} \int_{\phi=0}^{2\pi} SIOP_q(\theta, \phi) \cdot I(\theta, \phi) \sin(\theta) d\phi \, d\theta \quad (97)$$

While the obvious application for the SIOP is the determination of transmitted, reflected and absorbed power using incidence profiles containing radiance information, SIOPs can also be used for other purposes. Exemplary applications of radiometry and photometry are provided in Table 10. Further applications can, e.g., involve evaluation of glare parameters, degree of polarisation, angle of polarisation, photosynthetically active radiation (PAR), colour temperature etc.

<table><caption>Table 10 Exemplary list of potential SIOP target quantities and relevant units</caption><thead><tr><td>target quantity</td><td>Unit of SIOP</td><td>Incidence profile quantity</td><td>Units of incidence profile</td></tr></thead><tbody><tr><td>Transmitted, reflected or absorbed power [W]</td><td>W/(Wm<sup>-2</sup>)<br>= m<sup>2</sup></td><td>radiance</td><td>Wm<sup>-2</sup>sr<sup>-1</sup></td></tr><tr><td>Transmitted, reflected or absorbed energy [J]</td><td>J/(Jm<sup>-2</sup>)<br>= m<sup>2</sup></td><td>radiant energy</td><td>Jm<sup>-2</sup>sr<sup>-1</sup></td></tr><tr><td>Transmitted, reflected or absorbed flux density [Wm<sup>-2</sup>]</td><td>(Wm<sup>-2</sup>)/(Wm<sup>-2</sup>)<br>= 1</td><td>radiance</td><td>Wm<sup>-2</sup>sr<sup>-1</sup></td></tr><tr><td>Spectral flux density [Wm<sup>-3</sup>]</td><td>(Wm<sup>-3</sup>)/(Wm<sup>-3</sup>)<br>=1</td><td>spectral radiance</td><td>Wm<sup>-3</sup>sr<sup>-1</sup></td></tr><tr><td>Illuminance [lm m<sup>-2</sup>]</td><td>(lm m<sup>-2</sup>)/(lm m<sup>-2</sup>)<br>= 1</td><td>luminance</td><td>cd m<sup>-2</sup><br>= lm m<sup>-2</sup>sr<sup>-1</sup></td></tr><tr><td>Luminance flux [lm]</td><td>(lm)/(lm m<sup>-2</sup>)<br>= m<sup>2</sup></td><td>luminance</td><td>cd m<sup>-2</sup><br>= lm m<sup>-2</sup>sr<sup>-1</sup></td></tr></tbody></table>

Some remarks regarding the application of SIOPs:

**Time-integrated incidence profiles:** Due to the linearity of light (section 4.2.2), it is always possible to compute time-integrated target quantities when the SIOP is evaluated for a time-integrated incidence profile. This means that, e.g., total energy values for any time interval (hour, day, month etc.) can be computed within a single SIOP evaluation step if time-integrated incidence profiles for that interval are applied. The resulting very short evaluation times in the order of Milliseconds allow for applying optimisation algorithms that rely on a large number of iterations.

**SIOP configurations:** Apart from the different dimensionality, another important difference between the SIOP and Klems/BSDF calculation method lies in the fact that for the evaluation of any SIOP, always the initial, “undisturbed” incidence profile relevant to the entire system is used. In contrast, for evaluating a system comprising multiple BSDF layers, the undisturbed incidence profile is applied only on the first layer, whereas the further layers will be evaluated using the output of the preceding layer. Therefore, the BSDF/Klems approach allows using single-layer data individually or rearranging it into new configurations. However, the approach is linked to potential sources of error in the practical application (see section 3.3). In contrast, SIOPs of any layer should only be used in the context of the system they were calculated for.

**Level of detail target object:** the definition of the SIOP naturally does not contain any restrictions regarding the target objects. Since the SIOPs are not based on individual layers but are generated and evaluated within the context of layer systems, they are intended to contain detailed information
about the system. For example, the SIOPs for a window detail can and should be modelled based on a model containing details regarding its frame and/or reveal instead of only considering the glazing layers. Of course, this means that a large number of configurations has to be covered. However, by scaling the SIOP results based on a reference area, which should be contained in the description of the SIOP, it is possible to interpolate and extrapolate results using the “nearest” suitable model.

**SIOP database:** An intended open online database for SIOPs would allow users to share SIOPs for a wide variety of models and applications. If the models are provided with detailed information regarding their generation (model, reference area, material), they can straightforwardly be applied by any user and for any type of application.

**SIOP for lighting application:** on the first view, it seems that the limitation of SIOPs to a scalar value on the output side prevents the direct application of the concept in the lighting field, as directional information is usually required there. However, there are ways to overcome this limitation. The target quantity of the SIOP can, of course, also be subject to specific directional restrictions. For example, 145 SIOPs could be generated to represent the light transmitted in all directions of the KLEMS matrix. Alternatively, any other directional separation could be applied. Generally, any directional segmentation suitable for the used lighting software could be considered, e.g. separation into fewer directions, separation into diffuse and direct components only, etc. Likewise, the separation of wavelengths into any desired primary colour basis can be straightforwardly applied. However, lighting applications are not part of the present work.

## 6.1.2. SIOP definition and generation

### 6.1.2.1. Definition based on spherical harmonics

The SIOP representation relies on spherical harmonic expansion, see e.g. (Abramowitz et al., 1988; Ziegel et al., 2007). Spherical harmonics are well known for their use in quantum physics ("orbits"). As their German name “Kugelflächenfunktionen” (sphere-surface-functions) expresses, they represent a set of orthogonal functions defined on a sphere’s surface and are therefore well suited to serve as basis functions for the expansion of a general function defined on a sphere. The method is similar to the Fourier expansion of functions in 1D. Since spherical harmonics are intrinsically based on circular functions, they have proven to be very well suited for the application described here, as scattering processes are largely governed by trigonometric functions. In a general, complex-valued form, the spherical functions can be written as (Ziegel et al., 2007):

$$Y_{lm}(\theta, \varphi) = C_l^m \cdot P_l^m(\cos \theta) \cdot e^{im\varphi} \qquad (98)$$

$$C_l^m = \sqrt{\frac{2l + 1 (l - m)!}{4\pi (l + m)!}} \qquad (99)$$

$$l \ge 0 \quad \text{and} \quad -l \le m \le +l \qquad (100)$$

Where $P_l^m$ are the associated Legendre polynomials (Abramowitz et al., 1988; Kreyszig et al., 2011). Hence, for every positive integer $l$ there is a subset of $2l + 1$ independent functions that are iterated through by the sub-index $m$. Starting off with a constant function (i.e. a sphere) for $l=0$, the functions become increasingly complex and show an increasing number of local extrema with higher $l$ values. As mentioned, the set of spherical harmonics possesses the remarkable (almost astounding) property of being orthonormal, i.e. the product of any two different functions will vanish, and the integral over squared function on the entire hemisphere equals one for any spherical harmonics:
$$ \iint_{sphere} Y_{l_1 m_1} \cdot Y_{l_2 m_2} d\Omega = \begin{cases} 1 & \text{for } l_1 = l_2 \land m_1 = m_2 \\ 0 & \text{otherwise} \end{cases} \qquad (\text{orthogonality}) \qquad (101) $$

$$ \iint_{sphere} Y_{lm}^2 d\Omega = 1 \qquad (\text{normalisation}) \qquad (102) $$

For the definition of the real-valued SIOP function, it is sufficient to evaluate only the real part of the spherical harmonics, subsequently denoted here as YR. Utilizing a complex conjugate identity for negative m values, the function can be written as (Sommers, 2001):

$$ \begin{aligned} m = 0: & \quad YR(\theta, \varphi) = C_0^0 P_0^0(\cos \theta) \\ m > 0: & \quad YR(\theta, \varphi) = \sqrt{2} C_l^m P_l^m(\cos \theta) \cdot \cos(m\varphi) \\ m < 0: & \quad YR(\theta, \varphi) = \sqrt{2} C_l^{|m|} P_l^{|m|}(\cos \theta) \cdot \sin(m\varphi) \end{aligned} \qquad (103) $$

This formulation allows a straightforward software implementation. An efficient recursion formula (Ziegel et al., 2007) is used to evaluate the Legendre polynomials $P_l^m$ numerically. The scaling factors $C_l^m$ containing factorials (eq. (100)) are pre-computed to ensure fast evaluation. The shape of the first 12 functions is depicted in Figure 77.

For applying the functions as the basis for series expansions, it is helpful to enumerate them using a single index j. For this purpose, the two indices l and m are mapped onto a single series, using l as the first-order and m as the second-order enumerator (Bielewicz et al. 2004):

$$ j = l^2 + l + m \qquad (104) $$

This definition allows to represent the SIOP by using a single, one-dimensional vector $\vec{s}$, where the coordinates of the vector $s_j$ represent the coefficients of the series-expansion. The SIOP can then be evaluated as a weighted sum of the real-valued spherical harmonics functions:

$$ \begin{aligned} \text{SIOP}(\theta, \varphi) &= \vec{s}_j \cdot \overrightarrow{YR_j}(\theta, \varphi) \\ &= \sum_{j=0}^{jmax} s_j YR_j(\theta, \varphi) \end{aligned} \qquad (105) $$

<figure><img src="image_82.png" /><figcaption>Figure 77 polar plot of the first 12 real-valued spherical harmonic functions (j=0..11)
red: positive values; blue: negative values;
(only the relevant hemisphere is depicted;
rendering created by author)</figcaption></figure>
### 6.1.2.2. Sampling data for SIOP generation

To generate a specific SIOP for an application, such as modelling the solar transmittance of a shaded window, angularly resolved information in the form of a discrete set of sample values ($\vec{d}_i, v_i$) is used. The scalar values $v_i$ represent the desired quantities (e.g. the transmitted power) for light restricted to a single incidence angle, i.e. a parallel beam from the discrete direction $\vec{d}_i = (\theta, \varphi)_i$. Although the algorithm can be performed for any set of directions, a distribution of the sampling directions with a non-regular or non-uniform pattern in an angular sense is recommended. Unlike the angular distribution of the Klems matrix (see section 3.2), a distribution with less symmetry regarding the angles is less prone to aliasing effects. Considering these requirements, the Fibonacci sphere geometry seems to provide the ideal distribution of directions for this application, as discussed in section 5.3.

The target output values $v_i$ for the sampling directions will usually be generated based on virtual measurements, applying appropriate raytracing software as proposed in this work. However, provided a suitable measurement system is available, the input data could as well be of empirical nature.

### 6.1.2.3. Determination of SIOP coefficients

Even though the unique property of spherical harmonics of forming an orthonormal basis would allow direct calculation of the reconstruction coefficients by integration (see, e.g. Sommers (2001) or Liu et al. (2013)), it has proven that indirect determination of the coefficients employing parameter optimisation leads to better results for the application intended here. The reason for preferring optimisation is that the measured samples $v_i$ are subject to statistical errors. Even if these errors are minor, the direct method will lead to a precise reconstruction of all samples. This means that the locally contained noise will also be reconstructed. Therefore, direct reconstruction by integration tends to generate oscillating artefacts in between the given sampling points.

Finding the coefficients based on parameter optimisation instead of using direct integration technically turns the interpolation method into a fitting method, or - ideally - into a quasi-interpolation of the unknown underlying function without any statistical errors.

In order to find this smooth, continuous representation of the measurement using real-valued spherical harmonics, the non-linear parameter optimisation simulated annealing is used (see 5.2). The stochastic method is simple to implement but has proven to be very powerful in performing global optimisations on functions that depend on many parameters and have numerous local minima (non-convex functions).

For this application, the following cost function $f_{opt}$ is minimised during the optimisation process:

$$f_{opt}(s_j) = \left( \frac{1}{N} \sum_{i=1..N} [v_i - SIOP(\theta_i, \varphi_i)]^2 \right) \left( 1 + \frac{6k}{(j_{max} + 1)} \right) \rightarrow \min \quad (106)$$

$s_j$ ... spherical harmonics weighting coefficients

N ... number of measurement points

$v_i$ ... measured value at direction $d_i$

k ... number of non-zero coefficients

$j_{max}$ ... highest order allowed (predefined)

The first factor of the cost function describes the squared deviation of the SIOP reconstruction vs. the measurement data $v_i$ and represents the main optimisation target. The second factor represents a penalty in the optimisation process that increases with the number of non-zero expansion coefficients, denoted by $k$. The term is not essential to the optimisation process but has proven to be useful in reducing the amount of very small, and therefore negligible, coefficients during the optimisation. This ensures a more efficient (=faster) evaluation of the SIOP later on and tends to produce smoother representations if only little sampling information is available.
Finally, the resulting SIOP can be written and deployed in the form of a vector containing the real-valued coefficients $s_j$ of the expansion:

$$SIOP(\theta, \varphi) = \vec{s} \cdot \overrightarrow{YR}(\theta, \varphi) \qquad (107)$$

Since the vector usually contains mostly zero elements, it is practically more efficient to deploy (export) it in the form of an indexed array containing only non-zero coefficients $c_i$ and their indices $j_i$ (e.g. in XML or binary format):

$$\{SIOP\} = \{ (j_1, c_1), (j_2, c_2), \dots (j_n, c_n) \} \qquad (108)$$

### 6.1.2.4. SIOP reference system and polar plot representation

The data contained in the SIOP can best be depicted in the form of three-dimensional polar plots. For this purpose, the actual value of the SIOP determines the radius (measured from the origin) for the two spherical angles $(\theta, \varphi)$ for any point on the hemisphere. By convention, the right-handed reference coordinate system for the SIOP is defined in a way that the z-axis points towards the direction of normal incidence and the y-axis is aligned to the vertical axis of the object (or any other reference axis, if not available). The xy-plane, therefore, reflects the reference plane of the object (e.g. the plane parallel to glazing). For consistency, any SIOPs should be generated using this orientation.

The z-axis, which points towards the direction of normal incidence to the reference plane, also represents the main symmetry axis for the evaluation of the spherical harmonics (corresponding to $\theta = 0^\circ$), whereas the azimuthal angle $\varphi$ is defined on the xy-plane. The polar plot representation for the SIOP and the reference axis and angles are depicted in Figure 78 (the exemplary SIOP represents the transmission through a shaded window with 30° tilted blinds).

<figure><img src="image_83.png" /><figcaption>Figure 78 SIOP polar plot representation and reference system</figcaption></figure>

xy: reference plane (e.g. window plane),
$\theta, \varphi$: angles for spherical harmonics evaluation
### 6.1.3. Implementation

All data representing a *SIOP* is contained in the class `TRCSIOP` of the module `rc siop`. The class `TRCSIOP` contains an instance of the class `TFibonacciSphere` that is initialised with the number of directions (`ssNDir`) representing the sampled angular distribution. This means that the sampling directions are not stored in the file but are regenerated using `TFibonacciSphere`. The class `TRCSIOP` also includes an instance of the class `TSphericalHarmonics` that provides the algorithms necessary to evaluate the spherical harmonics according to equations (103). Further, it temporarily contains the azimuthal and tilt angle of the *SIOP* used for the most recent evaluation process, as well as the associated vector transformation matrix.

```python
type TRCSIOP = class(Orthogonalizer)(TRCSIOP)

private
    spHarm: TSphericalHarmonics;
    fiboHemi: TFibonacciSphere;
    tilt,azi: Double;
    rotMat: array[0..2] of Tvect3d;
    procedure precalcHemiMeansStdSYS;
    function performanceTest: String;
public
    ssiop: TSiopData;
    constructor Create;
    procedure loadFromFile(const fn: String);
    procedure save(const fn: String);
    destructor Destroy; override;
    // export functions
    procedure saveSamplesAsPNG(const fn: string; const dim: Integer; const minValue: Double;
        const maxValue: Double; const palId: Integer; const numcols: Integer);
    procedure saveFuncAsPNG(const fn: string; const dim: Integer; const minValue: Double;
        const maxValue: Double; const palId: Integer; const numcols: Integer;
        const logScale: Boolean);
    procedure saveSamplesAsLobe(const fn: string; const objName: String; const dim: Integer;
        const maxValue, wallDim, scalingFct: Double; const palId,
        palcols: Integer; const asPolar, withWall: Boolean);
    procedure saveSIOPAsLobe(const fn: string; const objName: String; const dim: Integer;
        const maxValue, wallDim, scalingFct: Double; const palId,
        palcols: Integer; const asPolar, withWall: Boolean;
        const tiltAngleForColoring: Double);
    // evaluation functions
    procedure initSIOPorientation(const tilt,azi: Double);
    function evalSIOPWorld(const dir: Tvect3d; var frontSideHit: Boolean): Real;
    function evalSIOPWorldforDebugging(const dir: Tvect3d; var locOut: Tvect3d;
        var frontSideHit: Boolean): Real;
    procedure precalcHemiMeans(const forceRecalc: Boolean);
    function evalSIOPStdSYS(const ph, cosTh: Double): Real;
    // aux. functions
    procedure createCOSlambertSIOP(const sys: Integer);
    function dumpStats: String;
    function dumpSamples: String;
End;
```

All permanent data defining a *SIOP*, which is usually stored in a file, is structured in the `TSiopData` class. It covers the following information:

```
Type
TSiopData = record
    { TSiopData }
    // sampling data
    ssNDir: Integer; // number of directions (Fibo.hemi)
    sslsiType: TLSIEvent; // type of data transm., abs..
    ssInciType: TinciType; // for abs. front/back/both
    ssSIOPname: String[48]; // name for SIOP
    ssGeoName: String[48]; // target name
    ssRefArea: Double; // reference area (0=none)
    ssCalcTime: Double; // timestamp of generation
    ssRaysPerm2: Double; // rays per m³ used for sampling
    ssDat: array of Double; // sample data
    ssSemDat: array of Double; // standard error of means for samples

    // functional data
    shCoeffCnt: Integer; // number of coefficients
    shCoeffs: array of TshCoeff; // coefficient list
    shMeanErr: Double; // mean error vs. samples
    shAbsMeanErr: Double; // abs mean error vs. samples
    shMaxErr: Double; // max error vs. samples
```
```R
shHemiCosRelErr: Double; // cos. weighted relative error
shHemiRelErr: Double; // average relative error
shItCnt: Integer; // no. opt. iterations used for generation

// cached data of most recent evaluation
tiltAngleofMeans: Double; // tilt angle used for recent calculation
meanLowerHemi: Double; // SIOP avg. val gnd
meanUpperHemi: Double; // SIOP avg. val sky
meanHorizBand: Double; // SIOP avg. val horizon
solidAngleDivPIhorizBand: Double; // vis. solid angle horizon
solidAngleDivPIupperHemi: Double; // vis. solid angle sky
solidAngleDivPIlowerHemi: Double; // vis. solid angle gnd
End;
```

The class contains three types of information:

*   **sampling data:** the original sampled data, that is usually determined in the virtual measurement, i.e. raytracing process or alternatively using empirical methods. This original data is always kept within the *SIOP*. The information can be used for extended analysis or regeneration of the functional form. Further, supplemental information regarding the source of the sampling data is contained here (e.g. target object name, reference area, type of *SIOP* etc.).

*   **functional data:** this part contains all information required to construct/evaluate the functional form of the *SIOP*.

*   **cached evaluation data:** variables necessary for evaluating a *SIOP* are stored here (see section 6.2.4). Since this data is mainly generated using numerical integration, it is cached in the *SIOP*. Reusing values computed on earlier occasions increases the performance of later evaluations with the same tilt angle. The data contains the orientation and coefficients of the most recent evaluation.

To actually evaluate a *SIOP* in a simulation run, the evaluation has to be initialised by calling `initSIOPorientation` (to define the *SIOP* orientation in the real-world coordinate system) and `pre-calcHemiMeans` to initialise the constants used for evaluation (see next section).

After that, the function `evalSIOPWorld` will return the value of the *SIOP* for any desired direction:

```R
function TRCSIOPSamples.evalSIOPWorld(const dir: Tvect3d; var frontSideHit: Boolean): Real;
  // note for solar position: dir points INTO direction of sun
var d: Tvect3d;
  ph,cosTh: Double;
begin
  d.x:=VDot dir,rotMat[0]); // rotate world direction to SIOP system
  d.y:=VDot dir,rotMat[1]);
  d.z:=VDot dir,rotMat[2]);
  cosTh:=d.z;
  // if backside of SIOP hit, exit and return 0:
  if d.z<=0 then Begin Result:=0;frontSideHit:=false;Exit;End;
  frontSideHit:=true;
  ph:=arctan2(d.x,d.y);
  Result:=evalSIOPStdSys(ph,cosTh);
end;
```

The evaluation of the spherical harmonics series is performed in the function `evalSIOPStdSYS`:

```R
function TRCSIOPSamples.evalSIOPStdSYS(const ph, cosTh: Double): Real;
  // if cosTh_ is already provided in valid range [0,+1] it is used instead of th (for performance)
  // VALID RANGE for Hemispherical evaluation
  // PHI = [-PI,+PI] or [0,2*PI]
  // THETA = [0,+PI/2] or COSTH [0,+1]

var c: Integer;
  s: double;
begin
  if cosTh_ < 0 then // BACKSIDE of SIOP hit, exit and return 0
    Begin
      Result:=0;Exit;
    End;
  s:=0;
  with ssiop do
    for c:=0 to shCoeffCnt-1 do // evaluate sph.harm. series with non-zero coefficients
```
```
Begin
    s:=s+shCoeffs[c].coeff*spHarm.Ysh(shCoeffs[c].jId,ph,0,costh_);
end;
if s<0 then s:=0;
Result:=s;
end;
```

The coefficients of the spherical harmonics (equation (103)) are provided by the function `Ysh` function, which is implemented in the module `TSphericalHarmonics`.

The determination of the coefficients based on the simulated annealing optimisation method is performed by the class `TRCSIOPanneal` in the module `rc SIOPanneal`.

The method `anneal` is called with the relevant parameters to convert any given `SIOP` samples into the functional form.

```python
Type
    TRCSIOPanneal = class(TObject) { TRCSIOPanneal }
    private
        spHam : TsphereHarm;
        ssiop : TRCSIOPSamples;
        fiboHemi : TFibonacciSphere;
        maxL : Integer;
        coeffs : array of Double;
        logMemo : TMemo;
    procedure plotReconstruction(const uptoL: IntegerR; const objfn, lobename: String);
    procedure log(const tx : String);
public
    constructor Create(var logMemo_ : TMemo);
    function anneal(const fn: String; // file-name of SIOP file containing sample data
                        const ldegree: Integer; // max. allowed degree of spherical harmonics
                        const targetDuration: double; // targetDuration: 1 delivers decent results already, 10 is approx. one hour for HQ result
                        const addZeroBase: Boolean; // optionally force to 0 for grazing angles
                        const NameAppendix : String): String; // appended to output filename of SIOP file with results
    destructor Destroy; override;
End;
```
### 6.1.4. Exemplary generation of SIOPs

In order to allow a better understanding of the method and to prove its efficiency, an exemplary application is provided here. The application does not focus on the quantitative results, as this is covered in later chapters. Instead, this example should demonstrate the process of SIOP generation.

Four SIOPs describing the direct transmission of a fenestration detail comprising a triple-glazed insulated glass unit (IGU) shaded by external Venetian blinds, are generated. In the example, the state “unshaded window”, as well as the states “fully-shaded” with three different slat angles, are analysed. Based on the nature of the method, the evaluation also includes any reflections with the reveal or on the window sill. To generate

the measurement samples for the SIOP generation, a (here unspecified) raytracing algorithm is applied to a detailed CAD model of the fenestration (Figure 79). Provided a suitable measurement method is available, the samples could also be of empirical nature. In this example, 1444 equally spaced sampling directions of the Fibonacci-hemisphere (see section 5.3) are used to measure the effect of different parallel-beam incidence profiles on the physical quantity of interest. In the case presented here, the power transmitted to the internal of the building is analysed.

<figure><img src="image_84.png" /><figcaption>Figure 79 3D model of target used</figcaption></figure>

<figure><img src="image_85.png" /><figcaption>Figure 80 Results of raytracing – sampled data for total transmitted power (unshaded; shaded 0°; shaded 30°; shaded 60°)</figcaption></figure>

The sampled measurements are converted into SIOPs using the method described above. The annealing process for each sampling set, which determines the optimal reconstruction coefficients, lasts for about 15 minutes each (on a single average CPU). Each optimisation iterates through approximately

30 millions possible configurations. The available base of spherical harmonics has been restricted to 440 potential base functions by setting $l_{max}=20$.

The quality of the resulting reconstructions is quantified based on the absolute mean errors for the individual samples and the mean errors weighted by the (almost identical) solid angles for each sample. Geometrically this also corresponds to the volume of the SIOP if it is represented in polar plot representation (see Table 11), as the solid angles represented by the samples are almost equal (see section 5.3).

<table><caption>Table 11 Key figures of generated SIOPs</caption><thead><tr><th></th><th>unshaded<br>window</th><th>shaded<br>slats 0°</th><th>shaded<br>slats 30°</th><th>shaded<br>slats 60°</th></tr></thead><tbody><tr><td>number of<br>non-zero coeff.</td><td>79</td><td>83</td><td>82</td><td>89</td></tr><tr><td>mean abs. error<br>(samples)</td><td>0.15%</td><td>1.08%</td><td>1.28%</td><td>1.47%</td></tr><tr><td>(signed) mean error<br>(volume)</td><td>0.0002%</td><td>-0.002%</td><td>-0.002%</td><td>-0.0005%</td></tr><tr><td>SIOP avg. value<br>upper hemisphere</td><td>0.2988</td><td>0.1704</td><td>0.0727</td><td>0.0191</td></tr><tr><td>SIOP avg. value<br>lower hemisphere</td><td>0.3400</td><td>0.1817</td><td>0.2408</td><td>0.1378</td></tr><tr><td>SIOP reconstr.<br>quality R²</td><td>0.9992</td><td>0.9992</td><td>0.9990</td><td>0.9978</td></tr></tbody></table>