# RadiCal - A novel method to model energy flows induced by solar radiation

Daniel Rüdisser

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

To cite this version:

Daniel Rüdisser. RadiCal - A novel method to model energy flows induced by solar radiation. Optics

/ Photonic. ibpsc - Graz University of Technology, 2023. English. NNT: . tel-04183376

HAL Id: tel-04183376

https://hal.science/tel-04183376

Submitted on 19 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
<figure><img src="image_3.png" /></figure>

# RadiCal

A novel method to model energy flows
induced by solar radiation

Dipl.-Ing. Daniel Rüdisser

## DOCTORAL THESIS

to achieve the university degree of
**Doktor der technischen Wissenschaften**

submitted to the
**Faculty of Civil Engineering**
**Graz University of Technology**

Supervisor and 1st reviewer
**Univ.-Prof. Dipl.-Ing. Dr.techn. Michael Monsberger**

2nd reviewer
**Prof. em. John Mardaljevic. BSc PhD**
## AFFIDAVIT

I declare that I have authored this thesis independently, that I have not used other than the declared sources/resources, and that I have explicitly indicated all material which has been quoted either literally or by content from the sources used. The text document uploaded to TUGRAZonline is identical to the present doctoral thesis.

Date

Signature

## Acknowlegment Funding Agency

This work has been promoted by the Austrian Federal Ministry for Climate Action, Environment, Energy, Mobility, Innovation and within the "Forschungspartnerschaften / Industrienahe Dissertation" program under grant agreement No. 862523 ("RadiCal").

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

FFG
Forschung wirkt.

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

Bundesministerium
Klimaschutz, Umwelt,
Energie, Mobilität,
Innovation und Technologie

## Copyright / attribution ShareAlike License (CC BY-SA 4.0)

The copyright of this thesis rests with the author. Unless otherwise indicated, its contents are licensed under a [Creative Commons Attribution-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-sa/4.0/) (CC BY-SA).

Under this licence, you may copy and redistribute the material in any medium or format for both commercial and non-commercial purposes. You may also create and distribute modified versions of the work. This on the condition that: you credit the author and share any derivative works under the same licence.

When reusing or sharing this work, ensure you make the licence terms clear to others by naming the licence and linking to the licence text. Where a work has been adapted, you should indicate that the work has been changed and describe those changes.

Please seek permission from the copyright holder for uses of this work that are not included in this licence or permitted under applicable copyright law.
Abstract

In light of global warming, it is essential to increase the energy efficiency of buildings and mitigate
the adverse effects of rising temperatures on building users. Solar-induced energy flows are a signifi-
cant, often even dominant, component of the building's energy balance, having a crucial impact on
the cooling and heating energy demand as well as on thermal comfort. However, in contrast to the
detailed and elaborate calculation methods applied to describe heat transmission through the building
envelope, the currently used methods to model solar-induced energy flows lack accuracy. The inac-
curacies result from the fact that the related energy flows are linked to a large number of complex
and multifaceted optical processes. The simplified models currently used for this purpose are based
on a high level of abstraction, involving potentially significant errors.

In order to overcome these shortcomings, a novel method for modelling solar-induced energy flows
for building applications is developed, implemented and tested. A fundamentally different approach
is proposed based on the analysis of the approaches and limitations of the currently applied methods.
It allows for deriving all energy-related transmission and absorption characteristics of target objects,
e.g. windows, façade elements or solar collectors, based on a detailed three-dimensional model in a
comprehensive way, considering the entire global radiation spectrum, ranging from ultraviolet to
near-infrared. In order to be able to model the energetically relevant optical processes accurately,
advanced models of physical optics, primarily applied in material science and astronomy, are adopted
for their use in the building science context. The models are integrated into a newly developed Monte
Carlo raytracer. In the algorithm, light propagation is still modelled based on the principles of geo-
metrical optics, while scattering processes, i.e. reflection, absorption and transmission, are modelled
based on the relevant electrodynamic models. The approach requires considering and tracking light's
polarisation state and spectral distribution throughout all processes. Applying physical optics allows,
in particular, for accurate modelling of coated glass surfaces. This is essential, as modern glazings
rely on thin-film coatings to improve their thermal and optical properties.

Beyond the physical optics approach, which fundamentally relies on solving Fresnel's equations in
their complex-valued form, a statistical roughness model and a generic subsurface-scattering model
allow the modelling of practically relevant surfaces. In order to solve the complex system of equations
efficiently, Monte Carlo techniques are applied throughout the entire implementation.

In addition to the novel method proposed for optical modelling, a new interface for deploying the
detailed optical data gathered in the raytracing process to building performance simulation is devel-
oped. The method is based on spherical harmonics and effectively enables simulation tools to perform
solar calculations with raytracing accuracy, but virtually at no computational cost.

Many validation steps and test cases are performed at different levels of modelling to prove and
demonstrate the validity, not only of the models and equations but also of their algorithmic imple-
mentation. At the highest level, the entire workflow is applied to perform virtual measurements of a
triple-glazed window equipped with four different shading devices. Real-world measurements per-
formed over several months with varying solar incidence profiles are compared with virtual measure-
ments performed on the modelled window. A dedicated measurement device and method were devel-
oped for this purpose.

The proposed alternative solar-modelling method could be successfully developed, implemented,
tested and applied in the course of the PhD. It can be used for a wide range of research and industrial
applications in the field of building performance simulations and beyond. Its full potential will be
harnessed and explored in subsequent studies.
# Kurzfassung

Vor dem Hintergrund des Klimawandels gilt es Gebäude so zu gestalten, dass sie einerseits energieeffizient sind und andererseits die Gebäudenutzer:innen vor der zunehmenden Hitze schützen. Eine besondere Bedeutung kommt hierbei solarinduzierten Wärmeströmen zu. Sie stellen eine wichtige, vielfach die dominante, Komponente der Gebäudeenergiebilanz dar und wirken bestimmend für den Kühlenergiebedarf, Heizenergiebedarf sowie den thermischen Komfort. Im Gegensatz zu den genauen und detaillierten Modellen, welche zur Beschreibung der Energietransmission durch die Gebäudehülle angewendet werden, sind die Methoden zur Berücksichtigung der Solarstrahlung oft mit hohen Ungenauigkeiten verbunden. Energieströme, welche durch Sonneneinstrahlung hervorgerufen werden, resultieren aus einer Vielzahl von komplexen, optischen Prozessen, welche die existierenden vereinfachten Modelle nur unzureichend abbilden können.

Um dieser Problemstellung zu begegnen, wird in der vorliegenden Dissertation eine neue Methode zur Modellierung von solarinduzierten Energieströmen entwickelt, implementiert und getestet. Die Methode ist als Alternative zu bestehenden Modellen zu sehen. Sie gestattet es, die energierelevanten Transmissions- und Absorptionseigenschaften von Gebäudekomponenten (Fenster, Fassadenelemente, Solarkollektoren) auf Basis detaillierter, dreidimensionaler Modelle und unter Berücksichtigung des gesamten solaren Spektrums zu ermitteln. Um die optischen Prozesse adäquat abbilden zu können, wurden moderne Modelle aus dem Gebiet der physikalischen Optik, welche aktuell vorwiegend in der Materialwissenschaft und Astronomie eingesetzt werden, ausgewählt und für den effizienten bauphysikalischen Einsatz adaptiert. Ein neu entwickelter Monte-Carlo-Raytracing-Algorithmus bildet die Plattform für die Anwendung der Modelle. Die Lichtausbreitung wird hierbei weiterhin gemäß Prinzipien der geometrischen Optik modelliert, die Streuprozesse werden jedoch vorwiegend auf Basis der relevanten elektrodynamischen Modelle abgebildet. Hierfür ist es notwendig, die Polarisation und die spektrale Verteilung des Lichtes zu berücksichtigen. Zentrale Bedeutung kommt der Lösung der Fresnel'schen Gleichungen in deren komplexwertigen Form zu. Der gewählte Ansatz ermöglicht insbesondere die Modellierung beschichteter Glasoberflächen. Dies ist essentiell, da bei modernen Verglasungen Dünnfilmtechniken zur Verbesserung der thermischen und optischen Eigenschaften eingesetzt werden. Neben den elektrodynamischen Modellen gestattet ein statistisches Oberflächenrauheitsmodell und ein generisches Subsurface-Streumodell die Abbildung praktisch relevanter Oberflächen. Zur effizienten Lösung der resultierenden Gleichungssysteme werden vorwiegend Monte-Carlo-Methoden angewendet.

Neben der neuartigen optischen Modellierung beinhaltet die Dissertation eine neue Methode zur effizienten Weitergabe der detaillierten Raytracing-Ergebnisse an thermische Simulationsplattformen. Die Methode basiert auf Kugelflächenfunktionen und gestattet Drittanwendungen, solare Energieströme mit Raytracing-Genauigkeit zu berücksichtigen, jedoch ohne den damit verbundenen Rechenaufwand. Die Korrektheit der Modelle und Implementierung wird mit vielen Validierungsbeispielen auf unterschiedlichen Ebenen demonstriert. Schließlich wird der gesamte Workflow zur Modellierung eines beschatteten und nicht-beschatteten Fensterdetails mit Dreifachverglasung eingesetzt. Hierbei werden Messungen am realen Fenster nachgebildet, welche über einen Zeitraum von mehreren Monaten, unter unterschiedlichsten Einstrahlungsbedingungen und mit vier unterschiedlichen Sonnenschutztypen durchgeführt wurden. Für diesen Zweck wurde eigens eine neuartige Messeinrichtung und Messmethodik entwickelt.

Die vorgeschlagene alternative Methode zur Modellierung solarer Energieströme konnte erfolgreich entwickelt, umgesetzt und getestet werden. Sie gestattet eine Vielzahl von Anwendungsmöglichkeiten, deren volles Potential in zukünftigen Arbeiten erst ausgelotet werden muss.
*We can scarcely avoid the inference that light consists in the transverse undulations
of the same medium which is the cause of electric and magnetic phenomena.*

James Clerk Maxwell, 1862

**Acknowledgements**

I would like to acknowledge the Austrian Funding Agency FFG and the Austrian Ministry for Transport, Innovation and Technology for funding large parts of this work.

I would like to express my gratitude to Prof. Michael Monsberger and Prof. John Mardaljevic for supervising and reviewing this dissertation. They always provided patient, kind and constructive support, even though the work significantly exceeded the initially targeted scope.

I would also like to extend a warm thank you to Prof. Christina Hopfe, Simon Weber and Kurt Iskra for proofreading many pages and providing valuable input. Coincidentally, Kurt already acted as my thesis supervisor many years ago.

I would further like to thank Anton Hammer and Robert Peißl of WO&WO Sonnenlichtdesign GmbH & Co KG for providing shading devices and their knowledge in a very unbureaucratic and kind way. Special thanks go to the measurement department of AEE INTEC, namely Waldemar Wagner, Roland Kerschenbauer and Reinhard Pertschy, for actively supporting and enabling my ambitious measurement ideas. Many thanks also go to my former "dissertation peer" and later boss Tobias Weiss and the entire team of AEE INTEC.

Last but not least, I want to thank my parents, my wife Marlies, and my daughters Emma and Martha for their continued love and support.
# CONTENTS

1. Introduction.................................................................................................................. 1

1.1. Motivation.................................................................................................................. 1

1.2. Research question and objective........................................................................ 3

1.3. Overview of thesis structure.................................................................................. 4

1.4. Overview of implemented modules........................................................................ 6

1.5. Scope........................................................................................................................ 7

2. Main methodological approaches................................................................................ 8

2.1. Physical optics approach........................................................................................ 8

2.2. Fallback approach – generic models........................................................................ 9

2.3. Solar incidence operators - interface to other methods........................................ 9

2.4. Validation.................................................................................................................. 10

2.5. Separation of the optical and thermal domain........................................................ 11

2.6. Implementation........................................................................................................ 12

2.7. Programming language ...................................................................................... 13

3. State-of-the-art analysis.................................................................................................. 14

3.1. General model and definitions............................................................................ 14

3.2. Modern approaches - BSDF and DSHGC ................................................................ 15

3.3. Limitations of BSDF and DSHGC approaches........................................................ 22

3.4. Approaches of international standards.................................................................... 25

3.5. Limitations of approaches in international standards................................................ 28

4. Light, sun and optics - applied principles, models and methods.............................. 29

4.1. A brief history of light .......................................................................................... 29

4.2. Nature of light ........................................................................................................ 31

4.3. Light-matter interaction overview........................................................................ 40

4.4. Global radiation........................................................................................................ 50

4.5. Solar position algorithm........................................................................................ 55

4.6. Complex-valued index of refraction functions........................................................ 58

4.7. Generalized Fresnel coefficients............................................................................ 62

4.8. Stokes vectors ........................................................................................................ 67

4.9. Müller matrix operations........................................................................................ 70

4.10. Roughness model and microfacet theory............................................................ 80

4.11. Subsurface scattering........................................................................................ 89

4.12. Inversion method to derive complex-valued refractive index functions of transparent materials ................................................................................................................ 93

4.13. Inversion method for a generic model of diffusely reflecting, opaque surfaces ............ 105

5. Additionally applied methods...................................................................................... 110

5.1. Monte Carlo method.................................................................................................. 110

5.2. Simulated annealing.................................................................................................. 114

5.3. A Fibonacci hemisphere........................................................................................ 117

6. Solar Incidence Operators Definition, generation and evaluation.............................. 123
<table>
  <tbody>
    <tr>
      <td>6.1. SIOP - Solar Incidence Operator</td>
      <td>123</td>
    </tr>
    <tr>
      <td>6.2. SIOP evaluation and incidence profiles</td>
      <td>134</td>
    </tr>
    <tr>
      <td><strong>7. A full spectral polarisation Monte-Carlo Raytracer</strong></td>
      <td><strong>154</strong></td>
    </tr>
    <tr>
      <td>7.1. Tracing a ray - overview</td>
      <td>154</td>
    </tr>
    <tr>
      <td>7.2. Rays and light particles</td>
      <td>155</td>
    </tr>
    <tr>
      <td>7.3. 3D model and geometry modules</td>
      <td>156</td>
    </tr>
    <tr>
      <td>7.4. Collision detection</td>
      <td>159</td>
    </tr>
    <tr>
      <td>7.5. Periodic planes</td>
      <td>163</td>
    </tr>
    <tr>
      <td>7.6. Smooth shading</td>
      <td>164</td>
    </tr>
    <tr>
      <td>7.7. Material-Interface class</td>
      <td>165</td>
    </tr>
    <tr>
      <td>7.8. Scattering classes (TLightSurfaceInteractionStokes)</td>
      <td>167</td>
    </tr>
    <tr>
      <td>7.9. The TLSISroughPol scattering class</td>
      <td>168</td>
    </tr>
    <tr>
      <td>7.10. Directional scanning process</td>
      <td>172</td>
    </tr>
    <tr>
      <td>7.11. Rendering and backward raytracing</td>
      <td>181</td>
    </tr>
    <tr>
      <td><strong>8. Full system empirical validation</strong></td>
      <td><strong>185</strong></td>
    </tr>
    <tr>
      <td>8.1. PyroScanner measurement method</td>
      <td>185</td>
    </tr>
    <tr>
      <td>8.2. Results</td>
      <td>204</td>
    </tr>
    <tr>
      <td><strong>9. Application case</strong></td>
      <td><strong>216</strong></td>
    </tr>
    <tr>
      <td>9.1. Target object</td>
      <td>216</td>
    </tr>
    <tr>
      <td>9.2. Model and workflow</td>
      <td>217</td>
    </tr>
    <tr>
      <td>9.3. Evaluation with temporally highly resolved data</td>
      <td>219</td>
    </tr>
    <tr>
      <td>9.4. Evaluation of energetic performance</td>
      <td>220</td>
    </tr>
    <tr>
      <td>9.5. Determination of effective performance parameters</td>
      <td>221</td>
    </tr>
    <tr>
      <td><strong>10. Discussion</strong></td>
      <td><strong>223</strong></td>
    </tr>
    <tr>
      <td>10.1. State-of-the-art methods and novel method</td>
      <td>223</td>
    </tr>
    <tr>
      <td>10.2. Physical optics integration</td>
      <td>224</td>
    </tr>
    <tr>
      <td>10.3. Required material data and inversion approaches</td>
      <td>225</td>
    </tr>
    <tr>
      <td>10.4. Results of validation steps</td>
      <td>226</td>
    </tr>
    <tr>
      <td>10.5. Limitations</td>
      <td>227</td>
    </tr>
    <tr>
      <td>10.6. Integration into current frameworks</td>
      <td>228</td>
    </tr>
    <tr>
      <td>10.7. Discussion with respect to research questions</td>
      <td>228</td>
    </tr>
    <tr>
      <td><strong>11. Conclusion</strong></td>
      <td><strong>230</strong></td>
    </tr>
    <tr>
      <td>11.1. Summary</td>
      <td>230</td>
    </tr>
    <tr>
      <td>11.2. Key findings</td>
      <td>231</td>
    </tr>
    <tr>
      <td>11.3. Further directions</td>
      <td>234</td>
    </tr>
    <tr>
      <td><strong>List of figures</strong></td>
      <td><strong>235</strong></td>
    </tr>
    <tr>
      <td><strong>List of tables</strong></td>
      <td><strong>239</strong></td>
    </tr>
    <tr>
      <td><strong>References</strong></td>
      <td><strong>240</strong></td>
    </tr>
    <tr>
      <td><strong>Appendix A: Full-system validation results</strong></td>
      <td><strong>250</strong></td>
    </tr>
  </tbody>
</table>
# 1. Introduction

## 1.1. Motivation

In light of global warming and considering that in developed countries, buildings account for almost half of the total energy consumption (U.S. Department of Energy, 2015; European Commission Energy dept., 2020), it is essential to increase the energy efficiency of buildings. Solar-induced energy flows through the building envelope, particularly solar inputs through transparent surfaces, are highly relevant in this context. Due to the constantly rising proportions of glazed surfaces in modern buildings, the significance of this energy component has dramatically increased in recent decades. Solar gains through transparent surfaces can considerably lower the demand for heating and lighting energy in winter (Zirnhelt et al., 2015; Marino, Nucara et al., 2017a; Ménard and Souviron, 2020) while often representing the predominant heat source in summer. Therefore, solar inputs are also directly linked to increased cooling energy demand, indoor discomfort or potentially hazardous overheating (Hamdy et al., 2017; Marino et al., 2017b; McLeod et al., 2017).

In order to increase the energy efficiency of buildings, it is essential to have suitable models for solar-induced energy flows available. Accurate and physically valid models can, on the one hand, be used to develop and optimise the energetic performance of glazings and shadings at the component level. On the other hand, they are essential to enhance the accuracy of building performance simulation (BPS) for the entire building. BPS is applied to improve the energy efficiency of buildings on various levels ranging from the early design phase over design development to operation. Hence, the quality of the models used to determine solar energy flows has impacts on many levels ranging from cooling and heating capacity design over energy demand determination and optimisation to the assessment of the building's user comfort.

However, in contrast to the detailed and elaborate calculation methods currently applied to describe heat transmission through the building envelope, the currently used methods to model solar-induced energy flows lack accuracy. The inaccuracies can be attributed to the fact that the related energy flows are linked to a large number of complex and multifaceted optical processes. Addressing this challenge with simple models requires a high level of abstraction that involves potentially significant errors.

Performing full-system measurements of windows, shading devices or façade elements instead of using models is no viable alternative. The available empirical methods are time-consuming and costly on the one hand and can only be carried out for specific, individual configurations on the other hand. Due to the variety and a large number of potential combinations of materials and systems, it is virtually impossible to cover all relevant products empirically.

The currently most applied methods to model the energetic performance of shaded or unshaded glazed components have been laid down in international standards. These calculation methods are based on severe simplifications. E.g., the commonly used Venetian blinds are still modelled assuming indefinitely thin, planar and ideally diffuse reflecting elements (ISO 15099, 2003; ISO 52022 part 1, 2018), performance values of glazings are calculated on the assumption of normal incidence and significant angular and spectral dependencies of light scattering processes are not considered (ISO 9050, 2003; EN 410, 2011). Based on these simplifications, calculating solar gains of shaded and unshaded glazings can involve significant errors (Platzer, 2001; Kuhn, 2006; Wermke, 2021).

Research performed in the last two decades focused mainly on directional dependence, which was disregarded in previous approaches, and led to novel, more complex approaches. The most commonly
used modern methods to determine the optical properties and energetic performance of glazings rely on scientific models that are implemented in the *Window* software by Lawrence Berkeley National Laboratory (LBNL Windows and Daylighting Group, 2022b). The method allows excellent modelling of multi-layered glazing components based on an extensive database containing product-specific, measured spectral data for normal incidence. However, it utilizes two simple models to approximate the directional behaviour for incidence at oblique angles (Curcija *et al.*, 2018) and features only a few specific and simplified models for considering shading devices. In order to model shading elements with higher accuracy, additional raytracing software, mostly *Radiance* (Radiance Community, 2022), is used to provide directional scattering information (BSDF) in the form of Klems' matrices (Klems, 1994b). The resulting combined optical model is then commonly applied in an energy balance model to derive directional solar heat gain coefficients (DSHGC); see, e.g., the two-layer model (Kuhn *et al.*, 2011; Bueno *et al.*, 2017). The derived DSHGC values can finally be used in energy performance simulation tools to determine the total transmitted power based on solar irradiance parameters. The modelled heat gains comprise directly transmitted radiation and secondary, transmissive heat flows. The entire workflow is complex, and its fragmentation causes inaccuracies that can be divided into two groups: on the one hand, the layered BSDF approach involves substantial limitations regarding its optical accuracy, particularly regarding spectral and angular selectivity, spatial information and polarisation. On the other hand, the derivation of DSHGC values requires the assumption of predefined thermal boundary conditions and thermal properties of all components. These parameters are mostly assumed as constant. However, they are, in fact, dynamic by nature and depend on the results of the final building simulation.

As a consequence of the limitations of the existing methods for modelling solar gains of transparent surfaces in the building envelope, the calculated energy performance of the entire building may be significantly inaccurate. This involves potentially wrong planning decisions in all phases of the building's life cycle, ranging from early design to operation. Beyond that, the limitations of the currently applied models have adverse impacts on the evaluation and optimisation of existing shading or glazing elements and the development of innovative, new systems.

The scientific work in the recent two decades has mostly focused on combining and enhancing the approaches described above and, in particular, on providing ways to integrate them into existing building performance simulation methods. In recent years, dedicated, stand-alone tools to dynamically simulate energy demand and daylighting aspects based on the approaches discussed above were added by research teams working in this field, see the *Comfen* tool by LBNL (Mitchel *et al.*, 2019) or the *Fener* tool (Bueno *et al.*, 2015). However, the fragmented DSHGC workflow and the underlying fundamental optical models, which mostly date back to the eighties, are rarely questioned, representing a current research gap. On the one hand, a more universal modelling approach must be able to consider the three-dimensional structure of building components with a high level of detail. On the other hand, it has to be able to model light scattering processes with the precision of microscopic laboratory tests, considering the spectral, angular and polarisation properties of the incidence light and the target. This high level of modelling is required to accurately model the interaction of light with objects combining coated or uncoated transparent and opaque materials. The thesis's main target is to develop a universal and unified method that meets these requirements to allow for accurate modelling of solar-induced energy flows for building science applications. The proposed method is based on a more fundamental level of light modelling, but is still able to provide its functionality in a comprehensible, robust and practical way.
## 1.2. Research question and objective

Based on the issues raised in section 1.1 and further discussed in the state-of-the-art analysis (chapter 3), the following main research question and four sub-research questions are formulated and addressed in this work:

1. How can the accuracy of modelling solar-induced energy flows in building science be enhanced?

1.1. Which optical processes, beyond simple geometrical optics, are relevant, and how can they be considered? (e.g. thin-film interference, polarisation, surface roughness)

1.2. Is it possible to create a physical-optics-based raytracer that can be used to efficiently evaluate the solar properties of relevant building components, in particular combinations of glazings and shading devices?

1.3. Can the raytracer be designed in a way that allows a seamless transition or downgrading to existing models?

1.4. If and how can the highly accurate raytracing be deployed to practical application in dynamic building performance simulation?

In order to answer these questions, the work was subdivided into subtasks. The subtasks are addressed by following these primary objectives:

### Objective 1: state of the art analysis

In order to understand how researchers and practitioners currently perform the modelling of solar-related energy flows in building science, the current research and the international standards in force are reviewed. The impact and root causes of the limitations associated with the current models and workflows are analysed. Further, the relevance of optical processes and advanced physical optics models are analysed, and their adaption for use in building science is evaluated.

### Objective 2: development of a physical-optics-based raytracer

Based on the knowledge and requirements gained in the first task, a suitable novel raytracing method is designed and implemented. The method shall be able to consider practically relevant effects of physical optics and cover the entire global radiation spectrum. Furthermore, the raytracer shall be implemented in an open design approach based on object-oriented programming methods in order to facilitate the seamless addition of new scattering and material models in the future. Further, the raytracer shall optionally be able to apply existing simple models in order to establish a consistent transition to current modelling approaches.

### Objective 3: development of a method to provide raytracing results for building performance simulation

A method to deploy the information contained in the accurate and detailed results of the raytracing analysis shall be established. The method should allow thermal simulation methods, particularly dynamic building performance simulation tools, to exploit the raytracing results. While the determination of the scattering (i.e. absorption and transmission) characteristic of the analysed target object requires an elaborate and potentially multi-hour-long raytracing process, its results shall be provided in a form that allows a fast and accurate evaluation. This is crucial for dynamic simulation, as sub-hourly evaluation steps for multi-month-long periods are required in building performance simulation, resulting in a large number of necessary evaluations.
**Objective 4: proof of concept and validation**

Due to its complexity and novelty, the developed method's validity and feasibility, as well as its correct implementation, must be demonstrated transparently and comprehensively. For this purpose, several validation examples and application cases shall be performed. These validation steps shall be performed at different levels: physical optics modelling, application for building simulation and a full-scale empirical validation of the developed workflow.

**1.3. Overview of thesis structure**

<figure><img src="image_6.png" /><figcaption>Figure 1 Graphical overview of the thesis structure</figcaption></figure>

Figure 1 provides a graphical overview of the thesis chapters with regard to the objectives formulated in section 1.2. Chapter 2 highlights the main methodological approaches followed in the course of the thesis. Due to the nature this PhD, the subsequent chapters also contain essential methodical elements. As the thesis comprises multiple, separate theoretical and empirical components (see Figure 2), it seemed more appropriate to include theory, methods, implementation and tests or validation as sub-chapters within the relevant chapters.

Chapter 3 discusses the current state-of-the-art. It focuses on the currently applied methods and models comprising modern approaches and the methods specified in international standards. The chapter also covers a discussion of the limitations that arise from the models used.

Chapter 4 covers all methods and models that are in this PhD used for modelling solar radiation, propagation of light and scattering of light. As mentioned above, the chapter contains not only the models but also the relevant theory, their implementation in the context of the developed method, as well as validation and testing results. In order to address the complexity of the field of optics and light, the chapter also provides a broader view of physical optics and discusses the relevance of scattering processes that are not part of the final method. Providing this broader context is essential to show the relevance and limitations of the applied light scattering models. The chapter starts with
a brief summary of the history of light. It primarily addresses objective 2 of the thesis, but also
objectives 1 and 4.

Chapter 5 covers three essential components of the developed method that 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. The chapter addresses objectives 2 and 4.

Chapter 6 addresses objective 4 by introducing a novel concept developed in the thesis and referred to as *solar incidence operator (SIOP)*. The operator represents an efficient method to deploy raytracing results to other applications, particularly dynamic building performance simulation. The first part of the chapter introduces the definition of the operator, whereas the second part focuses on its application, i.e. its evaluation based on different irradiance models.

In chapter 7, the physical models and methods of chapter 4 and the stochastic methods of chapter 5 are applied to establish a full-spectral Monte-Carlo raytracer. The raytracer is then used in an angular scanning mode to measure and extract all relevant optical characteristics of target objects (e.g. glazings) required to derive the corresponding accurate energetic models.

In chapter 8, the entire workflow of the developed method is validated against empirical data. The transmittance of a shaded triple-glazed window has been measured for several months in various configurations. The developed method is then applied to a detailed CAD model of the window to derive the corresponding *solar incidence operators*. In a digital-twin approach, these *SIOPs* are subsequently used to perform virtual measurements in the model. The computed transmittances are finally compared against the measured values.

Chapter 9 demonstrates the potential and performance of the method based on a typical application of dynamic thermal building performance simulation. A shaded window's energetic key performance figures are determined for various orientations and locations. All relevant energy flows are determined using a transparent simulation model designed explicitly for that purpose.

Important aspects of the developed method, its potential and limitations as well as the results of the
performed validations are discussed in chapter 10.

Chapter 11 finally provides a summary, key findings and an outlook regarding future work.
## 1.4. Overview of implemented modules

The method proposed here is called *RadiCal* as an acronym for the term *radiation calculator* but also to indicate its disruptive (radical) innovation character. The *RadiCal* method as well as the work performed in the course of the PhD, represent a framework comprising different functional modules. While the intended workflow of the method relies on most of these modules, the components are largely self-contained and can also be used independently or integrated into other frameworks. Figure 2 provides an overview of the main components implemented in this thesis. For the sake of clarity, many additional modules implemented for the purpose of validation, testing, import and export were omitted in this overview.

<figure><img src="image_7.png" /><figcaption>Figure 2 Overview of implemented modules and addressed objectives</figcaption></figure>

The dash-outlined areas indicate self-contained functional modules. The majority of modules were implemented for the physically based raytracer (covered in chapters 4 and 7). The main result of the raytracer scanning process is directionally sampled data. This data is stored and converted into a solar incidence operator *SIOP* (see section 6.1). The *SIOP* operator can be evaluated based on any provided irradiance model (see section 6.2). Finally, the feasibility and accuracy of the method are
demonstrated in a full-system validation case (chapter 8) and in an application case (chapter 9). The
measurement device and process developed for the validation step in the course of the PhD are also
covered in chapter 8.

1.5. Scope

The *RadiCal* method developed in the course of this PhD aims to provide an alternative approach that can overcome the shortcomings of currently applied methods to model solar-induced energy flows in field of building science. The novel method is theoretically founded, designed, implemented, tested and validated. It aims to provide a unified and enhanced method for modelling energy flows induced by solar radiation. In order to overcome the limitations of current methods, it is necessary to follow a more general, higher-level approach that allows for more specialized modelling. This can be achieved if light's electrodynamic nature is considered by applying models of physical optics. At this fundamental level, relevant and seemingly complex processes, which state-of-the-art methods cannot adequately address, arise naturally, e.g. thin-film interference of coatings, angular dependencies or polarisation effects.

However, beyond diving deeper into the physics of light, the method must also meet practical re-
quirements enabling its use in building science frameworks. In particular, this means that on the
"input side", it has to be able to import comprehensive geometries in standard CAD format effi-
ciently. On the "output side", a straightforward way to integrate the accurate optical information
determined within a raytracing process into thermal simulation frameworks is required. In order to
achieve this, a novel operator is introduced. It allows the deployment of comprehensive optical data
to other applications, particularly for dynamic building performance simulation. The operator's func-
tionality enables simulation tools to compute solar-related energy flows with raytracing accuracy at
virtually no computational costs.

In order to prove the functionality, validity and usability of the novel method, a full software imple-
mentation of the *RadiCal* method is required. This step is essential to demonstrate the feasibility
and validity of the developed method. The implemented modules are first applied to validation and
test cases at different modelling levels. Afterwards, the method is validated against full-scale meas-
urements on a typical shaded window with triple glazing. A dedicated measurement device and
method were developed for this purpose as part of the PhD. Finally, a typical building simulation
application case demonstrates the method's efficiency.

This PhD focuses on the energetic modelling of active and static shading systems, façades and glaz-
ings, or, more generally, on calculating all intended and unintended solar-related energy flows rele-
vant to buildings. However, the *RadiCal* method and the developed software tool can be used in
applications beyond that field, e.g. for solar thermal collectors or photovoltaics and, in particular,
also for lighting. Almost all methods and models applied here are relevant to the field of natural and
artificial lighting. Visible light is, of course, nothing but a distinct range of the entire global radiation
spectrum. Since the *RadiCal* method can consider distinct wavelengths or spectral bands, it is well
suited for modelling and studying lighting-related processes. Hence, the method also provides a novel
approach to addressing challenging topics in research and practical work in daylighting and artificial
lighting. However, due to the considerable effort required to cover the energy-related modelling com-
prehensively, any daylighting aspects had to be excluded from this work. It will, however, be inter-
esting to exploit this potential in subsequent research.
## 2. Main methodological approaches

This chapter gives an overview of the main methodological approaches and fundamental principles pursued in this work. Since the work covers elements from various fields, specific methodologies are elaborated further in the relevant chapters.

### 2.1. Physical optics approach

The term physical optics is generally used to refer to optical effects that cannot be explained based on simplified geometrical optics. Alternatively, perhaps more appropriately, the term wave optics is used, as this refers to the nature of light as an electromagnetic wave. The most fundamental approach to optics would require solving Maxwell's equations with high spatial and temporal resolution. Since this level is usually not achievable for practical applications, a significant degree of abstraction and simplification is necessary. The methods applied in geometrical optics are, in fact, a result of far-reaching abstraction. The present work attempts to improve the level of abstraction by enriching geometrical optics with the most relevant models of physical optics in a comprehensible and practically relevant way. The more fundamental level of physical optics allows more accurate models that potentially require significantly fewer input data than empirical models. Beyond that, more fundamentally based models are also more likely to meet future, yet unknown, requirements.

The field of physical optics naturally covers a wide range of effects ranging from non-linear optics at high field strengths in laser physics, over interferometry or Doppler-shifts in astronomy to nano-optics in material science. The present work identifies the processes with the most significant impact on solar-related energy flows. For this purpose, suitable models are adapted for implementation. Explicitly, these are:

* consideration of global radiation spectrum and spectral dependencies

* consideration of polarisation state

* consideration of fundamental material properties (refractive indices)

* consideration of angular and spectral dependencies

* consideration of surface roughness and subsurface scattering

* consideration of coherence effects on thin-flims (coatings)

* consideration of attenuation in bulk materials

The reason to follow this fundamental approach is two-fold. First, only physically based models are able to capture a satisfyingly wide range of practically relevant optical phenomena. Due to the intrinsically complex nature of light-matter interactions, simple models are only valid within strictly limited boundaries. Consequently, there is a constant need to apply additional ad-hoc models or workarounds to cover effects beyond the scope of the simplified model (e.g. wavelength dependences, polarisations effects etc.). Second, physical models require significantly fewer input parameters to achieve the same level of accuracy as empirical, data-driven models, as essential information on the processes is contained within the models in the form of fundamental relations. Therefore even limited inputs or constant approximations can potentially provide spectrally resolved results with higher accuracy than the models used in international standards.
## 2.2. Fallback approach – generic models

This PhD aims to provide models that are capable of capturing light scattering phenomena at a more fundamental physical level compared to geometrical optics. The implemented physical optics models mostly rely on fundamental physical quantities, particularly on the complex-valued index of refraction. Even though these quantities are accessible to measurement, they might not be readily available for many materials. In order to still be able to use the method with insufficient data, a fallback approach is pursued. The implementation of fallback cases allows a gradual downgrade of the models, finally ending at the level of common geometrical optics. In the case of missing material information, an inversion approach can be used to derive the information based on available empirical data (see section 4.12). If, however, the required spectrally resolved data is unavailable, the implementation of the fallback approach allows a downgrade to simpler, generic models. For example, if a material's spectral refractive index function is unavailable, a reasonable constant value for the global radiation range can still be applied to achieve spectrally resolved results. If this approach seems inappropriate, the refractive index function alternatively can be set to unity. In this case, the subsurface reflectance function can be used to model ideal diffuse reflection. Again, any available spectrally resolved information can be used in this model, either highly resolved or integrated for specific wavelength ranges (e.g. visible, UV, IR). If even these spectrally resolved reflectance values are unavailable, a constant diffuse reflectance value can be assumed for the entire spectrum. This final and most simple generic model corresponds to the geometrical optics model that most standard approaches are based on.

## 2.3. Solar incidence operators - interface to other methods

Beyond the application of physical optics, the concept of solar incidence operators (SIOP) is an essential component of the novel method presented here. SIOPs allow accurate determination of a scalar target quantity from any given solar incidence profile. Most often, this quantity will be the power transmitted through the target object or absorbed by specific parts of it. However, SIOPs can be generated for any quantities that are induced by solar radiation and can be described by a scalar value (e.g. reflected power, degree of polarisation, luminous exitance, gloss etc.). The formulation of the SIOP relies on spherical harmonics commonly known from quantum mechanics ("orbitals"). In the SIOPs, three-dimensional scattering information gathered in the raytracing process can be stored with a high level of data compression. SIOPs can finally be evaluated in a straightforward way by third-party tools, such as dynamic building simulation, allowing them to perform solar radiation calculations with raytracing accuracy but virtually at no computational costs.
<figure><img src="image_8.png" /><figcaption>Figure 3 RadiCal application workflow -</figcaption></figure>

Figure 3 gives an overview of how the developed *RadiCal* method is applied from a user's point of view:

1. A detailed 3D model of the target to be analysed is imported in a CAD format (e.g. a shaded window or solar collector).

2. Physically based material models are assigned to all surfaces of the model.

3. The raytracer scanning process is carried out with the desired angular resolution and precision.

4. The targeted quantities for all relevant components (e.g. transmission, absorption, reflection) are virtually measured for all sampling directions.

5. The sampled data is converted into the functional SIOP form.

6. The SIOP data containing the detailed optical scattering information of the target is deployed to third-party simulation tools.

## 2.4. Validation

Throughout the implementation, many validations and testing steps are performed to verify and demonstrate the validity of the implemented algorithms. Validation is given high priority for several reasons. First, the implemented algorithms are complex and beyond the standard algebra usually applied. Therefore they require thorough verification. Apart from the extensive three-dimensional vector calculations, four-dimensional vectors and matrices are applied, and a large part of the calculation is based on complex-valued algebra. Beyond that, the main algorithms rely on the stochastic Monte Carlo approach, which also adds complexity to testing and validation. Finally, validation cases are essential to demonstrate the validity and suitability of the different modules and workflow in a comprehensive way. The field of physical optics is a rather special topic. Therefore, the availability of suitable validation data is limited. For this reason, published results of fields like astronomy and material science were used to demonstrate the functionality and validity of the implemented algorithms in the first step. The validation cases contained in section 4.9.3, are considered most important as they show the validity of the algorithms with high accuracy.
The more comprehensive and elaborate empirical validation covered in chapter 8 demonstrates the feasibility and effectiveness of the entire approach based on a typical object (tripled glazed shaded window). However, due to a lack of precise spectral optical measurements for the glazing and blind material, some generic assumptions are required. In the validation cases that replicate lab measurements on a singular surface (section 4.9.3), empirical data with high accuracy is available. In the case of full-system measurements, the calculation results must be compared with measured data containing a significantly higher error level. This must be taken into account, as full-system measurements under varying real-world conditions are challenging. Therefore, the validation of the physical optics modules (section 4.9.3) proves that the algorithm works accurately and precisely on the material surface level, while the system validation demonstrates that the method can handle complex objects consisting of many different materials and geometrical shapes. The latter also demonstrates how results with good accuracy can be achieved, even if the required spectrally resolved information is unavailable.

Further, the comparison of in-plane irradiance computed with the *RadiCal/SIOP* approach versus classical calculations performed in a dynamic building performance simulation represents an essential validation step (see section 6.2.6.). It proves that, although the *RadiCal* approach allows the consideration of complex optical scattering processes and complex geometries, it seamlessly and precisely transitions to the classical approach if a simple object is modelled.

No intermodal comparisons with the commonly applied calculation methods or tools are included in this work. The current state-of-the-art comprises different methods that can be performed at varying levels of detail (see chapter 3). The methods have different limitations and are often based on severe simplifications. Hence, a comprehensive intermodel comparison requires additional considerations and discussion beyond the scope of this work. Nonetheless, intermodel comparisons are considered essential and will be covered in subsequent publications.

Beyond these quantitative validations, the generation of rendered images using the *RadiCal* algorithms for backward raytracing is an additional conclusive validation step. Renderings are a powerful method to identify model shortcomings, potential singularities or numerical instabilities in the equations used. They do not allow precise quantitative analysis. However, the eye can identify model shortcomings, as it is very sensitive to singular errors that might occur. The rendered image is only able to provide a flawless, realistic appearance if billions of light scattering events - involving a multitude of configurations regarding surface orientations, wavelengths and light paths - are computed correctly.

## 2.5. Separation of the optical and thermal domain

Any attempt to provide precomputed static directional solar heat gain coefficients (DSHGC) to calculate the total solar energy transmitted through a building component is linked to systematic errors. The total energy induced by solar radiation comprises optically transmitted energy and thermally transmitted secondary components. Deriving DSHGC values, therefore, implies solving a thermal balance equation with predefined thermal boundary conditions and thermal material properties. This means that in order to derive the commonly applied DSHGC values, it is necessary to anticipate the results and boundary conditions of the dynamic thermal simulations to which they will be applied. This actually constitutes a circular reference.

The DSHGC approach emerged from gradually adapting the constant total transmittance coefficient approach (g-value). It does, however, not seem reasonable to try to adapt this concept to assess the dynamic performance of glazings. The total transmittance coefficient was defined to provide a method
for simple performance characterisation of windows and glazing systems and relies on standardized and stationary boundary conditions. The method is not designed to accurately model dynamic effects related to thermal mass in the glazing or façade under real-world irradiance. While some errors occurring on dynamic simulations will cancel out when integrated heating or cooling energy demands are analysed, others will not, as they will effectively change the heating and cooling load profiles. Beyond that, dynamic effects significantly affect interior surface temperature profiles, which is relevant for thermal comfort assessment. In order to overcome these issues, the novel approach developed here follows a strict separation into solar radiation modelling (optical domain) and thermal modelling (thermal domain) (see Figure 4). The results provided by the *RadiCal* method allow for accurate calculation of all relevant absorbed and transmitted solar radiation components during the dynamic thermal simulation process, overcoming the shortcomings of the DSHGC approach and making its use obsolete.

<figure><img src="image_9.png" /><figcaption>Figure 4 Comparison of the workflow for modelling a shaded triple-glazed window DSHGC method (left) vs. RadiCal method (right) (red dotted: thermal domain, grey-dotted: optical domain)</figcaption></figure>

## 2.6. Implementation

As the objective of the PhD is to provide an accurate and efficient calculation method based on raytracing and physical optics, the work performed in the course of the thesis reaches beyond the theoretical formulation of the novel approach. The method's feasibility can only be shown based on its actual implementation and application to practical use cases. In total, eighty-two individual modules comprising more than 46.000 lines of Pascal source code were written for the PhD. This does, of course, comprise much code that is not relevant for the final, functional code of the *RadiCal* method, such as testing and validation routines as well as obsolete implementations that were replaced with alternative approaches. Hence, regarding the work performed, the finally selected and presented models and code snippets can only be considered the tip of the iceberg. Although applying measures to software development productivity is challenging, an accepted standard publication
(Jones, 2017) states a value of 615 lines of code per month as an industry average for pascal programming, which would total 75 months of programming for the code written during the thesis. The requirements for documentation in the software industry are higher, of course, and hardly comparable; however, the code implemented here was of a comparably complex nature, and apart from the actual coding, many test runs, validation cases and performance optimisation had to be performed. Due to the large extent of the written source code, presenting or discussing the entire implementation is impossible. Nonetheless, code snippets covering the implementation are essential elements of the thesis and are included in most sections. Their purpose is to provide insights regarding the implementation's structure and highlight essential algorithms or design choices. The validity of the implementation is proven based on carefully selected tests (see above).

## 2.7. Programming language

Most of this work was implemented in the programming language *Pascal* with the free compiler *FPC3.2*, embedded in the IDE (integrated development environment) *Lazarus 2.2*. The IDE allows cross-platform compilation for different operating systems and can be used freely under the *GNU GPL* license. Although the use of Pascal has been declining in recent years, several reasons led to the decision for this programming language.

The most significant reason is performance. Since Pascal code is compiled into the machine language of the desired target platform, the code is very performant. As the *RadiCal* method relies on a stochastic Monte-Carlo approach that involves a large number of iterations, it is essential to strive for the lowest possible execution time. Parts of the implemented Monte-Carlo algorithm were implemented in Python and found to be running 50-70 times slower than in Pascal. This performance gap would make any Python implementation virtually impossible, as the currently multi-hour-long ray-tracing scan would take several days to weeks, making the method impractical.

Another requirement for the choice of the programming language is its support of object-oriented programming (OOP). All relevant features of OOP, such as *polymorphism, inheritance, abstraction and encapsulation*, are available in Pascal and are crucial elements of the *RadiCal* implementation. In addition to these two mandatory requirements, further factors supported the choice. Pascal is considered to have good syntax readability in terms of clarity and consistency. Further, the operator overloading functionality allows straightforward implementation of the required algebra for complex-valued numbers as well as for three- and four-dimensional vectors and matrices. Finally, the author has been programming in Pascal-based languages (*Turbo Pascal, Delphi, Lazarus*) for several decades, significantly increasing the implementation's efficiency.

Beyond using Pascal for the main modules, C++, PHP, JavaScript, HTML and SQL were used to implement the control, monitoring and evaluation software for the measurement process covered in chapter 8.
# 3. State-of-the-art analysis

The following sections cover a state-of-the-art analysis of modelling methods for the energetic performance of shaded or unshaded fenestration systems. The *RadiCal* method can also be applied to other applications that rely on the accurate calculation of solar-induced energy flows, such as all types of solar thermal collectors, photovoltaic modules or building integrated applications of these. However, a comprehensive description of the currently applied modelling methods for these applications would go beyond the scope of this work. Therefore, the following analysis focuses on the presently applied calculation methods for fenestration systems, as the *RadiCal* method was initially developed for this purpose.

Due to its practical relevance, scientific research regarding the energy performance of windows and glazings is closely related to software tools and international standards. Consequently, this chapter contains many references to these tools and standards. Any research that involves the calculation of the thermal performance of windows has to rely on software, as the heat-balance equation has to be solved considering the glazing's specific spectral and thermal properties.

Slightly simplifying, it can be stated that the essential results of the research performed before the year 2000 are laid down in the relevant international standards. These methods have in common that directional properties of the incidence light are generally neglected, as the direction of light is mostly assumed to be of normal incidence.

As a result of research addressing the directional limitations, more sophisticated approaches have evolved after the year 2000. These methods are referred to as “modern approaches” in this chapter. The directional methods were gradually integrated into relevant tools, in particular into the commonly used *Window* software by Lawrence Berkeley National Laboratory (LBNL Windows and Daylighting Group, 2022b) and into *EnergyPlus* (U.S. Dept. of Energy, 2022) starting in the year 2011 (Lee and U.S. Dept. of Energy, 2017).

More recent developments in the research field, like the *Comfen* tool by LBNL (Mitchel *et al.*, 2019) or the *Fener* tool (Bueno *et al.*, 2015) concentrate on dynamic simulations to assess energy demand or daylighting aspects, allowing users to consider different devices and control strategies. However, all of these new methods are developed on the existing optical models and do not question or extend them.

## 3.1. General model and definitions

<figure><img src="image_10.png" /><figcaption>Figure 5 Components relevant for the total solar transmittance</figcaption></figure>

According to the applicable standards (see section 3.4), the total solar energy transmittance (or solar factor) $g_{tot}$ or (solar heat gain coefficient SHGC) describes the ratio of the total heat gains to the incident solar radiation. It can be written as the sum of two components:
$$g_{tot} = \tau_{e} + q_{i} \qquad (1)$$

The first component $\tau_e$, referred to as *solar direct transmittance*, reflects the fraction of solar radiation which is directly, i.e. "optically", transmitted into the interior of a building in the form of shortwave radiation.

The second component $q_i$, referred to as *secondary heat transfer factor* represents the fraction of solar radiation which is absorbed by the individual layers of the fenestration system, transmitted to the innermost surface of the system and subsequently released to the interior by convection and thermal radiation. Just like $\tau_e$, the quantity $q_i$ is dimensionless and not a unit of energy, as indicated by using the letter q.

While not explicitly stated in the formulation presented below, the relevant international standards (ISO 15099, 2003; ISO 9050, 2003; EN 410, 2011) implicitly define the secondary heat transfer factor $q_i$ as:

$$q_i = \frac{\Phi_{q,int}(E_S) - \Phi_{q,int}(0)}{E_S} \qquad (2)$$

Where $E_S$ is the incident solar radiation flux density, i.e. irradiance, while $\Phi_{q,int}(E_S)$ and $\Phi_{q,int}(0)$ represent the heat flux densities to the interior with and without solar irradiance. The considered heat fluxes comprise components of thermal radiation, convection and conduction. While there are some differences regarding the components' definitions and calculation methods, the applicable standards define secondary heat transfer as the amount of additional heat flux induced by solar irradiance. This also comprises non-linear effects, such as the change of thermal conductivity of gas layers due to an increase in absolute temperature and temperature gradients.

The relevant standards (discussed in section 3.4) provide detailed information regarding the thermal modelling of the multi-layer system (e.g. consideration of ventilation, buoyancy, change of gas viscosity etc.). However, the major determining factor for all secondary effects is the absorption of solar radiation in the individual layers of the fenestration system. Reliable thermal performance values can, therefore, only be obtained based on a model that is able to provide accurate layerwise absorption and transmission values for any given incident solar radiation profile.

## 3.2. Modern approaches - BSDF and DSHGC

### 3.2.1. The Klems approach

The modern approaches mainly derived in the recent two decades are based on concepts that date back to earlier work. In his work of 1994 (Klems, 1994b, 1994a), Joe Klems presented a model that allows to describing the solar-optical properties of complex fenestration systems by characterising the layers via matrices that contain spatially averaged but angularly resolved bidirectional reflectance and transmittance values. The work of Klems is based on earlier publications of Lawrence Berkely National Laboratory (LBNL), in particular on a paper in 1986 (Papamichael et al., 1986) introducing the matrix concept and a 1988 paper (Papamichael et al., 1988) that proposes a scanning radiometer measurement method to acquire the necessary optical information for the method.

In the Klems matrix approach, the actual spatial (i.e. locational) information is averaged; however, the directional information is kept in an angularly resolved way. For this purpose, Klems introduced an angular segmentation of the hemisphere into 145 patches covering approximately equal solid angles. In order to describe the directional scattering properties of a layer comprehensively, the transmittance or reflectance values for the 145 outgoing directions have to be provided for all of the
145 potential incident directions (see Figure 6). Therefore, the entire set comprises 145 x 145 = 21,025 elements. The elements forming the matrices are, of course, not independent, as the fundamental laws of the underlying optical principles have to be obeyed, particularly: positivity, Helmholtz reciprocity and energy conservation. The Klems matrix is, therefore, a discrete representation of a more general concept known as *bidirectional reflectance/transmittance distribution functions (BRDF/BTDF)* or their combination, referred to as *bidirectional scattering distribution functions (BSDF)*. The BSDF concept and terminology were introduced even earlier, e.g. by Nicodemus (Nicodemus, 1965) for analysis of longwave radiation and more comprehensively in 1980 by Bartell (Bartell, Dereniak and Wolfe, 1981). Klems mainly avoids this terminology but refers briefly to Nicodemus in his first work.

<figure><img src="image_11.png" /><figcaption>Figure 6 Klems matrix representation of LBNL BSDF Viewer.
*(Figure taken of radiance-online.org (McNeil, 2022))*</figcaption></figure>

However, in his two important publications in 1994, Klems derives a methodology for calculating solar-optical properties of multi-layer fenestration systems from their individual bi-directional layer properties. Klems begins with non-specular layers and gradually develops the model further to finally include specular layers and consider interreflection between the layers.

Klems notes that his approach is based on two significant approximations:

1. *Spectrally averaged optical properties:* The input and output reflectance and transmittance information are valid for integrated wavelength ranges, usually either the visible or global radiation spectrum.

2. *Spatial averaging:* While the angular information is tracked, the actual spatial information is lost, as an averaging over all points of a layer is assumed in the model.

Beyond these two approximations raised by Klems, the author of this work believes a third one should be added, although it is related to the second point:

3. *Infinite extension:* The layer model intrinsically assumes that all layers are infinitely extended, and any radiation transmitted by one layer will reach the next layer. Therefore, the model cannot cover any boundary effects.
In his work, Klems expressed most concerns regarding the first (spectral) issue and its practical limitations. In contrast, he argues that the second (spatial) issue has limited practical consequences, and he does not raise the third issue. The following section discusses some implications regarding the simplifications contained in the Klems method.

The angular transmittance (BTDF) and reflectance (BRDF) values necessary to define the Klems matrices of a fenestration component can either be determined by measurement or by modelling. In the empirical approach, measurements are mainly performed using spectrophotometers or goniophotometers (Kohler, 2014). While goniophotometers are frequently used to determine optical properties for daylighting applications, their use to measure integrated values for the entire global radiation spectrum is still a matter of research (Curcija et al., 2022; Ward et al., 2022). Since the commonly used silicon sensors cannot capture radiation with wavelengths beyond approx. 1000 nm, a combination with InGaAs or germanium sensors has to be used to extend the measurement capabilities into the near-infrared region. The modelling approaches to determine BSDF matrices are either based on analytical optical calculations, as applied in *LBNL Window* (LBNL Windows and Daylighting Group, 2022b) or raytracing. In the most commonly used method, however, a combination of measure and modelled data is used: the optical properties of the material's surfaces are determined empirically and subsequently used in a two- or three-dimensional model to determine the optical properties of the entire system. The reason for that is that, on the one hand, most goniophotometers cannot perform measurements on largely extended objects. On the other hand, the approach allows a more flexible and efficient workflow, as the measured material data can be used to model a variety of different systems.

### 3.2.2. The black-box or two-layer model

Independent of the BSDF approach, in 2006 Kuhn (2006) and, more specifically, in 2011 Kuhn, Herkel *et al.* (2011) propose a significant simplification for the modelling of a shaded glazing system by separating the fenestration system into a shading layer and a glazing system layer:

> The proposed model for $g_{tot}$ introduces an additional simplification level and treats the whole glazing unit as one layer and the shading device as another. One advantage is that the angular dependence of the glazing and shading device can be characterized externally with the method that best meets the requirements for the system under consideration. Another advantage is that manufacturers of glazing and shading devices have to guarantee only for the properties of their product, not for the combination. (Kuhn, 2006)

Providing accurate descriptions of the angular-dependent properties of the glazing and shading systems is still necessary for this approach. Kuhn proposes the use of a semi-empirical model suggested by Roos *et al.* for the glazing layer (2001). The Roos model and underlying models (Montecchi *et al.*, 2002) are, in fact, based on physical optics models and consider Fresnel's equations and the polarisation state of light. However, finally simplified semi-empirical fitting models are derived. Regarding the angle-dependent optical properties of the shading device layer, Kuhn suggests either (whole system) measurements or the use of 2D raytracing, as demonstrated in his earlier work (Kuhn *et al.*, 2001).

Finally, a virtual two-layer model is used to calculate direct transmittance and solar absorptance in the two virtual layers to derive the secondary heat transfer. For this approach, it is necessary to define thermal boundary conditions, as the relevant thermal transfers leading up to the directional total transmittance values must be considered. For different applications and configurations (e.g. internal or external shading) several equations are provided on how the layerwise models can be
integrated to finally derive angle-dependent relations for the total transmittance $g_{tot}$, also referred to as the directional solar heat gain coefficient (DSHGC). The precomputed dataset can then be integrated as ‘black-box’ model into existing dynamic building performance simulation tools. The proposed method allows for a consideration of angular effects and therefore represents a significant improvement over the ISO standard methods. However, the newly introduced simplification made by assuming two independent layers poses a new potential source for inaccuracy. The approach does not sufficiently allow the modelling of specific spectral and angular filtering or conditioning of the solar radiation by the first layer. The angular and spectral profile of the radiation incident to the second layer will be significantly different from the "neutral" radiation profile that is assumed for the calculation of its optical properties. To mitigate some of these issues, analytic relations that mix the layer's solar and visible optical properties are presented. Furthermore, interreflections between the two layers cannot be modelled accurately, as a simplified model for diffuse radiation is assumed.

### 3.2.3. BSDF and the two-layer model

Regarding the simplified angular information contained in the Kuhn model, Bueno and others (see Bueno et al., 2015; Bueno et al., 2017) propose an improved approach, as the two-layer model of Kuhn is ‘upgraded to accept BSDF data’. The method described is based on Kuhn's two-layer model; however, the directional solar heat gain coefficients DSHGC of the glazing unit and the BSDF of the glazing unit and the shading device are considered. The directional information for reflectance, transmittance, absorptance and total transmittance is used in the matrix form proposed by Klems. Regarding modelling the resulting energy flows, the ‘upgraded model’ described by Bueno et al. (2017) still applies the principles as suggested by Kuhn for the two-layer model. In turn, Kuhn's initial model is based on the simplified approach provided in the standard ISO 52022-1. In the equation proposed for modelling external Venetian blinds, three components are considered:

$$g_{sys} = \tau_{sd} g_{gu} + \alpha_{sd} \frac{\Lambda_{ext}}{\Lambda_2} + \tau_{sd} (1 - g_{gu}) \frac{\Lambda_{ext}}{\Lambda_1} \quad (3)$$

$$\text{with } \Lambda_{ext} = \left[ 1/U_{gu} + 1/\Lambda_1 + 1/\Lambda_2 \right]^{-1} \quad (4)$$

<table><tr><td>g<sub>sys/gu</sub> ... total solar transmittance of</td><td>τ<sub>sd</sub> ... transmittance of shading device</td></tr><tr><td>system/glazing unit</td><td>α<sub>sd</sub> ... absorptance of shading device</td></tr><tr><td>Λ<sub>ext</sub>·Λ<sub>1</sub>, Λ<sub>2</sub> ... thermal conductivities [W/(m<sup>2</sup>K)]</td><td>U<sub>gu</sub> ... thermal transmittance of glazing</td></tr></table>

The optical coefficients ($\tau, \alpha$) as well as the $g$ values in equation (3) are functions of the incident light direction; for reasons of clarity, the corresponding arguments were omitted in the above equations. The corresponding directional values are taken from the Klems matrices representing the BSDFs of the shading and glazing unit as well as the DSHGC values of the glazing unit.
The first term of equation (3) simply scales the *g* value of the unshaded glazing unit by the transmittance of the shading. The two subsequent terms represent additional heat gains due to the presence of the shading device. While not explicitly stated in the ISO standard or Kuhn's or Bueno's publications, it is apparent that the additional two terms result from a simple heat transfer model comprising three resistors (conductors) in series, see Figure 7. The middle term of equation (3) represents the fraction of heat absorbed in the shading device and transmitted to the interior, while the last term represents an additional fractional amount of heat trapped in the glazing unit due to the presence of the shading device. The equation proposed for internal blinds follows a very similar approach but contains additional terms to consider reflections.

<figure><img src="image_12.png" /><figcaption>Figure 7 Thermal model for the two additional heat gain components</figcaption></figure>

### 3.2.4. LBNL Window and Radiance

In the currently most advanced and probably most commonly used approach for modelling the energetic performance of complex fenestration systems, the free tools *LBNL Window* (LBNL Windows and Daylighting Group, 2022b) and *Radiance* (McNeil *et al.*, 2013) are used simultaneously. *LBNL Window* allows a quite accurate modelling of glazing systems comprising coated and/or uncoated glass panes. For normal or near-normal incidence, the results of the calculations can be considered very accurate, as *LBNL Window* uses the measured product data contained in the extensive *IGDB - International Glazing Data Base* (LBNL Windows and Daylighting Group, 2021) maintained by the same workgroup.

However, regarding the important directional dependencies of the optical coefficients for the glass panes, *LBNL Window* makes use of simplified models that consider unpolarised radiation only. There are currently two different angular models implemented in the *Window* software (Curcija *et al.*, 2018). The approach for uncoated glass utilizes the empirical transmittance and reflectance spectra, measured at normal incidence only, to estimate a spectrally resolved refractive index in the first step. Subsequently, this quantity is used in Fresnel's equations to derive the angular reflectance and transmittance of the glass. For coated glass panes, the angular properties are determined based on a regression fit. For the fitting process, two sets of fitting parameters are provided. A threshold value of 0.645 for the transmittance is used to determine if the coefficients for clear glass or bronze glass are used in the polynomial regression formula. The same approach for calculating angular properties is applied in the building performance simulation tool EnergyPlus (U.S. Dept. of Energy, 2022).

Regarding spectral resolution, *LBNL Window* allows various approaches. The spectral resolution of the calculation is mainly determined by the specular information available for the components. Systems containing exclusively glass layers usually allow spectrally resolved computation. If a shading layer is present, the calculations are primarily performed for integrated spectral ranges (e.g. the solar or visible range).

Shading devices are integrated into the calculation process in the form of shading layers. *LBNL Window* allows the following options (Mitchel *et al.*, 2019):

* *Homogeneous diffusing shade*

* *Perforated screen*

* *Woven shade*
* Venetian blind, horizontal

* Venetian blind, vertical

* Shade with XML data

* THERM file

The software performs calculations based on analytical or geometrical methods for the first five

<figure><img src="image_13.png" /><figcaption>Figure 8 Workflow for modelling with LBNL Window and Radiance/genBSFD
(figure taken of WINDOW Technical Documentation (Curcija, Vidanovic, et al., 2018))</figcaption></figure>

options. The materials used for the shades are generally considered to be diffusely reflecting (Lam-
bertian surfaces). For the simulation of Venetian blinds, *LBNL Window* has implemented a method
that complies with the standard ISO 15099. However, a set of options and additional features are
available that are beyond the scope of the standard and should be able to provide higher accuracy.
For example, the geometry of the slats can be assumed as non-flat by providing a curvature radius,
or the simplified modelling of diffuse radiation can be overruled by selecting a directional diffuse
calculation method.

In order to model more complex shading devices or to model shading devices in a more detailed way, *Window* allows the integration of externally computed BSDF data. The options *Shade with XML data* and *THERM file* are available for this purpose. The latter option allows the definition of the geometry and material using LBNL's *THERM* software (LBNL Windows and Daylighting Group, 2022a).

In this approach, the actual BSDF information is not generated by *LBNL Window*, but with the
external tool *genBSDF* (McNeil, 2015). This program is a three-dimensional raytracing representing
a spin-off of the powerful open-source daylighting tool *Radiance* (Radiance Community, 2022). An
overview of the workflow is depicted in Figure 8. The illustration is taken from the technical docu-
mentation for *LBNL Window* (Curcija et al., 2018).
The *genBSDF* tool offers no graphical user interface. It is a command line tool; the necessary input and output are provided as scripts and reports.

Since *Radiance* is a rendering software, it follows an RGB approach to achieve realistic visual representations. To still use the raytracer for energetic purposes, it is commonly recommended to apply *Radiance* material definitions of type *plastic* or *translucent plastic* and use identical values for the red, green and blue reflectance (or transmittance) values. The constant reflectance values provided should represent the measured spectrally integrated reflectance for either the solar or the visible band. The resulting BSDF files are then valid for the provided spectral range but do not contain any spectrally resolved information. This can be problematic for modelling systems that contain spectrally selective layers, as will be pointed out in the next section.
## 3.3. Limitations of BSDF and DSHGC approaches

### 3.3.1. Spectral accuracy

The most severe limitation that BSDF-based approaches share is probably the spectral simplification on which the methods are based. The values forming the Klems matrix are based on integrated spectral properties considering different bands, such as the global radiation spectrum or the visible light spectrum. In his work of 1994, Klems (1994b) himself already pointed out that the BSDF approaches focus on the ‘spatial and angular aspects of systems rather than the spectral properties’. He even showed in the paper’s appendix that his approach only works if the spectral selectivity is restricted to one layer at most. All other system layers have to be considered *spectrally flat*, i.e. having the same reflectance or transmittance value over the entire considered spectral range. Klems notes explicitly:

> As the method is applied to progressively more complicated systems, such as colored blinds or shades combined with tinted or selective glazing, or multiple blind or curtain systems of different colors, the treatment could become less accurate and eventually inapplicable. (Klems, 1994b, p.4)

Excluding rare special photonic applications that are able to perform wavelength conversion, generally, the wavelength of the incidence radiation will not change while propagating through the fenestration system, meaning that the initial spectral components will be absorbed, transmitted or reflected but not absorbed and reemitted with a different wavelength. Hence, the Klems approach can be multiplied into a series of independent, parallel calculations limited to appropriate spectral bands. However, this increases the complexity of the method significantly, and it is difficult or almost impossible to define an "all-purpose" segmentation of the global radiation spectrum that would be suitable for most applications. Internally, *LBNL Window* can use 15 different bands (17 wavelength points) to cover the global radiation spectrum. However, the approach is not fully implemented as the BSDF information of the shading device is not available in this format. The raytracing tool *Radiance* used to compute the BSDFs of shading uses only integrated values. Therefore, depending on the input data, either only the solar (=entire global radiation) or only the visible spectrum is usually considered in an integrated way. This separation into visible and solar radiation is also frequently found in standard building performance simulation tools. However, for energetic purposes, the use of such integrated values has to be questioned. Considering the entire radiation spectrum from ultraviolet to near-infrared, the assumption of 'flat' reflectance is not suitable, as even objects that appear white, grey, black or clear (transparent) in the visible light spectrum will usually be spectrally selective in the non-visible ranges. Hence, using the BSDF information of independent components to determine the properties of the entire system may imply a significant systematic error.

### 3.3.2. Spatial averaging

Generally, all currently applied methods are based on the Klems approach that captures angularly resolved information on the reflected or transmitted radiation but does not include any information regarding its spatial (i.e. locational) distribution. Therefore, the approach is intrinsically based on the assumption that each layer represented by BSDF information is homogeneous and infinitely extended. Consequently, the impact on any other layer caused by an inhomogeneously structured layer or geometric boundary effects cannot be modelled with this approach. Boundary effects are
generally relevant for virtually every glazing, as the glazing is either embedded in a window reveal or interrupted by framing structures. The significance of the boundary effect will, of course, decrease when the undisturbed area is large relative to its boundary.

The effect of the inhomogeneity of a layer is relevant in combination with boundary effects or any further inhomogeneous layer. Klems notes that the effect can be neglected as, typically, there will be no relevant spatial correlation between any two layers. This is true for small-scale structures, like fabric, grids or, to some extent, Venetian blinds. However, this will generally not be the case for devices with large-scale structures in the dimension of the glazing. Typical examples are horizontal louvres or overhangs mounted on top of a window, but also glazings or films with fritted or printed structures. More generally, it can be said that the effect is significant for all systems where the device's actual positioning relative to the glazing or any layer with structures of the same dimensions is relevant. Very apparent examples are provided by kinetic shading mechanisms where fritted or coated structures can be moved relative to each other to change the system's solar transmittance, e.g. of modern foil constructions (Flor et al., 2018). While the layerwise BSDF-based modelling would entirely break down for such applications, it can also lead to significant errors for more classical applications comprising large structures. This is the case for many vertically and horizontally oriented systems, as well as for fixed and movable devices, such as louvres, overhangs, hinged shading systems, canopies, fins, brise-soleis etc. In most of these applications, the spatial concentration of light shading or diffusing effects to specific regions is not an unintended or marginal side-effect but the actual intended effect. The geometric properties are usually specifically designed to exploit such behaviour. Therefore, it is essential that any methods used in the design or performance assessment process for such systems can capture such effects accurately.

### 3.3.3. Polarisation

At present, all commonly used methods for modelling fenestration systems that include shading devices do not consider the change of the polarisation state of the transmitted or reflected light beam. On the one hand, the ISO standard methods (see section 3.4) exclude polarisation by definition, as they are restricted to near-normal incidence. In the case of near-normal incidence, polarisation effects are irrelevant. On the other hand, the BSDF methods, which were specifically developed to model transmittance and reflectance at oblique angles, implicitly assume that the provided coefficients are valid for the intensity (or flux density) of light beams with a neutral polarisation state. This neutral state can only be assumed (to some extent) for the radiation incident on the first layer. After the first scattering process, regardless of whether transmission or reflection, any ray with an oblique incidence angle will show a significant polarisation that will affect the outcome of any further scattering process (see sections 4.2.4 and 4.7).

In principle, it would be possible to include polarisation information in the BSDF approach using dual calculation paths for perpendicular polarisation states; however, this has not been proposed in the current publications. The reason for this might be a lack of calculation tools that are able to create the BSDF matrices for polarised radiation. At least for unshaded glazings, *LBNL Window* might be adapted quite easily for such purpose.

### 3.3.4. Angular dependence

As discussed, the angular dependence of the transmittance and reflectance are derived based on simplified models. In fact, the angular dependencies, in turn, depend on the wavelength of the incidence light. While the effect is limited in the visible spectral range for bulk materials, there is still a
significant variation if the entire global radiation spectrum ranging from 300 nm to 2500 nm is considered. Furthermore, the effect is particularly pronounced if coatings are present, as spectral filtering is the designed function of these layers. Apart from these intended effects, exploited, e.g. for anti-reflective coatings, solar-protection coatings or low-E coatings, there can also be unintended effects that can, for example, lead to a colouring effect if coated materials are viewed at higher angles. All mentioned angular effects cannot be modelled appropriately based on simplified models. Valid results can only be achieved by solving Fresnel equations in a spectrally resolved way and by using the materials' accurate spectral properties.

### 3.3.5. Modelling of direct solar radiation

Another issue that is linked to the usually limited angular resolution of the Klems-matrix is the modelling of direct solar radiation. In the commonly applied 3-phase method, the contribution of direct radiation is included in the angular patch corresponding to the sun's incidence direction. This is also true for the outgoing Klems patch if the system does not diffuse the direct beam. Hence, on a clear day, the sun's radiance value is reduced by several orders of magnitude and instead distributed on a Klems patch covering a relatively large solid angle (see Figure 6 and Figure 75). As this averaging does not represent the true angular extent and radiance of direct solar radiation, it will cause significant inaccuracies if further BSDF layers are considered. This is especially true for the inhomogeneous structures discussed in the *Spatial averaging* section above.

### 3.3.6. Energetic model

The underlying thermal model of the two-layer method (Bueno *et al.*, 2017) and previous DSHGC approaches are based on necessary simplifications. The thermal model does not consider the layered structure of the glazing itself, which appears essential as modern tripled-glazed units show high thermal resistance values that will cause significant temperature differences among the individual glass panes. On the one hand, the temperature deltas will directly affect the energy balance model applied in the approach. On the other hand, it will have a secondary effect, as the absolute and relative temperature changes will significantly impact the gas layers' thermal conductivities. Beyond this issue, the choice of the constant thermal conductivity parameters $G_i$ that determine the additional heat flow from the two layers to the interior seems somewhat arbitrary. The recommended values are likely to have empirical origins, but no further explanations, options or models for adaption to different fenestration systems or climates are provided. Apart from these issues, assuming stationary thermal properties by applying constant thermal transmittance (U-value) to characterise the glazing system seems problematic. Though some fitting formulas to adjust the U-value for different boundary conditions (climates) are presented, these relations will usually strongly depend on the actual layering of the glazing systems and actual thermal boundary conditions. Hence, it is difficult to propose a generalised approach. All methods providing DSHGC values have to anticipate results and boundary conditions of the dynamic thermal simulations to which they will be applied. This constitutes a circular reference that can only be addressed with the assumption of constant standard values.

The integration of BSDFs created with *LBNL Window* (or otherwise) into *EnergyPlus*, referred to as 3-phase method, allows layer-wise energetic modelling similar to the method proposed here. However, the approach suffers from the BSDF-related issues discussed above, especially for fenestration systems containing shading layers.
## 3.4. Approaches of international standards

### 3.4.1. ISO 9050 and EN 410 for glazing systems

The standards ISO 9050 (2003) and EN 410 (2011) provide almost identical information on how to calculate the optical and energetic properties of glazings. In particular, they both provide equations and boundary conditions for the calculation of light transmittance and reflectance coefficients, as well as for the total solar energy transmittance $g_{tot}$. Apart from the calculation instructions, both standards include a set of tables that provides the relevant spectral distributions necessary to model the spectral power distribution of the global radiation spectrum and visible light. EN 410 additionally contains spectral distributions relevant to colour rendering, colour perception and a standard illuminant light source. ISO 9050 does not explicitly contain this information but refers to the applicable CIE standards for that purpose. Instead, the supplemental spectral distributions provided in ISO 9050 can be used to evaluate the damaging properties of light for artwork or general purposes as well as for human skin.

To determine the optical and energetic behaviour of a glazing system consisting of one or multiple layers, the shortwave radiation and heat balance of the system is modelled and solved. The shortwave radiation balance considers reflection, absorption and transmission in all relevant components. Consequently, interreflections between the individual interfaces are also considered in the model. All calculations are performed taking into account the relevant spectral distributions. For calculation of the total solar energy transmittance, the global radiation spectrum ranging from 300 nm to 2500 nm is considered, whereas, for the visual light properties, the appropriate spectral distribution between 380 nm and 780 nm is applied. The detailed consideration of the spectral distributions is essential for the accurate modelling of solar radiation.

However, concerning the angular profile of the incident solar radiation, both standards show a severe limitation as the scope of the model provided is restricted to 'quasi-parallel, near normal radiation incidence' (ISO 9050, 2003, p.2; EN 410, 2011, p.8). This poses a significant restriction as it reflects an idealised case that does provide benefits regarding calculation, spectral measurement or lab testing but has little relevance for the modelling of real-world applications.

Based on the proposed shortwave radiation model for normal incidence, a heat balance model is solved to determine the secondary internal heat gains caused by solar radiation. While the underlying equation systems can be solved for any number of layers, only the solutions for the cases of single-glazing, double-glazing and triple-glazing are provided in an explicit analytic form in both standards.

### 3.4.2. ISO 15099 approach for fenestration systems with a shading layer

The international standard ISO 15099 (2003) provides a methodology to calculate the thermal performance of fenestration systems. In order to consider the entire heat balance, the standard provides models for thermal transmission (i.e. conduction, convection and thermal radiation exchange), as well as for energy flows caused by solar irradiation. The models are based on the assumption of a layered quasi-parallel structure covering glass panes, air or gas gaps and shading layers. The thermal performance of the glazing is determined by setting up and solving the energy balance of the layered system. The standard provides many details regarding suitable calculation methods for modelling heat exchange via radiation (longwave), conduction and convection. However, regarding the solar (shortwave) irradiance components, it refers to ISO standard 9050 (2003) or applies its methods. As mentioned in the previous section, ISO 9050 and EN 410 do consider the spectral distribution of the
incident radiation but are restricted to the case of near-normal incidence, which poses a significant restriction for accurate calculations relevant to real-world applications.

While the ISO standard 15099 does not provide additional information regarding the shortwave model of the glazed layers, it does contain specific information regarding the modelling of shading devices. The model is, however, restricted to shading devices 'which are or which may, by proper approximation, be treated as a layer parallel to the pane(s) of the window' (ISO 15099, 2003, p.31).

<figure><img src="image_14.png" /><figcaption>Figure 9 Separation into disturbed and undisturbed components according to ISO 15099</figcaption></figure>

In order to integrate such a shading layer into the fenestration system, a separation into a *disturbed* and *undisturbed* component of the incident radiation is introduced. The disturbed part represents any radiation diffusely reflected or transmitted by the surfaces of a shading device (see Figure 9). This component is assumed to be perfectly diffusing (Lambertian). According to the model, both components, the direct (undisturbed) as well as the diffuse (disturbed), propagate through the fenestration systems independently. It is, therefore, necessary to provide diffuse transmittance and reflectance values for all layers contained in the system. The standard notes this fact; however, it does not state how these values should be calculated for glass panes or coatings, as the underlying standard ISO 9050, which it refers to, is restricted to 'quasi-parallel, almost normal radiation' (ISO 9050, 2003).

The standard ISO 15099 does, however, provide a detailed calculation method for the necessary coefficients for the case of idealised Venetian blinds. The model assumes an infinite array of tilted, parallel, infinitely extended, perfectly flat and ideally diffuse reflecting slats. Based on these assumptions, it is possible to derive the necessary reflectance, transmittance and absorption values for the shading layer for the three components *direct*→*direct*, *direct*→*diffuse* and *diffuse*→*diffuse*.

The view-factor-based calculation relies on a model defined by the input parameters slat angle $\gamma$, front-and back-side slat reflectance ($\rho_f(\lambda)$, $\rho_b(\lambda)$), slat width and slat distance.

By introducing a shading layer, characterised by these reflectance, transmittance and absorption coefficients, into the energy balance of the entire fenestration system, it is possible to approximate the energetic effect of the shading layer.

Apart from the substantial simplifications necessary for the idealized Venetian blind model, several other significant issues have to be considered in the approach:

* unspecified determination of glass or film properties for diffuse radiation component

* assumption of perfectly diffused radiation caused by the shading layer

* assumption of normal incidence for the direct radiation component for the involved glass layers

Despite these issues and the potentially significant errors associated with them, the ISO 15099 approach is widely used, even in modern building performance simulation tools.
### 3.4.3. Standard approach ISO 52022 series for fenestration systems with a shading layer

While ISO 15099 is widely used in the US, most European countries refer to the similar ISO 52022 (2018) standard series. The first part of the three-part series contains a method to estimate total solar transmittance, direct energy transmittance and light transmittance based on a simplified approach for various configurations. This simplified approach will not be discussed here.

Similar to ISO 15099, the third part of the standard series contains a section which specifies how to implement a shading layer into the energy balance of a fenestration system. Instead of referring only to ISO 9050 for modelling the shortwave radiation exchange, the standard primarily refers to the EN 410 (2011) standard. The information contained in EN 410 is very similar to ISO 9050, and, in particular, the severe restriction to 'quasi-parallel, near-normal radiation incidence' is explicitly included again.

Regarding the reflective and transmissive properties of shading layers, the standard refers to EN 14500 (2018). This standard covers measurement methods for 'curtain elements of solar protection devices'. The definition of curtain elements includes metallic or glass slats as well as fabrics. The standard is mainly intended for measuring the optical properties of the materials that form the device but can also be applied to the 'whole product' if the 'equipment is sufficiently large' (EN 14500, 2018, p.12). The testing methods distinguish between diffuse and direct reflectance and transmittance; however, for calculating the product properties, several references to methods covered by the ISO 52022 series or unspecified, 'more accurate methods' are included. The proposed calculation method for Venetian blinds, as detailed in ISO 52022-3, however, uses the provided information for the diffuse and direct optical properties in a very limited manner.

In analogy to ISO 15099, Annex D of ISO 52022-3 provides a geometrical calculation model to calculate the transmittance, reflectance and absorption properties of a shading layer formed by Venetian blinds. Again, the slats are assumed to be ideally flat, infinitely extended and to show no specular reflection (i.e. perfectly diffuse reflecting Lambertian surfaces).

The first significant deviation from the approach provided in ISO 15099 is that ISO 52022 restricts the calculation of shading layers to cases where no direct transmittance of solar radiation is possible. Consequently, only the *direct*→*diffuse* transmittance and reflectance are provided in section D.3 *Direct radiation* (ISO 52022 part 3, 2018, p.24) and the *diffuse*→*diffuse* transmittance and reflectance formulae are provided in section D.4 *Diffuse radiation*, whereas the *direct*→*direct* component is excluded from the calculation method. The calculation is again based on the diffuse reflectance values for the front- ($\rho$) and back-side ($\rho'$) of the blinds. The calculation method for the determination of the view-factors is not specified. Only view-factors for the specific case of slats adjusted perpendicular to the solar beam with a spacing ratio of slat distance to slat distance of $d/l = 1.0$ are provided in tabular form (ISO 52022 part 3, 2018, p.25).

Somewhat confusing, part two of the ISO 52022 standard series contains an 'informative' annex G (ISO 52022 part 2, 2018, p.23), that does consider *direct*→*direct* radiation and provides calculation adaption methods for pitched and vertical elements and additional view-factor tables. However, the Annex is not referred to in the standard and is denoted as 'informative' (instead of normative). While the model is based on these additional equations it extends the applicability of the method in a reasonable way, it is still based on the same model as presented in the actual standard method and, therefore, shares most of its limitations.

The calculation of view factors for any configurations, as it is implied in the standard, is a task that can be solved rather easily, best by using simple 2D raytracing algorithms. The significant issues of
the ISO 52022-3 method are, however, again related to the general assumption of normal incidence and the consideration of diffuse radiation. Nevertheless, there is a different approach regarding handling the diffuse component, representing the second significant difference between the two standards. Unlike ISO 15099, ISO 52022 does not suggest modelling the propagation of separate direct and diffuse components through the entire fenestration system. The additional diffuse component is solely considered for the radiation incident to the shading layer. Two individual calculation models are presented to model the effect of the diffuse and direct incidence radiation. A share of 85% direct and 15% diffuse radiation or other values 'suitable for the local climatic conditions' (ISO 52022 part 3, 2018, p.25) are then recommended as a ratio to calculate weighted averages of the two results. The finally generated "mixed" transmittance and reflectance values of the shading layer are then used in the radiation exchange model for direct normal incidence, as the general *Scope* section of the standard explicitly specifies: 'Diffuse irradiation or radiation diffused by solar protection devices is treated as if it were direct' (ISO 52022 part 3, 2018, p.1).

## 3.5. Limitations of approaches in international standards

The proposed method for modelling the performance of glazings with or without shading devices is unsatisfactory. Significant deviations to real-world performance values have to be expected due to the following issues:

* global radiation propagating through the system is modelled as quasi-parallel, near-normal incidence only (ISO 9050, EN 410)

* in combination with a shading layer, ISO 15099 suggests considering an additional perfectly diffused radiation component but does not specify how glazings and coatings should be modelled accordingly.

* ISO 52022 suggests the consideration of a diffuse component for shading devices but feeds the calculation results back into the direct-normal radiation model

* ISO 15099 and ISO 52022 provide calculation methods for shading devices only for the case of Venetian blinds.

* the Venetian blind model in both standards is oversimplified: perfectly flat, ideal diffusely reflecting, infinitely extended and infinitely thin slats are assumed.

Many issues may be caused by the fact that the initial purpose of the international standards was to provide a transparent, neutral and objective performance measure for windows and shading systems. The assumed boundary conditions, therefore, rather reflect a lab-like environment that does not reflect real-world conditions, as required for building performance simulation.
# 4. Light, sun and optics - applied principles, models and methods

The currently applied methods to simulate the propagation and impact of solar radiation are mainly based on geometrical optics and strongly simplified models. Hence, often without explicitly noting it, they are based on a strong level of abstraction and approximation. While this might be reasonable and often necessary from a practical point of view, it is essential to explicitly point out the underlying assumptions and simplifications made. To address this issue - and even if this may seem unconventional - this chapter includes very general parts addressing essential principles of light and physical optics. The field of light and optics is, in fact, very complex, and many models are only valid within strictly limited boundaries. Therefore, embedding the actual implementation within a broader context seems appropriate. After a brief historical overview, significant optical phenomena and models are introduced, and their practical relevance is discussed. Subsequently, more specific models and methods are introduced. Each of these later sections already contains details regarding the actual implementation and testing of the relevant models.

## 4.1. A brief history of light

Basic laws of geometric optics were known and exploited in ancient times in the oriental as well as in the occidental parts of the world. There is evidence (Sines and Sakellarakis, 1987) that glass was used as a magnifying lens or burning glass in ancient times already. The most common theories of light assumed that light was associated with particles or corpuscles that were either emitted by objects and perceived by the eyes (intromission theory) or, reversely, cast by the eyes towards the objects (emission theory). Euclid already formulated basic laws regarding the nature of reflection and the propagation of light in straight lines, whereas Ptolemy even detected the effect of atmospheric light refraction based on accurate astronomic observations. A remarkably detailed and accurate model of light was introduced by Ibn al-Haytham, an Arabic physicist and mathematician. Al-Haytham supported the intromission theory and identified the sun as the original source of light. He analysed and described the refraction and reflection of spherical and parabolic surfaces.

<figure><img src="image_15.png" /><figcaption>Figure 10 Timeline of essential persons influencing the current model and understanding of light</figcaption></figure>

While Isaac Newton was convinced to have found evidence that light consists of minute particles, Christiaan Huygens established a new wave theory, which was able to describe observed reflection and refraction phenomena based on the assumption that light would propagate with a lower speed in dense media. Thomas Young added supporting evidence to the wave theory by conducting interference experiments using pin-holes and double-slits. The observed interference patterns caused by
the diffraction of light could finally be explained by Auguste Fresnel, who extended and refined Huygens's wave theory, even though it is unclear if he was aware of Huygens's work.

James Clerk Maxwell established a new model which coupled electric and magnetic fields. His famous theory consists of four partial differential equations which allow a sinusoidal wave-like solution. Calculating the propagation velocity of these newly discovered electromagnetic waves, he found that it was almost equal to the speed of light. Remarkably, the speed of light had already been measured at that time. This led Maxwell to the significant conclusion stated in his publication in the year 1865:

> *The agreement of the results seems to show that light and magnetism are affections of the same substance, and that light is an electromagnetic disturbance propagated through the field according to electromagnetic laws.* (Maxwell, 1865, p.499)

The discovery of the electromagnetic nature of light did not yet set an end to the famous misconception that light would require a medium to propagate, which was called the ether (or luminiferous aether). In fact, Maxwell inspired the famous Michelson and Morley experiment in the year 1881, which was designed to determine the speed of the earth relative to a stationary ether. The seemingly paradoxical null result of the experiment, showing that light emitted on the moving earth propagates into all directions with the same constant speed, left physicists baffled. It took another 24 years and the genius of Albert Einstein to resolve this supposed paradox in his *Special Theory of Relativity*. However, in the very same year, Einstein also revived the particle theory in his explanation of the photoelectric effect by introducing the concept of photons. Einstein's theoretical concept was based on the experimental results of Hertz, Lenard and others but also on the theoretical work of Max Planck. In 1900, Planck presented a theory that could predict the spectral distribution of black-body radiation by introducing discrete (quantized) energy states. The resulting new theory of quantum mechanics by Schrödinger, Heisenberg and others put an end to the particle/wave discussion, however, leaving us with a big paradox beyond our everyday perception of the world, as A. Einstein noted:

> *But what is light really? Is it a wave or a shower of photons? There seems no likelihood for forming a consistent description of the phenomena of light by a choice of only one of the two languages. It seems as though we must use sometimes the one theory and sometimes the other, while at times we may use either. We are faced with a new kind of difficulty. We have two contradictory pictures of reality; separately neither of them fully explains the phenomena of light, but together they do.* (Einstein and Infeld, 1938, p.278)

Finally, the comprehensive field theory, called quantum electrodynamics and founded by Dirac, Bethe, Feynman and others, was able to combine special relativity, quantum mechanics and electrodynamics to provide a unified model, which is able to describe most light-matter interactions in a very accurate way. Richard Feynman's book *QED-The Strange Theory of Light and Matter* (Feynman, 1985) provides an introduction to the essentials of quantum electrodynamics in an accessible language.
**4.2. Nature of light**

**4.2.1. Relevant spectral ranges and definition of “light”**

For technical applications, it is a widespread convention to use the term “light” to refer to the visible light spectrum and “solar radiation” to refer to the entire global radiation spectrum. For the sake of readability, the more general convention is followed here: the term “light” is used to refer to the entire global radiation spectrum comprising ultraviolet, visible and near-infrared radiation (NIR), while the term “visible light” is used to exclude the non-visible spectral ranges explicitly.

<figure><img src="image_16.png" /><figcaption>Figure 11 Overview of the global radiation spectrum and relevant spectral ranges</figcaption></figure>

Following another commonly used convention, the global radiation spectrum is limited to the spectral range from 300 nm to 2500 nm. While extraterrestrial radiation still comprises approx. 3% of the total power above the threshold of 2500 nm, most of it is absorbed in the atmosphere by water vapour and CO2 molecules, and, therefore, is not relevant. Compliant with the applicable standards (see e.g. EN 410, 2011), *visible light* is defined in the present work as radiation with a wavelength ranging from 380 nm to 780 nm.

Figure 11 gives an overview of the lower half of the global radiation spectrum and the relevant wavelength range definitions. It also shows the spectral tristimulus functions relevant to the perception of colours, as defined within the standard CIE 1931 (CIE - International Commission on Illumintation, 1931). The bottom chart depicts the photopic luminous efficiency function that describes the overall sensitivity of the eye to light. Further, it contains the melanopic efficiency function (International Well Building Institute, 2022). It describes the sensitivity of melanopsin photoreceptors to light, which is highly important for the body's circadian (day/night) rhythm. More detailed information on the global radiation spectrum is contained in section 4.4.
### 4.2.2. Linearity of light

Throughout this work, linearity of light is assumed, implying that light waves will not interact, i.e. light beams can cross without influencing each other. The superposition principle is a direct consequence of this linearity. It implies that two wave-functions $E_1$ and $E_2$ satisfying Maxwell's equations can simply be added to find their superposition $E_{tot}$, which will itself be a solution of Maxwell's equations, i.e. it will satisfy the derived wave equation, see (31):

$$E_{tot}(x,t) = E_{1}(x,t) + E_{2}(x,t) \qquad (5)$$

Even though this might sound trivial, it is a rather particular property of Maxwell's equations (as well as of other wave equations). It means that the wave's propagation speed is entirely independent of its amplitude. If this does not hold, as is the case for many natural processes, then the principle of superposition cannot be applied. This, in turn, means that each particular wave or wave configuration has to be modelled individually. Such a model can, of course, only be defined for specific configurations considering only a limited number of wave components. Hence, linearity of light is the condition sine qua non, or better, a necessary assumption to perform raytracing on complex objects using a large number of rays. Linearity is assumed in most optical methods and models without being

explicitly mentioned. The related field of optics is referred to as *linear optics*.

There is, however, a field of optics called *nonlinear optics*, which implies that light waves can interact. While Maxwell's equations mathematically lead to perfectly linear wave solutions, the contained electric and magnetic field constants within a medium are, in fact, (weakly) dependent on the field strengths. This means that the induced polarising effects of electromagnetic waves in a medium will alter its dielectric and magnetic susceptibility. Therefore, the presence of multiple electromagnetic waves, or generally stronger electromagnetic fields, will alter the electromagnetic properties of the material, which will usually lead to increased attenuation and some scattering effects. Hence in the real world, light waves do interact. However, usually, i.e. considering typical field strengths of global radiation, these non-linear effects are weak and can therefore be neglected.

For completeness, it should be added that nonlinear effects can even be observed in the absence of matter, i.e. in a vacuum. However, these effects are again significantly weaker and can only be explained based on quantum electrodynamics.

For the practical application of the presented method, linearity and the resulting principle of superposition is essential as it allows decomposing radiation processes in space, angularly, spectrally, as well as in time. Conversely, the individual components can later be superimposed in the evaluation process, considering any specific angular irradiance profiles or time frames.

<figure><img src="image_17.png" /><figcaption>Figure 12 Principle of superposition; "collision" of two wavelets with different wavelengths</figcaption></figure>
### 4.2.3. Coherence and interference

Coherence is an ideal property of waves allowing them to form stationary, or quasi-stationary, interference patterns in space or in time. In the present work, coherent as well as non-coherent superposition of light beams are assumed for modelling different processes. While this may sound inconsistent at first, it is not, as both processes play a role in the real world. A short introduction to coherence is given below to understand the criteria that decide which model has to be applied for a specific process.

It is often assumed that a light beam consists of many largely extended, or infinitely extended, monochromatic sinusoidal electromagnetic plane waves with different wavelengths. If that were the case, numerous interference effects, where light intensity is either amplified or suppressed, could be observed when light is reflected on surfaces or propagates through transparent materials. In fact, such behaviour can be observed when laser light is diffusely reflected off a surface. The micro-roughness of the surface causes specific phase shifts that result in a random granular pattern called speckles. The principle of superposition is applied, as mentioned in the last section, to analyse the electric field formed by two identically polarised waves having the same wavelength and travelling along the same axis. Its superposition is written as:

$$E_{tot}(x, t) = E_{1}(x, t) + E_{2}(x, t) \qquad (6)$$

Using the complex representation, the sinusoidal waves are expressed as:

$$ 
\begin{aligned}
E_{tot}(x,t) &= A_1 \exp[i(kx - \omega t + \delta_1)] + A_2 \exp[i(kx - \omega t + \delta_2)] \\
&= \dots = A_{tot} \exp[i(kx - \omega t + \delta_{tot})]
\end{aligned}
\qquad (7) \\
\begin{array}{l}
A \text{ ... amplitude } [V/m] \\
k \text{ ... wavenumber } [m^{-1}] \\
\omega \text{ ... angular frequency } [s^{-1}] \\
\delta \text{ ... phase offset}
\end{array}
$$

The resulting superposition is, therefore, itself a sinusoidal wave travelling at the same speed. Considering that the intensity (radiant flux) of the waves is proportional to the square of their amplitudes ($I \propto A^2$), the intensity of the superposition can be written as (see, e.g., Hecht, 2015; Bergmann, Schaefer, et al., 2019):

$$I_{tot} = I_1 + I_2 + 2\sqrt{I_1 I_2} \cos(\Delta\delta) \qquad (8)$$

Where $\Delta\delta = \delta_2 - \delta_1$ is the phase difference between the two initial waves. Hence, the resulting intensity of two interfering waves of the same wavelength and intensity $I_1 = I_2 = I_0$ can have any value ranging from 0 (destructive interference) to $4 \cdot I_0$ (constructive interference) depending on the actual phase difference $\Delta\delta$. Equation (8) can further be used to derive the intensity of non-coherent superpositions: for a large number of interfering waves with random phase differences, the average of the cosine in the cross-term of equation (8) will be equal to 0, and therefore the intensities will simply add up.

Hence, to interfere, waves have to maintain a constant phase difference. In reality, they do this only over a certain specific length referred to as coherence length $\Delta l_c$. Since the “ideal” sinusoidal wave with a constant phase only exists within this spatial extent, the term wave packet is used to refer to this elementary light component. The spatially limited coherence is a direct consequence of the temporally limited atomic processes involved in the generation of light. Light in the spectral range of global radiation is created by transitions of the outer electrons.

The wave's spatial and temporal extent, referred to as coherence time $\Delta t_c$, is, on the one hand, connected simply by the speed of light:

$$\Delta l_c = c \Delta t_c \qquad (9)$$
On the other hand, Fourier analysis can be used to show that a spatially limited wave packet can only be formed by a superposition of frequencies limited to a narrow band around the central frequency $v_0$. In fact, the wave packet again consists of numerous photon wave trains having the same spectral distribution. The frequency distribution, which can be measured using a spectrometer, is of Gaussian shape. Its Fourier transform is a sinusoidal wave packet which is as well modulated by a Gaussian envelope. The coherence length corresponds to the width of this Gaussian distribution (see Figure 13).

<figure><img src="image_18.png" /><figcaption>Figure 13 Wave packets in space and time domain and corresponding frequency spectra</figcaption></figure>

Based on Fourier analysis, it is also apparent that the frequency bandwidth $\Delta v$ is reciprocally linked to the coherence time:

$$ \Delta v \propto \frac{1}{\Delta t_c} \qquad (10) $$

Hence, light sources with a narrow frequency distribution, like Lasers, exhibit long coherence times, and their wave packets extend over large distances.

The actual coherence time of a light source is determined by the duration of electronic transitions and is affected by thermal effects. Table 1 gives an overview of coherence lengths and times for different light sources. Since the coherence of global radiation is relevant to the present work, values of three different sources are included. While the widely used values of Hecht et al. (Hecht, 2015) and Saleh et al. (Saleh and Teich, 1991) are simply calculated based on equation (10) and by considering the range of the visible spectrum only, Ricketti et al. (2022) recently have established an empirical method. This evaluation takes into account the entire global radiation spectrum.

<table><caption>Table 1 Overview of coherence time and corresponding coherence lengths of different light sources</caption><thead><tr><th>light source</th><th>coherence time<br>Δt<sub>c</sub></th><th>coherence length<br>Δl<sub>c</sub></th></tr></thead><tbody><tr><td>"White light (400-700 THz)" (Hecht, 2015)</td><td>3 fs</td><td>900 nm</td></tr><tr><td>"Filtered sunlight 400-800nm" (Saleh and Teich, 1991)</td><td>2.67 fs</td><td>800 nm</td></tr><tr><td>"atmospherically filtered sunlight" (Ricketti et al., 2022)</td><td>1.12±0.04 fs</td><td>336 nm</td></tr><tr><td>Light-emitting diode<br>(λ<sub>o</sub> = 1 µm, Δλ<sub>o</sub> = 50 nm) (Saleh and Teich, 1991)</td><td>67 fs</td><td>20 µm</td></tr><tr><td>Low-pressure sodium lamp (Saleh and Teich, 1991)</td><td>2 ps</td><td>600 µm</td></tr><tr><td>Single-mode He-Ne laser (λ<sub>o</sub> = 633 nm)<br>(Saleh and Teich, 1991)</td><td>1 µs</td><td>300 m</td></tr></tbody></table>
As can be seen, the coherence length of global radiation lies in the range of a few hundred nanometers up to a few micrometres if only a limited spectral range is considered, i.e. filtered. This is the essential threshold that determines if a coherent or non-coherent model has to be used. Consequently, in the present work, coherent processes are used for thin-film modelling (see section 4.7) and, of course, implicitly for all Fresnel calculations. In contrast, non-coherent superposition of light is assumed in the thick-slab model (see section 4.3.7), for scattering in the surface-roughness model and generally for the rest of the raytracing process.

Further, it should be noted that the relations above cannot explain how the very distant atomic oscillators spread on the sun's surface, which have to be assumed to be locally incoherent, can create a coherent field at far distances. The phenomenon, which might seem somewhat counterintuitive, requires a deeper and more elaborate model of spatial coherence that cannot be covered here. In particular, the *Cittert-Zernike theorem* (van Cittert, 1934) describes how the covariance of waves coming from individual incoherent sources increases during their propagation to a remote location. Following the theorem, a simple relation can be derived for the coherence of a distant, extended, incoherent light source. In the model, the measure *coherence area* $A_c$ is used. It specifies the surface area in a plane normal to the optical axis in which points show high coherence. For a circular light source, like the sun, the relation is:

$$A_c = \frac{D^2 \lambda^2}{\pi d^2} \qquad (11)$$

$D$ ... distance to light source [m]

$\lambda$ ... wavelength [m]

$d$ ... diameter of light source [m]

By evaluating equation (11), a coherence area of the sun is in the range of a few $10^{-3}$ mm² is obtained. Nonetheless, the value can reach a few square meters for distant stars. The corresponding high coherence of such light is extensively exploited in astronomical research.

### 4.2.4. Polarisation

The electric field of a monochromatic plane wave travelling in the z direction can always be expressed by a vector sum of two orthogonally oriented waves:

$$\vec{E}(z,t) = \vec{E}_x(z,t) + \vec{E}_y(z,t) \\ \vec{E}_x(z,t) = \vec{x} \cdot E_{0,x} \cos(kz - \omega t + \delta_x) \qquad (12)$$

$$\vec{E}_y(z,t) = \vec{y} \cdot E_{0,y} \cos(kz - \omega t + \delta_y)$$

$\vec{E}$ ... electric field vector [V/m]

$\vec{E}_x, \vec{E}_y$ ... x, y components of field vector [V/m]

$E_{0,x}, E_{0,y}$ ... amplitudes [V/m]

$k$ ... wave number ($\omega/c$) [rad/m]

$\omega$ ... angular frequency [rad/s]

$\delta_{x,y}$ ... phase offsets [rad]

Following geometrical considerations, it is obvious that, depending on the phase difference (or phase shift) of the two orthogonally oriented waves ($\Delta\delta = \delta_y - \delta_x$), different propagation modes of the resulting electric field vector $\vec{E}(z,t)$ can be observed. These are referred to as polarisation states. The linear polarisation case is found for $\Delta\delta = 0$ or $\Delta\delta = \pi$. With this phase shift, the harmonic waves $\vec{E}_x(z,t)$ and $\vec{E}_y(z,t)$ will move in phase. Consequently $\vec{E}(z,t)$ will be bound to a single plane, where the angle of the plane is defined by the ratio of $E_{0,x}$ to $E_{0,y}$.

A phase shift of $\Delta\delta = \pi/2$ and $\Delta\delta = 3\pi/2$ will lead to a rotation of the resulting electric field vector $\vec{E}(z,t)$ around the z-axis, forming a helix. If the amplitudes $E_{0,x}$ and $E_{0,y}$ are equal, the projection of this helix will form a circle. Therefore, this polarisation state is referred to as circular polarisation. If the amplitudes are unequal, the z-projection will form an elliptical shape, with $E_{0,x}$ and $E_{0,y}$
defining the dimensions of the ellipse. This state is therefore referred to as *unrotated elliptical polarisation*.

In the general case, covering any other phase difference angles, a rotated elliptical pattern shape is formed. This case is referred to as *elliptical polarisation*. It can be represented as a superposition of linear and circular polarisation. Analysing the effect of the phase shift more closely, it is evident that that any angles in-between $0$ and $\pi$ will lead to a left-hand oriented helix, while angles in-between $\pi$ and $2\pi$ will cause a right-hand rotation of the electric field vector along the z-axis. An overview of all polarization states resulting from the different phase shifts is depicted in Figure 14.

<figure><img src="image_19.png" /><figcaption>Figure 14 Different polarization states resulting from different phase shifts.</figcaption></figure>

While this model helps to understand the concept of polarization, it will, in general, not be applicable for practical modelling. Two essential aspects have to be considered additionally. First, the model relates to the electromagnetic field vectors of monochromatic waves. These fields, determined by the electric amplitudes, cannot be observed directly. Only the associated intensities, being proportional to the squared value of the amplitudes, are accessible for measurement. Second, light consists of a multitude of individual wave packets and wave trains. Hence, a “statistical” approach that uses average values and allows the description of partial polarization is required. Already in the nineteenth century, Sir George Gabriel Stokes developed a model that addresses both issues. This Stokes model is applied throughout the present work. Further details of the model and definitions are presented in section 4.8.

For most light sources, like the sun or gas discharges, the light is emitted by electronic transitions of randomly oriented atomic emitters. While a specific polarisation can be assigned to each wave train, the resulting polarisation components of their superposition, forming the macroscopic light beam, will vanish and leave no preferential orientation. The term *natural light* or *unpolarised light* is commonly used to refer to this kind of radiation.

The ideal “boundary” states of *fully polarised* or *unpolarised light* are important concepts, but practically most often, partially polarised light will be observed. While light that is artificially polarised by a polarisation filter shows a very high degree of polarisation, a number of more natural processes can also lead to a significant partial polarisation of light.
4.2.4.1. "Natural" polarization

Considering solar radiation and following its path starting at its creation, the main mechanisms of “natural” polarisation can be identified. Shortly after the generation of the radiation on the sun, strong magnetic fields on the solar surface cause a partial polarisation of the solar radiation by the so-called *Zeeman effect*. While this effect is weak and has no practical significance for the radiation reaching the earth, it is essential for astronomy, as it allows the determination of strong magnetic field variations (e.g. at sunspots). While solar radiation reaches the earth and travels through its atmosphere, it is subject to *scattering processes*. This means that the radiation excites molecules in the atmosphere, which behave as dipole antennas that radiate secondary fields. While the molecules are randomly orientated, the dipole profile of the reemitted radiation will cause a polarisation effect, as the radiation of vertically oriented dipoles will not reach the earth’s surface. Generally, the scattering profiles will lead to a polarisation maximum at sky regions oriented 90° towards the direction of solar incidence and to polarisation minima in the direction of the sun (Schott, 2009).

Following the path of the solar radiation further, significant polarisation will be observed once the light is reflected off a surface. If the material is transparent, the portion of the beam that is not reflected will be refracted into the medium. This transmitted beam will again be polarised depending on its incidence angle. The significant polarisation of the reflected and the transmitted beam is of high practical relevance, as it will strongly affect the outcome of any further reflection or refraction event. While these effects are disregarded in most raytracing applications, they are captured by the method developed in this work.

4.2.4.2. Polarisation caused by reflection and refraction

As stated above, the polarisation state
of a light beam can be projected on any
arbitrarily oriented orthogonal coordi-
nate system aligned perpendicular to
the direction of propagation. To effi-
ciently model the interaction of a light
beam at a material interface, however,
it is necessary to specify its polarisation
state in a reference system aligned to
the surface. The orthogonal directions
forming this reference system are re-
ferred to as “s” and “p”. The terminol-
ogy originates from the German words
“senkrecht” (perpendicular) and “par-
allel” and refers to the direction of the
electric field vector relative to the sur-
face normal vector.

If unpolarised light hits a surface at an oblique angle, the s-polarised components will be reflected significantly stronger. Conversely, the refracted beam will be predominantly p-polarised. The effect is often explained by the fact that the incoming radiation excites dipole oscillators in the material that cannot emit radiation along their axis. However, this can be shown more generally by solving Maxwell’s equation for material interfaces demanding continuity of the electric and magnetic field at the boundary.

<figure><img src="image_20.png" /><figcaption>*Figure 15 Polarisation caused by reflection and refraction*</figcaption></figure>
The calculation shows that for any dielectric material, there is a specific angle at which the p-polarised component of the reflected beam will completely vanish. Correspondingly, and taking into account energy conservation, the entire p-polarised radiation will be refracted, i.e. transmitted into the medium. This angle is called the *Brewster angle*. Natural, unpolarised light will be entirely *polarised* on reflection at *Brewster angle*. If, however, the incoming beam is entirely p-polarised, no reflection on the surface will occur, and the beam will be transmitted without any loss. The latter effect is exploited in photonics when the output windows of lasers are tilted at Brewster angles to suppress reflection on the windows.

Considering reflection on metallic surfaces, similar relations can be observed, though with significant quantitative deviations. Since the electrons' movement is less restricted in metals, they generally show stronger reflection and absorption but less polarisation. Instead of the *Brewster angle*, a pronounced maximum for the degree of polarisation at a particular angle can be found. However, the reflected beam never reaches total polarisation.

As stated above, the reflectance and transmittance coefficients are always calculated individually for the s- and p-polarised components. To get the values for unpolarised radiation, where the s- and p-components of the incident radiation are equal, the intensities of the reflected and transmitted polarisation components have to be averaged. The total reflectance *R* and transmittance *T* for unpolarised light are therefore simply determined by:

$$R = \frac{R_s + R_p}{2}, \quad T = \frac{T_s + T_p}{2} \qquad (13)$$

*R<sub>s,p</sub> ... reflectance for s/p - polarized wave*

*T<sub>s,p</sub> ... transmittance for s/p - polarized wave*

The reflectance functions for s-/p- and unpolarised incident light are depicted in Figure 16. A window glass (soda-lime glass) surface represents the dielectric case, whereas gold is chosen as a metallic surface. The calculations were performed for a wavelength of 600 nm. As can be seen, the reflection of the dielectric glass surface causes a significant degree of polarisation, reaching a value of 100% at the Brewster angle of 57 degrees. The polarisation of the light refracted/transmitted into the glass is generally lower but increases with the incidence angle to reach a value of 40% at grazing angles. The reflection off the metallic (golden) surface is stronger but generally induces less polarisation. However, the relation of the three components is similar, and while there is no actual Brewster angle for metallic materials, it can be seen that the degree of polarisation still shows a pronounced maximum at an angle of 73 degrees.

<figure><img src="image_21.png" /><figcaption>Figure 16 Angular Fresnel reflection/transmission coefficients for s-/p- and unpolarised incident light (blue dotted, dashed and solid lines) and degree of polarisation for unpolarised light after reflection/refraction (red dashed line)</figcaption></figure>
The polarising effect of reflection can easily be demonstrated if a photo of a glass surface is taken at an angle close to the Brewster angle. Since the reflected light will be almost entirely polarised, a linear polarisation filter in front of the camera can be used to filter these reflections. The effect is illustrated in Figure 17. Since the picture was taken at a close distance, the viewing angle relative to the glass changes; therefore, the filtering is pronounced at the centre of the window only.

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

Figure 17 Photo of polarised reflection on window glass taken behind vertically aligned polarisation filter (bottom) and horizontally aligned polarisation filter (top)
## 4.3. Light-matter interaction overview

### 4.3.1. Scattering

When light, being an electromagnetic wave, reaches or penetrates matter, it starts to interact with the electrical charges of the material. The electromagnetic wave induces dipoles which reradiate an electromagnetic field. The induced field, or secondary wave, superposes with the incident field. Depending on the strength of the interaction and its geometry, the superposition of the fields will cause either reflection, transmission (and refraction) or diffraction of the incident light. Therefore, the most general definition of scattering is a redirection of light that is not absorbed. This definition, of course, includes the effects of refraction and reflection on surfaces. However, the term is often used to only refer to the processes that cause a redirection of light when

it passes through a medium. The most important cases for this type of scattering are referred to as Mie- and Rayleigh-scattering.

Rayleigh scattering occurs when light is scattered on particles that are significantly smaller than its wavelength. In this case, light is scattered equally in forward- and backward directions. Due to the strong inverse wavelength dependence ($I \propto 1/\lambda^4$), Rayleigh scattering is responsible for the blue appearance of the sky but also for the red, yellow and violet colours frequently visible in the sky and clouds, mainly during sunset and sunrise.

Mie scattering occurs when light is scattered on particles having a size close to its wavelength. Unlike Rayleigh scattering, Mie scattering shows dominant forward scattering characteristics and has only a limited wavelength dependence. As the dimensions of water droplets are roughly in the range of the visible light wavelengths, Mie scattering occurs there. Due to the limited wavelength dependence, all colours of the visible light spectrum are scattered equally; therefore, clouds usually appear white or grey.

Since the size of the scattering particle relative to the wavelength $\lambda$ is relevant, it is helpful to define a non-dimensional size parameter $x$ to determine the effective scattering mechanism. For spherical particles with radius $r$, this can simply be defined as (Petty, 2006):

$$x = \frac{2\pi r}{\lambda} \qquad (14)$$

An overview of the relevance of the differenz scattering regimes based on this relation is provided in Figure 18.

<figure><img src="image_23.png" /><figcaption>Figure 18 Effective scattering regimes by particle size and wavelength figure based on (Petty, 2006)</figcaption></figure>
### 4.3.2. Specular reflection

If the incident light interacts with electrons that are relatively free to move, e.g. with unbound electrons in a metal, the electrons will begin to oscillate with a phase difference of 180°, inducing a secondary field. The interference of the incoming light with this secondary field will cause the extinction of the transmitted field and cause a distinct reflection of the incoming light on the surface. Following geometric considerations, it can be shown that the outgoing radiation will leave the surface at the same angle as the incidence occurred. This results in the law of reflection which states that the angle of incidence equals the angle of reflection. Consequently, pronounced specular reflections can be found mainly on the smooth (flat) surfaces of metals or glass. However, there is an important exemption to that rule, as many rough surfaces will exhibit specular reflections for higher incidence angles, known as gloss. A commonly known case is the reflection of sunlight on the rough street pavement close to sunset or sunrise. As a simple explanation, it can be stated that the “effective roughness” of the surface, “seen” from the incident light, decreases with higher angles of incidence for some surface topologies.

While the intensity of the specular reflection is determined by the Fresnel coefficients (see section 4.7), the law of reflection has to be applied to determine the outgoing vector $\vec{v}_o$ of the incoming vector $\vec{v}_i$ and the surface normal $\vec{n}$. Simple vector algebra can be used to derive a computationally very efficient relation. It can easily be comprehended if the tangential ($\vec{v}_{i,t}$) and normal component ($\vec{v}_{i,n}$) of the incoming vector is considered:

$$ 
\begin{aligned}
\vec{v}_o &= \vec{v}_{i,t} - \vec{v}_{i,n} \\
&= (\vec{v}_i - (\vec{v}_i \cdot \vec{n}) \cdot \vec{n}) - (\vec{v}_i \cdot \vec{n}) \cdot \vec{n} \\
&= \vec{v}_i - 2 \cdot (\vec{v}_i \cdot \vec{n}) \cdot \vec{n}
\end{aligned} 
\quad (15) $$

No additional normalisation is required if the incoming vector is normalised. This means that the algorithm is very efficient, as only a few basic arithmetic operations are required.

### 4.3.3. Diffuse reflection

The visible appearance of most objects is dominated by the diffuse reflection of light from their surfaces. While specular reflection relies on the fact that the secondary emitters are regularly aligned so that the induced radiation forms distinct interference patterns with the incident radiation, diffuse reflection represents the more general case of randomly arranged oscillators. The irregularity of the material can either be caused by a rough surface topology or by microscopic irregularities inside the material if the light is able to penetrate below the surface. The latter case is also known as subsurface scattering. It is the relevant mechanism

to describe how paints work: generally, paint is formed of a binder containing microscopic pigments. The binder usually has a low index of refraction to allow light penetration below the surface. The randomly arranged pigments in the binder usually have higher refractive indices in specific wavelength ranges and scatter the incoming radiation in all directions. Usually, multiple

<figure><img src="image_24.png" /><figcaption>Figure 19 Specular and diffuse reflection on car paint (photo by Señor Codo/Wikimedia, used under CC BY 2.0)</figcaption></figure>
scatterings will occur within the paint layer so that the direction of the light exiting the surface is randomly distributed. Corresponding to the random orientation of the directions, also the polarisation axis of the reflected wave trains are randomly oriented so that diffusely reflected light is largely unpolarised. Diffuse reflection usually is accompanied by specular reflection. While the specular reflection occurs in the uppermost regions of the binder or on additional coating layers, the diffuse reflection occurs on the submerged pigments. This behaviour can be modelled in the *RadiCal* method (see Figure 104).

The specular reflections on the binder or coatings are usually not spectrally selective. This means that all colours are reflected with almost equal intensity. In sharp contrast, many pigments exhibit a pronounced spectrally selective behaviour. The combination of the “colourful” specular reflection and “monochromatic” diffuse reflection can be seen in the photo depicted in Figure 19. As discussed in section 4.2.4, a polarisation filter can be used to suppress the specular reflection, whereas the mainly unpolarised diffuse reflection will remain unaffected.

### 4.3.3.1. Lambert reflectance

The Lambertian model is a practically relevant and commonly used definition for an ideal diffusely reflecting surface. It represents the case in which the reflectance of a surface is constant for any given incidence and reflectance angle, implying that objects with Lambertian surfaces appear equally bright from all viewing angles. This is a good approximation for many surfaces and allows a significant level of simplification in many optical models. The characteristics of a Lambertian surface can be shown by deriving the corresponding bi-directional reflectance distribution (BRDF) function $f_r(\vec{\omega}_t, \vec{\omega}_o)$ of the hemispherical-directional reflectance $\rho_d$, commonly also defined as the albedo of a surface. $\rho_d$ is defined as the fraction of light scattered into a direction by a surface illuminated uniformly by hemispheric light (or vice versa):

$$ \rho_d = \int_{\Omega} f_r(\vec{\omega}_t, \vec{\omega}_o) (\vec{\omega}_t \cdot \vec{n}) d\Omega (\vec{\omega}_t) \quad (16) $$

$\vec{\omega}_t$ ... direction of incoming light

$\vec{\omega}_t$ ... direction of outgoing light (viewing)

$\vec{n}$ ... surface normal

Since $f_r(\vec{\omega}_t, \vec{\omega}_o) = \text{const.}$ is demanded, the term can be excluded from the integration so that only the product of the incoming direction with the surface normal vector has to be integrated over the hemisphere.

$$ \rho_d = f_r \int_{\Omega} (\vec{\omega}_t, \vec{\omega}_o) (\vec{\omega}_t \cdot \vec{n}) d\Omega (\vec{\omega}_t) = f_r \pi \quad \rightarrow \quad f_r = \frac{\rho_d}{\pi} \quad (17) $$

The constant *BRDF function* of any Lambertian surface can, therefore, be found by dividing its albedo by $\pi$.
### 4.3.4. Refraction

Refraction represents the change in direction that a light wave suffers on a boundary between two media with different indices of refraction. Since the refraction index measures the slowing of the phase speed of light in a medium, refraction can be explained using the Huygens-Fresnel principle, in which the wavefronts of light are described as a superposition of spherical wavelets.

Figure 20 illustrates how the slower propagation speed (i.e. higher index of refraction) in the second medium will

cause the refracted beam to bend towards the surface normal. The same effect can be observed at beaches where the slower speed of waves in the shallow water will cause a swing of the wavefront towards the shoreline so that the waves will mostly reach them at a perpendicular angle.

Snell's law defines the ratio of the angle of incidence to the angle of refraction:

$$ \frac{\sin \theta_1}{\sin \theta_2} = \frac{n_2}{n_1} \qquad \sin \theta_1 n_1 = \sin \theta_2 n_2 \qquad (18) $$

For polychromatic light, refraction usually leads to a separation of colours as the refraction index is generally a function of wavelength. Typically, shorter wavelengths will suffer stronger refraction.

Based on equation (18), it can be seen that the transmitted beam will be refracted towards the surface normal. If, in contrast, the incidence medium is optically denser ($n_1 > n_2$) the transmitted beam will be refracted away from the surface normal. Consequently, there is a critical angle $\theta_c$ at which $\theta_2$ will reach 90°. At this angle, the transmitted beam will vanish, and all energy will be reflected off the interface.

$$ \sin \theta_c = \frac{n_2}{n_1} \qquad (19) $$

This case is referred to as total internal reflection. Although the reflection behaviour changes here, total internal reflection does not represent a discontinuity, as the intensity of the transmitted beam decreases continuously with the angle to reach 0 at the critical angle.

For a more detailed analysis of refraction and reflection at interfaces, it is necessary to consider the polarisation state of light, as well as the intensities of the beam components. A comprehensive model covering these aspects implicitly has been implemented in the RadiCal method and is presented in detail in section 4.7.

<figure><img src="image_25.png" /><figcaption>Figure 20 Huygens-Fresnel principle</figcaption></figure>
### 4.3.5. Diffraction

Similar to the case of refraction, the effect of diffraction can be explained using the Huygens-Fresnel principle. However, diffraction occurs in the same medium if a light wave passes nearby an obstacle. Diffraction will affect the linear propagation of light and generally causes a bending of the light beam around a barrier. Figure 21 shows how the bending effect can be explained based on the superposition of wavefronts of spherical wavelets.

Diffraction cannot easily be integrated into a raytracing process, as it causes a “bending of light rays”, or, respectively, of fractions of it. This is, of course, beyond the scope of geometrical optics, which is based on the linear propagation of light. Diffraction is, therefore, not contained in the *RadiCal* method. To assess the impact of diffraction or, respectively, its neglect, a short estimation of the diffraction angle for relevant wavelengths and dimensions is provided below.

In order to analyse the impact of diffraction, two different cases have to be distinguished. In the first case, the diffraction caused

by a small opening, a hole in 2D or a slit in 1D, is analysed. In the second case, the bending of light around an edge, i.e. the diffraction pattern caused by a semi-infinite screen is investigated (“knife edge diffraction”).

<figure><img src="image_26.png" /><figcaption>Figure 21 Diffraction explained by Huygens-Fresnel principle</figcaption></figure>

#### 4.3.5.1. Case 1 – Fraunhofer diffraction behind a slit

For this case, the Fraunhofer diffraction model has to be applied, as it describes the far-field diffraction pattern caused by a small finite opening.

In order to estimate the relevance of diffraction, it is sufficient to perform the calculation in one dimension, i.e. for the diffraction at a slit-like opening. The Fraunhofer diffraction model (Hecht, 2015) states that the angular intensity profile of the diffraction pattern behind a slit of dimension *b* can be modelled according to:

$$ \frac{I(\beta)}{I(0)} = \left(\frac{\sin \beta}{\beta}\right)^2 \qquad (20) $$

While $I(0)$ is the central intensity and $I(\beta)$ is the oblique intensity, if the dimensionless length measure $\beta$ is defined as

<figure><img src="image_27.png" /><figcaption>Figure 22 Intensity profile of Fraunhofer diffraction</figcaption></figure>

$$ \beta = \frac{\pi b}{\lambda} \sin\theta \qquad (21) $$

The corresponding intensity profile of the diffraction pattern can be seen in Figure 22. As can be seen, roughly 90% of the total intensity is contained in the central maximum.

In order to obtain a meaningful analysis of the bending angle, an angle $\theta_h$ is defined at which the intensity drops to half the central intensity value.
Equation (20) is numerically solved for $\frac{I(\beta_h)}{I(0)} = 0.5$ and substituted into (21) to find that:

$$ \beta_h = 1.39156 = \frac{\pi b}{\lambda} \sin\theta_h \rightarrow \sin\theta_h = \frac{1.39156 \lambda}{\pi b} \qquad (22) $$

Numerical integration of (20) up to $\beta_h$ reveals that 72% of the total intensity lies within normal incidence and this angle.

Equation (22) is used to calculate the bending angles for half intensity for different wavelengths and slit sizes. The calculation is performed for the boundary wavelengths of the global radiation spectrum, as well as for the boundaries of the visible spectrum. The tangent of the angles, which, of course, coincides with the deflection $d$ of the rays at a distance of 1 m, is depicted in Figure 23. As can be seen for openings of 1 mm, this deflection will range below 1 mm for most of the global radiation spectrum and will only reach 2.2 mm for the longest wavelength of the global radiation spectrum.

<figure><img src="image_28.png" /><figcaption>Figure 23 Tangent of half intensity diffraction angle $\theta_h$ for various slit sizes and wavelengths</figcaption></figure>

### 4.3.5.2. Case 2 – Fresnel diffraction behind a semi-infinite screen

While the Fresnel diffraction model is usually relevant for the near-field, in close proximity to the obstacle, it can be shown that it is also applicable to determine the diffraction at an edge, as in this case, the dimensions of the opening and screen are (semi-)infinite.

The derivation of the relation is elaborate and not covered here; it can be found in (Hecht, 2015). The result of a plane wave hitting an extended obstacle is analysed by calculating the intensity profile behind a semi-infinite screen located in the half-plane $x \le 0$. The lights intensity on a screen at a distance $r_0$ can be described by:

$$ I(x) = \frac{I_0}{2} \left[ \left( \frac{1}{2} - C(v(x)) \right)^2 + \left( \frac{1}{2} - S(v(x)) \right)^2 \right] \\ \text{with } v(x) = x \sqrt{\frac{2}{\lambda r_0}} \qquad (23) $$

$C(v), S(v)$ ... Fresnel integrals
$I(x)$ ... intensity on screen at location $x$

<figure><img src="image_29.png" /><figcaption>Figure 24 Intensity profile resulting from diffraction at edge of semi-infinite screen at 1 m distance</figcaption></figure>

Again, the interference pattern on the screen is calculated for a screen at a distance of 1 m, so that the displacement on the screen corresponds to the tangent of the “bending angle” measured from the
edge at $x = 0$. The resulting intensity distribution is again calculated for the boundary wavelengths of the global radiation spectrum. The result is depicted in Figure 24.

As can be seen, a characteristic intensity profile with two specific features develops: Firstly, it is obvious that some light is bent behind the obstacle into the geometrically shaded area (visible at $x < 0$). Precisely at the edge in between the illuminated and shaded region, the intensity reaches a value $1/4$, since the half-amplitude can be expected there. In the illuminated area subsequent to the edge, the intensity rises to reach a maximum value at a distance of approx. 0.5 mm for the UV radiation (300 nm) and approx. 1.5 mm for the infrared radiation (2500 nm). Subsequently, the intensity profile shows oscillations to finally converge to the expected value of 1 for unobstructed light incidence.

#### 4.3.5.3. Coherence considerations

In the two cases evaluated above, coherent superposition of light waves is assumed. As specified in section 4.2.3, sunlight's temporal and spatial coherence is limited. Hence, the interference patterns of both cases will usually not fully develop as presented in the above charts when sunlight is being used as the source. While the effect can be observed, see, e.g., Greco *et al.* (1997), it is usually limited to minimal angles or induced by very small openings only. Divitt *et al.* (2015) have empirically shown that due to the limited spatial coherence of sunlight, strongly developed interference effects can only be observed for slit distances in the range of a few micrometres and below.

#### 4.3.5.4. Supposed diffraction

Many of the observed optical effects in which small openings (e.g. keyholes) cause significantly broadened light beams are attributed to light diffraction, even in scientific literature. A typical example is shown in Figure 25. In fact, these effects are merely caused by the spatial extension of the sun. The angular solar diameter of 0.5 degrees will cause a beam broadening of 3.5 cm on a perpendicularly oriented screen at a distance of 4 m and significantly larger values on any tilted screen. Complementary, the extended blurry edge of shadows is a consequence of this effect and is as well not related to diffraction. The observation that no particular colouring effect is visible in these blurry edges further supports this: as can be seen in the two cases above, diffraction strongly depends on wavelength; consequently, real diffraction of sunlight will

lead to a dispersive effect and, hence, usually result in a colourful appearance (see Figure 26).

<figure><img src="image_30.png" /><figcaption>Figure 25 Illumination pattern on the floor not caused by diffraction, but by the finite angular diameter of the sun ($\delta_{sun} \approx 0.53^{\circ}$).</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>Figure 26 Examples for true diffraction

left: sunlight diffracted at the side-rail of a blind (linear); centre: view through a textile screen (rectangular openings of approx. 300 µm); right: view through bird feather;

photos by Prof. Dietrich Zawischa (Zawischa, 2023)</figcaption></figure>

4.3.5.5. Conclusions regarding the consideration of diffraction

For light with wavelengths in the range of the global radiation spectrum, “bending effects” caused by diffraction at openings larger than approx. 1 mm are weak for coherent light sources (see Figure 23). This is especially true for the shorter wavelengths, which are energetically more relevant (see Table 2). If the limited coherence of sunlight is additionally considered, the relevance of diffraction effects is further reduced. Generally, it is relevant only for openings with dimensions in the micro-metre range. The potential level of detail in modelling the raytracing targets will usually not allow a significant number of openings with widths or diameters below 1 mm; however, this will not be relevant for most practical applications.

Textile materials or fine metallic grids might, however, pose a significant exception to this. To overcome this issue, macroscopic surface scattering models, which intrinsically model the microscopic diffraction process, can be developed and implemented in the *RadiCal* method.

Note that even if a significant amount of diffraction starts to occur, the direction of maximum intensity for light passing any one- or two-dimensional openings will still correspond to the central, undeflected path as modelled in the geometrical optics approach. Only the laterally scattered fractions are neglected by the linear model.
### 4.3.6. Absorption and attenuation

While light propagates through a medium, it can be observed that its intensity decreases along its linear path. The terms extinction and attenuation are usually used to refer to this process; however, the terminology is ambiguous and depends on the scientific subject.

The reason for the attenuation of light is two-fold. First, the electromagnetic wave interacts with matter and transfers energy to it, known as absorption. Second, light can be partially scattered in the medium due to irregularities in the electromagnetic properties of the medium. The irregularities can be either of microscopic nature with particle sizes below 1µm, e.g. aerosols, or larger, as water droplets in clouds or fat globules in milk (see section 4.3.1).

<figure><img src="image_32.png" /><figcaption>Figure 27 Attenuation in medium according to Beer-Lambert's law</figcaption></figure>

Apparently, there is a fundamental difference between these two processes, as absorption leads to an actual weakening of the electromagnetic field, whereas pure scattering only redirects fractions of the initial light beam. The reason why both processes are often considered together as attenuation is of empirical nature, as the attenuation of light plays an essential role in spectroscopy or astronomy. In such applications, attenuation is simply a measure of the amount of light emitted by the source but not reaching the detector on a linear path, regardless of whether it is scattered away or absorbed. In this work, only attenuation by absorption is considered. There are two reasons for that: firstly, the effect of scattering is negligible for most materials relevant in building science, while absorption ranges over several orders of magnitude and defines whether the material is transparent, partly transparent or opaque. Secondly, to appropriately consider scattering effects, they have to be integrated into the raytracing algorithm, as the energy of the scattered light is not lost but redirected. Hence, energy conservation can only be ensured if these redirection processes are modelled accordingly. While currently not implemented, the necessary scattering algorithm can quite easily be added in future, again based on the Monte-Carlo approach.

In order to quantify attenuation in a medium, a simple relation can be applied. By observation, Pierre Bouger (1760) and Johann Heinrich Lambert (1760) found that equal parts of the same absorbing medium will absorb equal fractions of light:

$$ \frac{dI}{dx} = -\alpha I(x) \qquad (24) $$

I(x) ... irradiance at depth x [W/m²]

α ... attenuation coefficient [1/m]

x ... distance [m]

After integration, the fundamental equation for absorption is found. It was later formulated by Beer (1852) and is now commonly known as Beer-Lambert or Lambert-Beer-Bourguer law:

$$ I(x) = I_0 \cdot e^{-\alpha x} \qquad (25) $$

I₀ ... initial irradiance (at interface) [W/m²]
Sometimes the reciprocal value of the attenuation coefficient $\alpha$ is used. It is referred to as attenuation length $\mu = 1/\alpha$ and defines the length at which the intensity has decreased to the fraction $1/e$ (~37%) of the initial intensity (see Figure 27).

If attenuation by absorption is considered only, Beer-Lambert's law can be derived from an electro-dynamic point of view, i.e. by solving Maxwell's equation for an absorbing medium. This is demonstrated in section 4.6.1. However, looking ahead, the fundamental equation that links the “macroscopic” or empirical quantity attenuation coefficient $\alpha$ to the “atomic” extinction coefficient $\kappa$ (or imaginary part of the refractive index) is already presented here. Since generally $\kappa$ is a function of wavelength $\lambda$, it is stated as:

$$ \alpha(\lambda) = \frac{4\pi}{\lambda} \kappa(\lambda) \qquad (26) $$

### 4.3.7. Thin-films and thick slabs

The term thin-film is generally used to refer to one or more thin, planar and parallel layers in-between two bulk materials, usually air and optically denser material. Typical examples are low-E coatings on glazings, antireflective coatings on sunglasses or oil films on water. The defining characteristic of thin-films is the occurrence of interference effects. The varying refractive indices of the individual layers cause a series of partial reflections at the interfaces. The incoming light is therefore split into multiple paths of light. The superposition of the separate beams determines the thin-film stack's total transmittance or reflectance.

In the thin-film model, it is assumed that the relevant path differences of the beams range below the coherence length of the light, which is in the range of a few micrometres for sunlight (see section 4.2.3). Hence, models of wave optics that consider the waves' phases must be applied to model such thin-film stacks (see section 4.7.1.). Since the optical path length relative to the wavelength determines if constructive or destructive interference occurs, thin-films show distinct spectrally selective properties.

In contrast, the thick slab model has to be applied if the thickness of the layer ranges significantly above the coherence length. In this case, individual beams are no longer able to interfere. Hence, the resulting intensities of the transmitted or reflected beams are simply determined by summation of the intensities of the individual beams.

In the field of glazing, which is of particular importance for this work, the determination of thin-films and thick slabs is relatively straightforward. Coatings on glazings usually consist of layers having only a few nanometers of thickness, while thicknesses of typical glass panes are in the range of a few millimetres. Therefore, coatings must be modelled as thin-films, considering interference effects, while the thick slab model has to be applied for glass panes (see section 4.12).

<figure><img src="image_33.png" /><figcaption>Figure 28 Coherent vs. non-coherent processes</figcaption></figure>
## 4.4. Global radiation

### 4.4.1. Theoretic background and model

The commonly used distinction between longwave radiation (thermal radiation) and shortwave radiation (global radiation or solar radiation) becomes obsolete if the source of the global radiation is considered. Global radiation is, of course, nothing but the thermal radiation of the sun reaching the earth's surface.

Therefore, the irradiance of the global radiation reaching the top of the earth's atmosphere (known as extraterrestrial radiation) can be calculated with surprising accuracy by only using Stefan-Boltzmann's law and simple geometric considerations.

Using Stefan-Boltzmann's law and applying the effective temperature of the sun $T_{solsurf} = 5778$ K (NSSDC, 1998) the radiant exitance of the sun's surface $M_{solsurf}$ can be determined as:

$$M_{solsurf} = \sigma \cdot T_{solsurf}^4 = 6.3206 \cdot 10^7 \text{ W/m}^2 \qquad (27)$$

Considering that the irradiance decreases with $1/r^2$ the radiant exitance at the sun's surface can be scaled by multiplying it by the ratio of the radius of the sun $r_{sol}$ (CNEOS NASA, 2022) to the average earth-sun distance, known as astronomical unit AU (CNEOS NASA, 2022), to determine the irradiance reaching the earth's atmosphere:

$$M_{exaterr.} = M_{solsurf} \cdot \left(\frac{r_{sol}}{AU}\right)^2 = 6.3206 \cdot 10^7 \cdot \left(\frac{6.96608 \cdot 10^8}{1.49598 \cdot 10^{11}}\right)^2 = 1366 \text{ W/m}^2 \qquad (28)$$

In order to get the spectral radiant exitance $M(\lambda, T)$ of extraterrestrial radiation, Planck's law for black-body radiation can be applied:

$$M(\lambda, T) = \frac{2 h c^2}{\lambda^5} \frac{1}{\exp\left(\frac{hc}{k_B \lambda T}\right) - 1} \qquad (29)$$

$T$ ... temperature [K]

$\lambda$ ... wavelength [m]

c ... speed of light [m/s]

$k_B$ ... Boltzmann constant [J/K]

h ... Planck constant [J/Hz]

Relevant for solar calculations, however, is the radiation reaching the earth's surface. In the earth's atmosphere, solar radiation is attenuated, i.e. absorbed and scattered by atmospheric gases. Mostly water vapour, carbon dioxide, ozone and oxygen molecules absorb in this wavelength range (Finlayson-Pitts and Pitts, 2000), leading to the characteristic absorption bands of the global radiation spectrum (see Figure 29). The optical path length through the atmosphere also has to be considered. This is generally done by referring to the air mass coefficient AM, whereby AM1 is used for the path length of a strictly vertical incidence angle (zenith angle). The quantity AM0 is sometimes used to refer to extraterrestrial radiation. Finally, AM1.5 is the commonly used case to characterise average irradiance at a typical incidence angle (48° zenith angle).
<figure><img src="image_34.png" /><figcaption>Figure 29 Components of the ASTM G173 reference spectrum and black-body radiation at 5778 K</figcaption></figure>

The actual spectral distribution and total irradiance fluctuate strongly and depend on many parameters, such as weather, pollution, solar angle etc.

For reference purposes, standard spectra are used. The most commonly used spectrum is provided by NREL and the *American Society for Testing and Materials* and is called ASTM G-173 (ASTM, 2020). The spectrum reflects the irradiance on a 37° tilted surface, facing directly into the sun under clear sky conditions. Alternatively, the international standard EN 410 (2011) used to calculate the total solar transmittance of shaded and unshaded glazings provides the spectral distribution in normalised form. A comparison between the two spectral distributions can be seen in Figure 31.

Using the developed modules (see section 4.4.3), some statistical key figures characterising the global radiation spectrum are derived and presented in Table 2. The analysis was performed using the ASTM G-173 reference spectrum (ASTM, 2020).

<table><caption>Table 2 Statistical analysis of global radiation reference spectrum ASTM G-173</caption><tr><td>median wavelength:</td><td>715 nm</td><td>fraction of UV spectrum (&lt;380 nm):</td><td>3.3 %</td></tr><tr><td>mean wavelength:</td><td>834 nm</td><td>fraction of visible spectrum (380..780 nm):</td><td>53.8 %</td></tr><tr><td>mode wavelength:</td><td>481 nm</td><td>fraction of near-infrared spectrum (&gt;780 nm):</td><td>42.9 %</td></tr><tr><td colspan="4">upper limit wavelength covering 66.6% of total energy:</td><td>875 nm</td></tr><tr><td colspan="4">upper limit wavelength covering 90% of total energy:</td><td>1462 nm</td></tr></table>

It can be seen that the relatively narrow visible wavelength range contributes more than half of the total energy of the entire spectrum and, therefore, also contains the wavelengths with the highest spectral irradiance values. Ultraviolet radiation plays a minor role, as wavelengths within this range are significantly scattered away or absorbed by atmospheric ozone and water molecules (Finlayson-Pitts and Pitts, 2000).
4.4.2. Implementation

For further use, only the normalised spectral distribution of the global radiation is relevant. For later
use in the Monte-Carlo raytracer, it is essential to be able to generate random wavelengths with a
*probability distribution function (PDF)* matching the desired normalised spectrum.

<figure><img src="image_35.png" /><figcaption>Figure 30 Inversion of CDF to provide random wavelengths following the desired spectrum schematic overview (PDF as relative weight of specific wavelength interval)</figcaption></figure>

In order to achieve this, the spectrum's *PDF(λ)* function is generated by normalising the spectrum,
maintaining the discrete interval steps as provided in the specific standard. Also, based on the defined
interval steps, the cumulated distribution function *CDF(λ)* is found, simply by summation. Since this
discrete function cannot be inverted to explicitly achieve *CDF*⁻¹, the inversion is carried out during
the sampling algorithm: firstly, a random number *r* between 0 and 1 is created. The random number
reflecting the *CDF*-function value is used to determine the function's argument, i.e. the associated
wavelength interval. This is performed using a binary search algorithm. Finally, a random wavelength
is selected within the interval, assuming an even distribution. For efficiency reasons, only one random
number generation can be used to cover both steps. The implementation is contained in the function
*getSampleWL* of module *TGlobRad* of the library *rc_globalradiation*.

function TGlobRad.getSampleWL : Double;
var r : Double;
    mid : Integer;
    i1,i2 : Integer;
Label fnd;
Begin
    r:=random; // random number 0..1
    Begin
        i1:=0;
        i2:=maxBin;
        repeat // binary search loop
            mid := (i1+i2) div 2;
            if (r<=CDF[mid])and((mid=0)or(r>CDF[mid-1])) then
                begin // found interval
                    i1:=mid;
                    goto fnd;
                end
                else if (r>CDF[mid]) then i1 := mid+1
                    else i2 := mid-1;
            until i2 < i1;
        end;
        fnd;
        // use remaining fraction of random number as random number to save time
        if i1=0 then r:=(r-0)/(cdf[0]-0) else r:=(r-cdf[i1-1])/(cdf[i1]-cdf[i1-1]);
        Result:=(1-r)*wlmin[i1] + r*wlmax[i1];
    end;
end;
The module can be initialised with different spectral distributions. Currently, the ASTM G173-03 global tilt, ASTMG173-03 direct+circum solar and the EN 410 spectra are available.

```python
Type TGRType = (gsASTMtilt, gsASTMcircDir, gsEN 410);

type TGlobRad = class(Object) { TGlobRad }
private
    PDF : array of Double;
    CDF : array of Double;
    wlMin : array of Double;
    wlMax : array of Double;
    grType : TGRType;

    maxBin : Integer;
    cln : String;
procedure readFromResFile(const fn : String);
procedure AddEntry(const wmin,wmax,pdf_,cdf_ : Double);
function splitDouble(const sept: String; const dmin, dmax: Double; var d: double): boolean;
public
    constructor Create(const grType_ : TGRType);
    function getSampleWL : Double;
    function getSampleWLat(const cdfV: Double): Double;
    destructor Destroy; override;
end;
```

The desired spectrum (`gsASTMtilt`, `gsASTMcircDir`, `gsEN 410`) is specified in the call of the constructor method `Create` when an instance of the class is created.

The additionally included method `getSampleWLat` function allows specifying a value for the CDF parameter instead of using a random number. Hence, the function can be used to pull a set of energetically equidistant sampling wavelengths of the global radiation spectrum. The method is applied in the inversion algorithm to determine refractive index functions (see section 4.12).
4.4.3. Testing and validation

In order to validate the implementation, one billion random wavelength samples are generated using
the `getSampleWL` function (the calculation time required for this is approx. one minute on a single
CPU). The distribution of the sampled data is visualised as a histogram with 5 nm intervals and
compared against the normalised reference spectrum ASTM G173-03. Figure 31 shows the retrieved
results for the `gsEN 410` and `gsASTMtilt` spectrum, as well as the normalised ASTM G173-03 spectrum
for reference. In general, the spectra show similar distributions; however, the ASTM spectra offer a
higher resolution, as the EN 410 spectral distribution is provided in intervals ranging from 20 nm to
100 nm only. For most practical cases, the higher resolution will not be relevant, however.

<figure><img src="image_36.png" /><figcaption>Figure 31 Validation result for global radiation sampling module.
grey: spectral distribution according to ASTM G173-03,
red: sampled data using gsEN 410 mode, , blue: sampled data using gsASTMtilt mode.</figcaption></figure>

As can be seen in Figure 31, the implemented code is able to generate wavelength samples follow-
ing the distribution of the global radiation spectrum. This is essential for its later application and
proves the validity of the algorithm.
# 4.5. Solar position algorithm

## 4.5.1. Theoretic background and model

Numerous algorithms are available for the purpose of calculating the solar position. The accuracy and computation time of the available algorithms vary greatly, however. In the present work, the SG2 algorithm of Ph. Blanc and L. Wald (2012) has been implemented. The SG2 algorithm is a numerically significantly optimised version of the SPA algorithm provided by NREL (Reda and Afshin, 2005). The SPA algorithm, in turn, is based on the astronomical algorithms provided by Jean Meeus (1998) in his book first published in 1991. The SPA algorithm, based on Meeus work, has long been the gold standard as Reda and Afshin specify a maximum uncertainty of ±0.0003° for the years ranging from 2000 b.c. to 6000 a.d.

Based on Blanc and Wald's numerical optimisations, the SG2 algorithm is several orders of magnitude faster than the original SPA algorithm and provides a maximum uncertainty of the solar vector of less than ±0.0025°. This value still exceeds by far the accuracy required for most practical applications, whereas computational performance is crucial for many practical purposes. For this reason, the SG2 algorithm has been implemented here.

## 4.5.2. Implementation

For the actual implementation in Pascal, the original publication (Blanc and Wald, 2012) and a validated implementation in C and Python provided by PSL university (O.I.E. Mines-ParisTech, 2022) have been used.

The code has been implemented in the library `rc Sunpos sg2`. It encapsulates the class `TSunposIrradiance` containing the following methods:

```pascal
Type
  TSunposIrradiance = class(igoObject) { TSunposIrradiance }
    private
      seasonalSolarConst: array[1..366] of Double;
      approx_R_1: TApproxSet;
      approx_L_10: TApproxSet;
      approx_Dpsi_1: TApproxSet;
      approx_Epsilon_1: TApproxSet;
      approx_nu_0_0: TApproxSet;
      approx_M_0: TApproxSet;
      approx_Dtau_cst: Double;
    function SetApprox10(v1,v2,v3,v4,v5,v6,v7,v8,v9,v10: Double): T10DoubleArray;
    function SetApprox10reci(v1,v2,v3,v4,v5,v6,v7,v8,v9,v10: Double): T10DoubleArray;
    function sg2_CALCCHAOS(j: Double; approx: TApproxSet): Double;
    function sg2_jd_to_ymd(jd: Double): Tymd;
    function sg2_setRange(x: Double): Double;
    function sg2_setRange(x: Double;mode: Integer): Double;
    function sg2_setRange(x: Double;mode: Integer; K: Double): Double;
    function sg2_refractCorrection(gamma_S0, P, T: Double; method: TRefractMethod)
      : Double;
    function isEmpty(d: Double): Boolean;
    function sg2_CALCDELTA(y: Double): Double;
    function sg2_ref_ellps(eT: TELLPSType): TELLPS;
    function airMassKastenSandia(const costhetaz_, thetaz, locPressHpa: Double): Double;
    procedure preCalcSeasonalSolarConst;
    procedure InitApprox;
  public
    constructor Create;
    destructor Destroy; override;
    // solar position functions (SG2)
    function calcSunPos(j_, lat, lon, h: Double; var param_opt: TSG2optParams)
      : TsolarPos;
    function JulianD(y, m, d: Integer; daytime: Double; timezone: Double): Double;
    function getJulianDiff(time0, timezone: Double): Double;

    // irradiance functions
    function getSolarConstGround(const elevationInRad: Real): Real;
    function getSeasonalExaSolarConst(const dt: TDateTime): Double;
```
```javascript
procedure getPerezF1F2(const tim, locPressHPa: Double;
                       const DHI, DNI, thetaz: Double; var F1, F2, b: Double);
procedure getPerezF1F2interpol(const tim, locPressHPa: Double;
                                 const DHI, DNI, thetaz: Double; var F1, F2, b: Double);
function getPerezIrradiance(const dt: Double; const DNI_, DHI_, GHI_: Double;
                             locAlbedo, locPressHPa, locLat, locLon, locAlt: Double;
                             refPlaneV: Tvect3d; const refPlaneCosBeta, refPlaneSinBeta: Double;
                             UsePerezInterpolation: boolean; sunOpts: TSG2optParams;
                             irrDat: TIrradianceData): Boolean;
function getReindlIdiffHorizFractionfd(const kt, sunEl, TaCelsius, humi: Double): Double;
function getReindlIdiffHorizFractionSimplefd(const kt, sunEl: Double): Double;
end;
```

The key method is the `calcSunPos` function, which returns the solar angles for a specified date. The function `JulianD` allows converting standard date and time values into a single numeric Julian-date value (of type `double`). Apart from the SG2 code, the class also contains functions for irradiance evaluation. These are covered in section 6.2.5.

### 4.5.3. Testing and validation

The validation has been done with multiple comparisons against free online tools using the SG2 algorithm. In particular, however, the validation example provided in the appendix of the SPA paper (Reda and Afshin, 2005) is used here to demonstrate the validity of the implementation.

<caption>Table 3 Input parameters for the validation example
taken of appendix A of NREL/TP-560-34302 (Reda and Afshin, 2005)</caption>

#### A.5. Example

- Date = October 17, 2003.

- Time zone(TZ) = -7 hours.

- Latitude = 39.742476°.

- Elevation = 1830.14 m.

- Time = 12:30:30 Local Standard Time (LST).

- Longitude = -105.1786°.

- Pressure = 820 mbar.

- Temperature = 11°C.

The code used to perform the solar position calculation for the given validation parameters in the present implementation is:

```javascript
// Default values for optional SG2 parameters
Const DefaultSG2opts: TSG2optParams =
  (TZ: 0;
   P: 1000;
   T: 10;
   ellps:etSPA;
   delta_t: 1.7E300;
   alpha: 0;
   beta: 0;
   atmos_refract: 0.5667;
   sun radius: 0.2667;
   gamma_S_sunrise: -0.8333;
   corr refract: true;
   corr refract_method: rmSAE;
  );

// TEST CASE of Appendix A NREL NREL/TP-560-34302 "SPA PAPER"
lat:=39.742476; lon:=-105.1786;
h:=1830.14;
specialOpts:=DefaultSG2opts;
specialOpts.P:=820;
specialOpts.T:=11;
sp:=sunPos calcSunPos(sunPos.JulianD(2003,10,17,EncodeTime(12,30,30,0), -7 {tz}),
                  lat,lon,h,specialOpts);
Log('Solar position'+#9
+format('azimuth:'+'#9'%'0.5f'%'#9'elevation:'#9'%'0.5f'%',[sp.alpha/Pi*180,sporrHO/Pi*180]]);
```

Table 4 shows the calculation result of the *RadiCal* implementation compared to the solar position stated in the SPA paper. As can be seen, the deviation ranges well below the specified uncertainty.
<table><caption>Table 4 Comparison of solar position: SG2 RadiCal implementation vs. SPA algorithm</caption><thead><tr><th></th><th>solar azimuth angle</th><th>solar elevation angle</th></tr></thead><tbody><tr><td>SG2<br>(RadiCal implementation)</td><td>194.33969 °</td><td>39.88815 °</td></tr><tr><td>SPA</td><td>194.34024 °</td><td>39.88838 °</td></tr><tr><td><strong>difference</strong></td><td><strong>-0.00055 °</strong></td><td><strong>-0.00023 °</strong></td></tr></tbody></table>
# 4.6. Complex-valued index of refraction functions

## 4.6.1. Theoretic background and model

The following paragraphs cover some fundamentals of physical optics relevant to the implementation. The author chose the form and order of the relations presented so that their application in the later sections seems logical. There are, however, many approaches and ways to derive and formulate the presented relations. For further information, please refer to the standard works of optical theory, e.g. Hecht (2015), Bergmann-Schaefer (2019) or Saleh and Teich (1991).

Reformulating Maxwell's equations and applying them to an uncharged medium, the following differential equation is found:

$$ \frac{\partial^2 \vec{E}}{\partial x^2} + \frac{\partial^2 \vec{E}}{\partial y^2} + \frac{\partial^2 \vec{E}}{\partial z^2} = \mu_r \epsilon_r \mu_0 \epsilon_0 \frac{\partial^2 \vec{E}}{\partial t^2} + \mu_r \mu_0 \sigma \frac{\partial \vec{E}}{\partial t} \quad (30) $$

* $\vec{E}$ ... electric field vector [V/m]
* $x, z, y$ ... cartesian coordinates [m]
* $t$ ... time [s]
* $\sigma$ ... electrical conductivity [S/m]

* $\mu_r$ ... relative magnetic permeability [1]
* $\epsilon_r$ ... relative electrical permittivity [1]
* $\mu_0$ ... absolute/vacuum magnetic permeability [H/m]
* $\epsilon_0$ ... absolute/vacuum electric permittivity [F/m]

The last term of the equation represents a damping force that leads to an exponential attenuation in the medium if the conductivity is not equal to zero. Following Clark Maxwell's approach leading to the pivotal discovery (see section 4.1), it is obvious that rest of the equation represents the classical wave equation:

$$ \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} = \frac{1}{v^2} \frac{\partial^2 u}{\partial t^2} \quad (31) $$

* u ... scalar function in time and space

By comparison of the two equations, the phase speed of light $v$ in the medium can be found as:

$$ v = \frac{1}{\sqrt{\mu_r \epsilon_r \mu_0 \epsilon_0}} \quad (32) $$

Introducing the (measured) constant $c_0$ for the speed of light,

$$ c_0 = \frac{1}{\sqrt{\mu_0 \epsilon_0}} = 2,99792458 \cdot 10^8 \text{ m/s} \quad (33) $$

a definition for the *index of refraction* $n$ is established:

$$ n = \frac{c_0}{v} = \sqrt{\mu_r \epsilon_r} \quad (34) $$

The index of refraction represents, therefore, simply the ratio of the vacuum speed of light to the speed of light in the medium.

In a very efficient - but surprisingly quite rarely applied - approach, this scalar definition of the index of refraction $n$ can be extended to the more general concept of a *complex index of refraction* $\tilde{n}$, which also includes an imaginary part, often referred to as the *extinction coefficient* $\kappa$. Since these coefficients show, in general, a dependence on the wavelength $\lambda$ [nm], it is stated as:

$$ \tilde{n}(\lambda) = n(\lambda) - i\kappa(\lambda) \quad (35) $$

In order to show this definition's usefulness but also to link the extinction coefficient to an empirically accessible quantity, it is substituted into the plane wave solution. It can be shown that a plane wave
is a solution to equation (30). In exponential form, a plane wave travelling in the x-direction can be written as:

$$ \vec{E} = \vec{E}_0 e^{i(\omega t - kx)} \qquad (36) $$

ω... angular frequence [rad/s]
k ... wavenumber

A complex wave-number $\tilde{k}$, complex phase speed $\tilde{v}$ is introduced, and equations (34) and (35) are used to bring the complex index of refraction into the wave equation:

$$ \tilde{v} = \frac{c_0}{\tilde{n}} \quad \rightarrow \quad \tilde{k} = \frac{\omega}{\tilde{v}} = \frac{\omega \tilde{n}}{c_0} = \frac{\omega}{c_0} (n - i\kappa) \qquad (37) $$

$$ \vec{E} = \vec{E}_0 e^{i(\omega t - \tilde{k}x)} = \vec{E}_0 e^{i\left(\omega t - \frac{\omega}{c_0}nx\right)} - \frac{\omega}{c_0 \kappa x} \qquad (38) $$

The real part of this equation equals:

$$ E(x) = E_0(x) \cdot e^{-\frac{\omega}{c_0 \kappa x}} \cdot \cos\left(\omega\left(t - \frac{n}{c_0}x\right)\right) \qquad (39) $$

The solution now reflects a sinusoidal plane wave with a phase speed of $c_0/n$ (cosine term), which is exponentially attenuated in the medium (exponential term).

At this point, the concept of attenuation and Beer-Lambert's law (see section 4.3.6) is recalled:

$$ I(x) = I_0 \cdot e^{-\alpha x} \qquad (25) $$

I(x) ... irradiance at depth x [W/m²]

I₀ ... initial irradiance (at interface) [W/m²]

α ... attenuation coefficient [1/m]

Considering that irradiance is proportional to the square of the field amplitude $E$, and by comparing the exponential decay term of (39) with (25), it becomes obvious how the imaginary part of the complex index of refraction (or extinction coefficient) $\kappa$ can be determined empirically by attenuation measurements:

$$ \begin{gathered} I(x) \propto E(x)^2 \propto e^{-2\frac{\omega}{c_0}\kappa x} \rightarrow \\ \alpha = \frac{2\omega}{c_0}\kappa = \frac{4\pi}{\lambda}\kappa \rightarrow \kappa = \frac{\alpha\lambda}{4\pi} \end{gathered} \qquad (40) $$

For further practical application in this work, equations (35), (25) and (40) are most relevant. In the *RadiCal* approach, all materials are described by their spectrally resolved, complex-valued index of refraction functions. The two components of this function, the real part, representing the *index of refraction* $n(\lambda)$ and the imaginary part reflecting the *extinction coefficient* $\kappa(\lambda)$, are not independent of each other as both are results of the electromagnetic interaction of light with the medium, affecting its “speed $n$” and its “loss $\kappa$”. The fundamental link between them is provided in form of the *Kramers-Kronig relation* (de L. Kronig, 1926; Kramers, 1927; Bohren, 2010). The connection between them also leads to *dispersion relations* that are essential for measuring and describing these properties. Dispersion relations for many types of materials can e.g. be found in the *Handbook of optical materials* (Weber M. J. et al., 2003).
### 4.6.2. Implementation

The implementation in the library `rc_refractivityIndex` features a base class `TRIbaseclass` that provides the $n$ and $\kappa$ values for any given wavelengths through the function `getnk(const wl: Double; var n, k: Double)`. By applying polymorphism, an abstract base class is inherited to child classes, which carry the specific definitions. Table 5 provides an overview of the currently implemented classes and the definitions used to model the $n$ and $\kappa$ values as function of wavelength.

<table><caption>Table 5 Currently implemented descendant classes of TRIbaseclass</caption><thead><tr><th>class name</th><th>n - definition</th><th>κ - definition</th></tr></thead><tbody><tr><td>TRIconstNK<br>"simple"</td><td>n ≡ n<sub>0</sub></td><td>k ≡ k<sub>0</sub></td></tr><tr><td>TRIsplineNK<br>"splines"</td><td>cubic spline defined with N<sub>n</sub> knots<br>(λ<sub>i</sub>, n<sub>i</sub>)</td><td>spline defined with N<sub>k</sub><br>knots<br>(λ<sub>i</sub>, κ<sub>i</sub>)</td></tr><tr><td>TRIsimpleDispersion-AndKspline<br>"simple glass"</td><td>n = C<sub>1</sub> + C<sub>2</sub>λ<sup>2</sup> + C<sub>3</sub>/λ<sup>2</sup></td><td>spline defined with N<sub>k</sub><br>knots(λ<sub>i</sub>, κ<sub>i</sub>)</td></tr><tr><td>TRIformula1<br>"Sellmaier<br>dispersion"</td><td>n<sup>2</sup> - 1 = C<sub>1</sub> + (C<sub>2</sub>λ<sup>2</sup>/λ<sup>2</sup> - C<sub>3</sub><sup>2</sup>/λ<sup>2</sup> - C<sub>5</sub><sup>2</sup>/λ<sup>2</sup> - C<sub>7</sub><sup>2</sup>/λ<sup>2</sup> - C<sub>9</sub><sup>2</sup>/λ<sup>2</sup> - C<sub>11</sub><sup>2</sup>/λ<sup>2</sup> - C<sub>13</sub><sup>2</sup>/λ<sup>2</sup> - C<sub>15</sub><sup>2</sup>/λ<sup>2</sup> - C<sub>17</sub><sup>2</sup>/λ<sup>2</sup>)</td><td>tabular with<br>linear interpolation</td></tr><tr><td>TRIformula2<br>"Sellmaier-2<br>dispersion"</td><td>n<sup>2</sup> - 1 = C<sub>1</sub> + (C<sub>2</sub>λ<sup>2</sup>/λ<sup>2</sup> - C<sub>3</sub><sup>2</sup>/λ<sup>2</sup> - C<sub>5</sub><sup>2</sup>/λ<sup>2</sup> - C<sub>7</sub><sup>2</sup>/λ<sup>2</sup> - C<sub>9</sub><sup>2</sup>/λ<sup>2</sup> - C<sub>11</sub><sup>2</sup>/λ<sup>2</sup> - C<sub>13</sub><sup>2</sup>/λ<sup>2</sup> - C<sub>15</sub><sup>2</sup>/λ<sup>2</sup> - C<sub>17</sub><sup>2</sup>/λ<sup>2</sup>)</td><td>tabular with<br>linear interpolation</td></tr><tr><td>TRIformula3<br>"polynomial"</td><td>n<sup>2</sup> = C<sub>1</sub> + C<sub>2</sub>λ<sup>2</sup>C<sub>3</sub> + C<sub>4</sub>λ<sup>2</sup>C<sub>5</sub> + C<sub>6</sub>λ<sup>2</sup>C<sub>7</sub> + C<sub>8</sub>λ<sup>2</sup>C<sub>9</sub> + C<sub>10</sub>λ<sup>2</sup>C<sub>11</sub> + C<sub>12</sub>λ<sup>2</sup>C<sub>13</sub> + C<sub>14</sub>λ<sup>2</sup>C<sub>15</sub> + C<sub>16</sub>λ<sup>2</sup>C<sub>17</sub></td><td>tabular with<br>linear interpolation</td></tr><tr><td>TRIformula4<br>"RefractiveIndex.INFO"</td><td>n<sup>2</sup> = C<sub>1</sub> + (C<sub>2</sub>λ<sup>2</sup>C<sub>3</sub>/λ<sup>2</sup> - C<sub>4</sub><sup>2</sup>C<sub>5</sub>/λ<sup>2</sup> - C<sub>6</sub>λ<sup>2</sup>C<sub>7</sub>/λ<sup>2</sup> - C<sub>8</sub><sup>2</sup>C<sub>9</sub>/λ<sup>2</sup> + C<sub>10</sub>λ<sup>2</sup>C<sub>11</sub> + C<sub>12</sub>λ<sup>2</sup>C<sub>13</sub> + C<sub>14</sub>λ<sup>2</sup>C<sub>15</sub> + C<sub>16</sub>λ<sup>2</sup>C<sub>17</sub></td><td>tabular with<br>linear interpolation</td></tr><tr><td>TRIformula5<br>"Cauchy"</td><td>n = C<sub>1</sub> + C<sub>2</sub>λ<sup>2</sup>C<sub>3</sub> + C<sub>4</sub>λ<sup>2</sup>C<sub>5</sub> + C<sub>6</sub>λ<sup>2</sup>C<sub>7</sub> + C<sub>8</sub>λ<sup>2</sup>C<sub>9</sub> + C<sub>10</sub>λ<sup>2</sup>C<sub>11</sub></td><td>tabular with<br>linear interpolation</td></tr><tr><td>TRIformula6<br>"Gases"</td><td>n - 1 = C<sub>1</sub> + (C<sub>2</sub>/C<sub>3</sub> - λ<sup>-2</sup>/λ<sup>-2</sup>) + (C<sub>4</sub>/C<sub>5</sub> - λ<sup>-2</sup>/λ<sup>-2</sup>) + (C<sub>6</sub>/C<sub>7</sub> - λ<sup>-2</sup>/λ<sup>-2</sup>) + (C<sub>8</sub>/C<sub>9</sub> - λ<sup>-2</sup>/λ<sup>-2</sup>) + (C<sub>10</sub>/C<sub>11</sub> - λ<sup>-2</sup>/λ<sup>-2</sup>)</td><td>tabular with<br>linear interpolation</td></tr><tr><td>TRIformula7<br>"Herzberger"</td><td>n = C<sub>1</sub> + (C<sub>2</sub>/λ<sup>2</sup> - 0.028/λ<sup>2</sup>) + C<sub>3</sub>(1/(λ<sup>2</sup> - 0.028)<sup>2</sup>) + C<sub>4</sub>λ<sup>2</sup> + C<sub>5</sub>λ<sup>4</sup> + C<sub>6</sub>λ<sup>6</sup></td><td>tabular with<br>linear interpolation</td></tr><tr><td>TRIlinTabN</td><td>tabular with<br>linear interpolation</td><td>tabular with<br>linear interpolation</td></tr><tr><td>TRIlinTabNK</td><td>tabular with<br>linear interpolation</td><td>tabular with<br>linear interpolation</td></tr></tbody></table>
The class `TRIconst` is used to create simplified generic materials; the two spline classes are used to approximate refractive index functions of transmittance and reflectance spectra by inversion (see section 4.12). The rest of the classes are used to import data provided by data sources in various formats. In the flexible implementation, only terms with defined values for the constants $C_i$ are evaluated. Therefore, the omission of constants allows an even wider range of potential dispersion relations. The primary source for refractive index data is the excellent website and database `RefractiveIndex.INFO` (Polyanskiy, 2022). The author M.N. Polyanskiy collects and provides data for a wide range of materials and spectral ranges. The information is provided for free access using YAML (Ben-Kiki, Evans and döt Net, 2021) format. The implemented `TRIclasses` allow a direct import of the provided data in the YAML files for use with the `RadiCal` method.

### 4.6.3. Testing and validation

Each of the implemented classes has been tested by the import of sample data. The spectrally resolved $n$ and $\kappa$ values were then calculated based on the functions implemented in the specific classes. The results were plotted and compared with the charts and data provided on `RefractiveIndex.INFO`.

Spectra in the global radiation wavelength range of relevant materials, imported to `RadiCal` in various formats, are presented in Figure 32.

<figure><img src="image_37.png" /><figcaption>Figure 32 Overview of refractive index functions in the global radiation spectral range imported and exported with the implemented classes

data sources: `RefractiveIndex.INFO` (Polyanskiy, 2022)</figcaption></figure>