9 The Exchange Factor Transformation

9.1 Introduction

Similarly to Hottel’s method, outlined in chapter Chapter 2, the exchange factor transformation of the present work starts from the exchange area matrix. But unlike Hottel’s method, which works directly with the exchange area matrix, the proposed transformation uses its non-dimensional counterpart, the exchange factor matrix, where each entry of the exchange area matrix has been divided by its equivalent area of emission. Before solving a specific problem, the proposed method apply a matrix transformation to augment the exchange factor matrix with additional information about the steady state of reflection and scattering in the domain.

The transformation of the exchange factor matrix into an augmented form is partly inspired by the “F-hat” method developed by Beckman (Beckman 1968), which extends the view factor matrix for a transparent enclosure to include reflections. Further analysis revealed that the ‘F-hat’ method only includes single reflections and does not account for multiple reflections. The effective reflection-scattering coefficients of multiple reflection-scattering processes on general domains are functions of extinction and geometry, which seemingly made the problem impossible to solve. However, as will be demonstrated in this book, the problem is neither intractable nor impossible to solve. In fact, it has a closed form analytical solution.

Given the availability of an analytically derived matrix of exchange factors, generally referred to as the matrix \mathbf{F}, the entire solution is analytical. Thanks to Narayanaswamy (Narayanaswamy 2015), who found a general analytical solution to equation (Eq. 2.1) with \overline{t}(S)=1, this is possible for any convex enclosure composed of polygon surfaces and a transparent medium, giving perfect solutions to multiple reflection systems, which would otherwise be represented by coupled Fredholm equations of the second kind (Minkowycz et al. 1988), the single assumption being that angular distributions of emission and reflection are identical.

For problems involving participating media (\overline{t}(S)<1), analytical solutions to the multi-dimensional integrals of equations (Eq. 2.1)-(Eq. 2.3) have yet to be found. In this case, one might employ Monte Carlo methods to obtain near-analytical estimates of \mathbf{F}. To ensure that \mathbf{F} satisfies both The First and The Second Law of Thermodynamics, smoothing should be employed, such as for example the method proposed by Daun et al. (Daun, Morton, and Howell 2005), or preferably the method proposed in this work.

The solution method proposed in this work is general and, furthermore, it is flexible. In fact, a key strength of exchange factor-based methods for solution of the RTE, is that they allow the formidable monolithic problem of energy transfer by an electromagnetic wave field, to be broken down into two conceptually separate and more manageable parts: The geometry and extinction-based ray propagation part and the system-wide energy interaction part. If \mathbf{F} is obtained through Monte Carlo ray tracing methods, any number of relevant optical phenomena may be modelled, such as refraction and diffraction, refraction being significant in atmospheric research, due to the variable refractive index. This can be achieved by solving differential equations for each ray (Born and Wolf 1980). Since the proposed transformation records only ray trajectory endpoints, it is in some sense blind to the phenomena which govern the trajectories, which explains why the matrix transformations of the present work will preserve any modelled ray trajectory properties.

For the proposed transformation, since \mathbf{F} refers only to the first interaction, Monte Carlo methods used to obtain \mathbf{F} can be truncated after this initial interaction, making them computationally less demanding than basic Monte Carlo implementations that must track multiple reflection-scattering events (Howell and Daun 2021). This will be increasingly true for higher single reflection-scattering coefficients due to the increased length of the ray paths. Instead, the proposed transformation performs the multiple reflection-scattering ray tracing analytically using solution of linear systems. This approach offers key advantages:

It significantly improves computational efficiency compared to basic ray tracing methods for multiple reflection-scattering,
It provides analytical ray tracing accuracy with quantifiable bounds on the necessary numerical precision,
The precision requirements are expressed as a function of both the geometry and extinction (described by \mathbf{F}) and the reflection and scattering properties of each element.

This enables problem-specific numerical precision tuning which is readily implementable in modern numerical programming languages (IEEE Standard for Floating-Point Arithmetic 2019).

9.2 Methodology

The methodology is split into three parts: First, the exchange factor matrix is formally defined, next, the exchange factor transformation is defined, and lastly, the transformed exchange factors are applied to solution of the RTE.

9.2.1 Definition of the Exchange Factor Matrix

The proposed transformation requires as an input a general description of the radiative connectivity of the domain in the form of the exchange factor matrix \mathbf{F}. The exchange factor matrix \mathbf{F} may be obtained analytically or numerically. All of the subsequent operations of the proposed transformation are analytical. The (m+n)\times(m+n) row-stochastic blocked first interaction exchange factor matrix \mathbf{F} for a system of m bounding surfaces and n gas volume elements is defined as:

\mathbf{F} = \begin{bmatrix} \mathbf{F}_\mathrm{ss} & \mathbf{F}_\mathrm{sg} \\ \mathbf{F}_\mathrm{gs} & \mathbf{F}_\mathrm{gg} \end{bmatrix}

where \mathbf{F} is row stochastic to ensure that The First Law of Thermodynamics is satisfied and the entries of \mathbf{F} are the fractions of energy emitted by the emitter of row i which has its first interaction with the element of column j. The entries of \mathbf{F} can be defined from equations (Eq. 2.1)-(Eq. 2.3), by normalizing with the equivalent total capacity for emission and reflection-scattering of the emitter zone (and using the extinction coefficient \beta in place of the absorption coefficient \kappa):

F_{s_j s_k} = \frac{\overline{s_j s_k}}{(\rho_j+\varepsilon_j)A_j}, \quad F_{s_k g_\gamma} = \frac{\overline{s_k g_\gamma}}{(\rho_k+\varepsilon_k)A_k}, \quad F_{g_\gamma s_k} = \frac{\overline{g_\gamma s_k}}{4\beta V_\gamma}, \quad F_{g_\gamma g_{\gamma*}} = \frac{\overline{g_\gamma g_{\gamma*}}}{4\beta V_\gamma}.

where \rho+\varepsilon=1 where \rho is the single reflectivity and \varepsilon is the emissivity and for volumes of participating media, \kappa+\sigma_s=\beta where \sigma_s is the single scattering coefficient. Additionally, \mathbf{F} satisfies The Second Law of Thermodynamics (reciprocity) (Daun, Morton, and Howell 2005):

\mathbf{E} \mathbf{F} = \mathbf{F}^T \mathbf{E} \tag{9.1}

where \mathbf{E} is a diagonal matrix of equivalent emission areas with \mathbf{E}_{ii}=(\rho_i+\varepsilon_i) A_i for i\leq m and \mathbf{E}_{ii}=4 \beta_i V_i for i>m, where the full capacity to emit and reflect-scatter is used.

Furthermore, in the general case, \mathbf{F} will have some zero entries, i.e., physical elements without direct view of each other, meaning \mathbf{F} is not merely a positive matrix but a non-negative matrix. Therefore, for the proofs of this book to guarantee convergence, \mathbf{F} must be irreducible (Horn and Johnson 2013), where irreducibility corresponds to the physical requirement that the domain contains no isolated regions, that is, for any two elements i and j, there exists a sequence of direct radiation exchanges i \to k_1 \to k_2 \to \ldots \to j connecting them. This requirement can always be satisfied by treating isolated regions as separate radiative transfer problems.

For the present transformation, \mathbf{F} should be derived exactly as if dealing with pure absorption, but always using the total capacity to emit and reflect-scatter instead of the absorption coefficient. Note that extinction \beta does not need to be uniform across the domain, but should be uniform for each element. In the case of non-uniformity, \mathbf{F} must be obtained using the same non-uniform distribution of \beta. Then after \mathbf{F} has been obtained, any part of the capacity to absorb may be converted into a capacity for reflection-scattering, meaning a single \mathbf{F} covers all combinations of \kappa+\sigma_s which sum to \beta.

9.2.2 The Exchange Factor Transformation

The column constant (m+n)\times(m+n) blocked single reflection-scattering matrix \mathbf{B} is defined as:

\mathbf{B} = \begin{bmatrix} \mathbf{B}_\mathrm{ss} & \mathbf{B}_\mathrm{sg} \\ \mathbf{B}_\mathrm{gs} & \mathbf{B}_\mathrm{gg} \end{bmatrix} = \begin{bmatrix} \rho_1 \quad \dots \quad \rho_m \quad \omega_1 \quad \dots \quad \omega_n \\ \vdots \quad \ddots \quad\quad \vdots \quad\quad \vdots \quad \ddots \quad \vdots \\ \rho_1 \quad \dots \quad \rho_m \quad \omega_1 \quad \dots \quad \omega_n \end{bmatrix}

where \rho_i is the reflectivity of boundary element i and \omega_i=\sigma_{\mathrm{s},i}/\beta_i is the single scattering albedo of volume element i. The single reflection-scattering interaction matrix \mathbf{K} and the steady state path matrix \mathbf{S}_\infty are defined as:

\mathbf{K} = \mathbf{F} \circ \mathbf{B} \quad , \quad \mathbf{S}_\infty = (\mathbf{I}-\mathbf{K})^{-1} \mathbf{F}

where \circ is the Hadamard product and \mathbf{K} combines first interaction probabilities with single reflection-scattering probabilities, and \mathbf{S}_\infty describes reflection-scattering paths. The absorption matrix \mathbf{A} and the multiple reflection-scattering matrix \mathbf{R} are expressed as:

\mathbf{A} = (\mathbf{1} - \mathbf{B})^T \circ \mathbf{S}_\infty \circ (\mathbf{1} - \mathbf{B}) \tag{9.2}

\mathbf{R} = (\mathbf{1} - \mathbf{B})^T \circ \mathbf{S}_\infty \circ \mathbf{B} \tag{9.3}

where \mathbf{1} is a matrix of ones. The first matrix, (\mathbf{1} - \mathbf{B})^T, accounts for the probability of emission, the middle matrix, \mathbf{S}_\infty, accounts for all possible reflection-scattering paths at the steady-state, and the last matrices, \mathbf{1} - \mathbf{B} or \mathbf{B} account for absorption or reflection at the destination.

9.2.3 Application of the Transformed Exchange Factors

Next, based on fundamental radiative energy balances, which state that an energy source is equal to the total radiant power, minus the absorbed power, minus the reflected incident power, and that the emissive power equals the total radiant power, minus the reflected incident power, the matrices \mathbf{C} and \mathbf{D} are defined as:

\mathbf{C} = \mathbf{I} - \mathbf{A}^T - \mathbf{R}^T \quad , \quad \mathbf{D} = \mathbf{I} - \mathbf{R}^T \tag{9.4}

Now, given a vector \mathbf{e} of the emissive powers and a vector \mathbf{q} of the source terms of each element, where only one entry of these can be specified for each element, the mixed-boundary matrix \mathbf{M} can be assembled row by row:

\mathbf{M}[i,:] = \mathbf{C}[i,:] \quad , \quad \mathbf{M}[j,:] = \mathbf{D}[j,:]

using i for known entries of the source vector \mathbf{q}_i and using j for known entries of the emissive power vector \mathbf{e}_j, where i and j are mutually exclusive. Next, the linear system:

\mathbf{M}\mathbf{j} = \mathbf{h} \tag{9.5}

is solved for the total radiant power \mathbf{j} of each element, where \mathbf{h} is a vector of known emissive powers or source terms, depending on which is known, corresponding to the i and j rows of \mathbf{M}.

Knowing the total radiant power of each element, the absorbed incident power \mathbf{g}_\mathrm{a} and the unknown terms of the emissive power vector \mathbf{e} and the source term vector \mathbf{q} can be determined from:

\mathbf{g}_\mathrm{a} = \mathbf{A}^T \mathbf{j} \quad , \quad \mathbf{e}_i = \mathbf{q}_i + \mathbf{g}_{\mathrm{a},i} \quad, \quad \mathbf{q}_j = \mathbf{e}_j - \mathbf{g}_{\mathrm{a},j} \tag{9.6}

where i and j are again mutually exclusive. The reflected-scattered power \mathbf{r} and the total incident power \mathbf{g} can be determined from:

\mathbf{r} = \mathbf{R}^T \mathbf{j} \quad , \quad \mathbf{g} = \mathbf{g}_\mathrm{a} + \mathbf{r} \tag{9.7}

The total intensity of volumes and of Lambertian surfaces can be determined from:

\mathbf{i}_{\mathrm{V},i} = \frac{\mathbf{j}_{\mathrm{V},i}}{4\pi V_i} \quad , \quad \mathbf{i}_{\mathrm{S},i} = \frac{\mathbf{j}_{\mathrm{S},i}}{\pi A_i}

Finally the temperature of each element can be determined from its emissive power using \varepsilon and \kappa and the Stefan-Boltzmann law:

T_{\mathrm{s},i} = \left( \frac{\mathbf{e}_{\mathrm{s},i} }{ \varepsilon_{\mathrm{s},i}\sigma A_{\mathrm{s},i} }\right)^{1/4} \quad , \quad T_{\mathrm{g},i} = \left( \frac{\mathbf{e}_{\mathrm{g},i} }{4 \kappa\sigma V_{\mathrm{g},i} n_i^2 }\right)^{1/4} \tag{9.8}

where \sigma is the Stefan-Boltzmann constant and n_i is the refractive index of the medium of element i.

As mentioned in the introduction, for a medium with varying refractive index, \mathbf{F} should be obtained by solving differential equations for each ray trajectory. While differential equations require a continuously changing refractive index, the proposed transformation uses a discrete mesh of the domain. For this reason, a mean value of the refractive index in each volume element should be used in equation (Eq. 9.8) to calculate the temperature. The ray trajectories in a medium with variable refractive index are found by solving the following vector differential equation for each ray (Born and Wolf 1980):

\frac{d}{dS}\left(n \frac{d \mathbf{x}}{dS}\right) = \nabla n

where \mathbf{x} is the position and S is the distance along a ray path and $$ is the gradient operator. Since the proposed method uses only the points of emission and first interaction, everything regarding the geometric ray trajectories is handled during the ray tracing step.

# The Exchange Factor Transformation {#sec-exchange_fact_trans} ## Introduction Similarly to Hottel's method, outlined in chapter @sec-matrix_methods, the exchange factor transformation of the present work starts from the exchange area matrix. But unlike Hottel's method, which works directly with the exchange area matrix, the proposed transformation uses its non-dimensional counterpart, the *exchange factor* matrix, where each entry of the exchange area matrix has been divided by its equivalent area of emission. Before solving a specific problem, the proposed method apply a matrix transformation to augment the exchange factor matrix with additional information about the steady state of reflection and scattering in the domain. The transformation of the exchange factor matrix into an augmented form is partly inspired by the "F-hat" method developed by Beckman [@Beckman1968], which extends the view factor matrix for a transparent enclosure to include reflections. Further analysis revealed that the 'F-hat' method only includes single reflections and does not account for multiple reflections. The effective reflection-scattering coefficients of multiple reflection-scattering processes on general domains are functions of extinction and geometry, which seemingly made the problem impossible to solve. However, as will be demonstrated in this book, the problem is neither intractable nor impossible to solve. In fact, it has a closed form analytical solution. Given the availability of an analytically derived matrix of exchange factors, generally referred to as the matrix $\mathbf{F}$, the entire solution is analytical. Thanks to Narayanaswamy [@Narayanaswamy2015], who found a general analytical solution to equation (@eq-excfact_ss) with $\overline{t}(S)=1$, this is possible for any convex enclosure composed of polygon surfaces and a transparent medium, giving perfect solutions to multiple reflection systems, which would otherwise be represented by coupled Fredholm equations of the second kind [@Minkowycz1988], the single assumption being that angular distributions of emission and reflection are identical. For problems involving participating media ($\overline{t}(S)<1$), analytical solutions to the multi-dimensional integrals of equations (@eq-excfact_ss)-(@eq-excfact_gg) have yet to be found. In this case, one might employ Monte Carlo methods to obtain near-analytical estimates of $\mathbf{F}$. To ensure that $\mathbf{F}$ satisfies both The First and The Second Law of Thermodynamics, smoothing should be employed, such as for example the method proposed by Daun et al. [@Daun2005], or preferably the method proposed in this work. The solution method proposed in this work is general and, furthermore, it is flexible. In fact, a key strength of exchange factor-based methods for solution of the RTE, is that they allow the formidable monolithic problem of energy transfer by an electromagnetic wave field, to be broken down into two conceptually separate and more manageable parts: The geometry and extinction-based ray propagation part and the system-wide energy interaction part. If $\mathbf{F}$ is obtained through Monte Carlo ray tracing methods, any number of relevant optical phenomena may be modelled, such as refraction and diffraction, refraction being significant in atmospheric research, due to the variable refractive index. This can be achieved by solving differential equations for each ray [@BornWolf1980]. Since the proposed transformation records only ray trajectory endpoints, it is in some sense blind to the phenomena which govern the trajectories, which explains why the matrix transformations of the present work will preserve any modelled ray trajectory properties. For the proposed transformation, since $\mathbf{F}$ refers only to the first interaction, Monte Carlo methods used to obtain $\mathbf{F}$ can be truncated after this initial interaction, making them computationally less demanding than basic Monte Carlo implementations that must track multiple reflection-scattering events [@HowellDaun2021]. This will be increasingly true for higher single reflection-scattering coefficients due to the increased length of the ray paths. Instead, the proposed transformation performs the multiple reflection-scattering ray tracing analytically using solution of linear systems. This approach offers key advantages: (1) It significantly improves computational efficiency compared to basic ray tracing methods for multiple reflection-scattering, (2) It provides analytical ray tracing accuracy with quantifiable bounds on the necessary numerical precision, (3) The precision requirements are expressed as a function of both the geometry and extinction (described by $\mathbf{F}$) and the reflection and scattering properties of each element. This enables problem-specific numerical precision tuning which is readily implementable in modern numerical programming languages [@IEEEfloat2019]. ## Methodology The methodology is split into three parts: First, the exchange factor matrix is formally defined, next, the exchange factor transformation is defined, and lastly, the transformed exchange factors are applied to solution of the RTE. ### Definition of the Exchange Factor Matrix The proposed transformation requires as an input a general description of the radiative connectivity of the domain in the form of the exchange factor matrix $\mathbf{F}$. The exchange factor matrix $\mathbf{F}$ may be obtained analytically or numerically. All of the subsequent operations of the proposed transformation are analytical. The $(m+n)\times(m+n)$ row-stochastic blocked first interaction exchange factor matrix $\mathbf{F}$ for a system of $m$ bounding surfaces and $n$ gas volume elements is defined as: $$ \mathbf{F} = \begin{bmatrix} \mathbf{F}_\mathrm{ss} & \mathbf{F}_\mathrm{sg} \\ \mathbf{F}_\mathrm{gs} & \mathbf{F}_\mathrm{gg} \end{bmatrix} $$ where $\mathbf{F}$ is row stochastic to ensure that The First Law of Thermodynamics is satisfied and the entries of $\mathbf{F}$ are the fractions of energy emitted by the emitter of row $i$ which has its first interaction with the element of column $j$. The entries of $\mathbf{F}$ can be defined from equations (@eq-excfact_ss)-(@eq-excfact_gg), by normalizing with the equivalent total capacity for emission and reflection-scattering of the emitter zone (and using the extinction coefficient $\beta$ in place of the absorption coefficient $\kappa$): $$ F_{s_j s_k} = \frac{\overline{s_j s_k}}{(\rho_j+\varepsilon_j)A_j}, \quad F_{s_k g_\gamma} = \frac{\overline{s_k g_\gamma}}{(\rho_k+\varepsilon_k)A_k}, \quad F_{g_\gamma s_k} = \frac{\overline{g_\gamma s_k}}{4\beta V_\gamma}, \quad F_{g_\gamma g_{\gamma*}} = \frac{\overline{g_\gamma g_{\gamma*}}}{4\beta V_\gamma}. $$ where $\rho+\varepsilon=1$ where $\rho$ is the single reflectivity and $\varepsilon$ is the emissivity and for volumes of participating media, $\kappa+\sigma_s=\beta$ where $\sigma_s$ is the single scattering coefficient. Additionally, $\mathbf{F}$ satisfies The Second Law of Thermodynamics (reciprocity) [@Daun2005]: $$ \mathbf{E} \mathbf{F} = \mathbf{F}^T \mathbf{E} $$ {#eq-reciprocity} where $\mathbf{E}$ is a diagonal matrix of equivalent emission areas with $\mathbf{E}_{ii}=(\rho_i+\varepsilon_i) A_i$ for $i\leq m$ and $\mathbf{E}_{ii}=4 \beta_i V_i$ for $i>m$, where the full capacity to emit and reflect-scatter is used. Furthermore, in the general case, $\mathbf{F}$ will have some zero entries, i.e., physical elements without direct view of each other, meaning $\mathbf{F}$ is not merely a positive matrix but a non-negative matrix. Therefore, for the proofs of this book to guarantee convergence, $\mathbf{F}$ must be *irreducible* [@HornJohnson2013], where irreducibility corresponds to the physical requirement that the domain contains no isolated regions, that is, for any two elements $i$ and $j$, there exists a sequence of direct radiation exchanges $i \to k_1 \to k_2 \to \ldots \to j$ connecting them. This requirement can always be satisfied by treating isolated regions as separate radiative transfer problems. For the present transformation, $\mathbf{F}$ should be derived exactly as if dealing with pure absorption, but always using the total capacity to emit and reflect-scatter instead of the absorption coefficient. Note that extinction $\beta$ does not need to be uniform across the domain, but should be uniform for each element. In the case of non-uniformity, $\mathbf{F}$ must be obtained using the same non-uniform distribution of $\beta$. Then after $\mathbf{F}$ has been obtained, any part of the capacity to absorb may be converted into a capacity for reflection-scattering, meaning a single $\mathbf{F}$ covers all combinations of $\kappa+\sigma_s$ which sum to $\beta$. ### The Exchange Factor Transformation The column constant $(m+n)\times(m+n)$ blocked single reflection-scattering matrix $\mathbf{B}$ is defined as: $$ \mathbf{B} = \begin{bmatrix} \mathbf{B}_\mathrm{ss} & \mathbf{B}_\mathrm{sg} \\ \mathbf{B}_\mathrm{gs} & \mathbf{B}_\mathrm{gg} \end{bmatrix} = \begin{bmatrix} \rho_1 \quad \dots \quad \rho_m \quad \omega_1 \quad \dots \quad \omega_n \\ \vdots \quad \ddots \quad\quad \vdots \quad\quad \vdots \quad \ddots \quad \vdots \\ \rho_1 \quad \dots \quad \rho_m \quad \omega_1 \quad \dots \quad \omega_n \end{bmatrix} $$ where $\rho_i$ is the reflectivity of boundary element $i$ and $\omega_i=\sigma_{\mathrm{s},i}/\beta_i$ is the single scattering albedo of volume element $i$. The single reflection-scattering interaction matrix $\mathbf{K}$ and the steady state path matrix $\mathbf{S}_\infty$ are defined as: $$ \mathbf{K} = \mathbf{F} \circ \mathbf{B} \quad , \quad \mathbf{S}_\infty = (\mathbf{I}-\mathbf{K})^{-1} \mathbf{F} $$ where $\circ$ is the Hadamard product and $\mathbf{K}$ combines first interaction probabilities with single reflection-scattering probabilities, and $\mathbf{S}_\infty$ describes reflection-scattering paths. The absorption matrix $\mathbf{A}$ and the multiple reflection-scattering matrix $\mathbf{R}$ are expressed as: $$ \mathbf{A} = (\mathbf{1} - \mathbf{B})^T \circ \mathbf{S}_\infty \circ (\mathbf{1} - \mathbf{B}) $$ {#eq-A} $$ \mathbf{R} = (\mathbf{1} - \mathbf{B})^T \circ \mathbf{S}_\infty \circ \mathbf{B} $$ {#eq-R} where $\mathbf{1}$ is a matrix of ones. The first matrix, $(\mathbf{1} - \mathbf{B})^T$, accounts for the probability of emission, the middle matrix, $\mathbf{S}_\infty$, accounts for all possible reflection-scattering paths at the steady-state, and the last matrices, $\mathbf{1} - \mathbf{B}$ or $\mathbf{B}$ account for absorption or reflection at the destination. ### Application of the Transformed Exchange Factors Next, based on fundamental radiative energy balances, which state that an energy source is equal to the total radiant power, minus the absorbed power, minus the reflected incident power, and that the emissive power equals the total radiant power, minus the reflected incident power, the matrices $\mathbf{C}$ and $\mathbf{D}$ are defined as: $$ \mathbf{C} = \mathbf{I} - \mathbf{A}^T - \mathbf{R}^T \quad , \quad \mathbf{D} = \mathbf{I} - \mathbf{R}^T $$ {#eq-C_and_D} Now, given a vector $\mathbf{e}$ of the emissive powers and a vector $\mathbf{q}$ of the source terms of each element, where only one entry of these can be specified for each element, the mixed-boundary matrix $\mathbf{M}$ can be assembled row by row: $$ \mathbf{M}[i,:] = \mathbf{C}[i,:] \quad , \quad \mathbf{M}[j,:] = \mathbf{D}[j,:] $$ using $i$ for known entries of the source vector $\mathbf{q}_i$ and using $j$ for known entries of the emissive power vector $\mathbf{e}_j$, where $i$ and $j$ are mutually exclusive. Next, the linear system: $$ \mathbf{M}\mathbf{j} = \mathbf{h} $$ {#eq-solve_for_j} is solved for the total radiant power $\mathbf{j}$ of each element, where $\mathbf{h}$ is a vector of known emissive powers or source terms, depending on which is known, corresponding to the $i$ and $j$ rows of $\mathbf{M}$. Knowing the total radiant power of each element, the absorbed incident power $\mathbf{g}_\mathrm{a}$ and the unknown terms of the emissive power vector $\mathbf{e}$ and the source term vector $\mathbf{q}$ can be determined from: $$ \mathbf{g}_\mathrm{a} = \mathbf{A}^T \mathbf{j} \quad , \quad \mathbf{e}_i = \mathbf{q}_i + \mathbf{g}_{\mathrm{a},i} \quad, \quad \mathbf{q}_j = \mathbf{e}_j - \mathbf{g}_{\mathrm{a},j} $$ {#eq-emissive_source} where $i$ and $j$ are again mutually exclusive. The reflected-scattered power $\mathbf{r}$ and the total incident power $\mathbf{g}$ can be determined from: $$ \mathbf{r} = \mathbf{R}^T \mathbf{j} \quad , \quad \mathbf{g} = \mathbf{g}_\mathrm{a} + \mathbf{r} $$ {#eq-r_and_g} The total intensity of volumes and of Lambertian surfaces can be determined from: $$ \mathbf{i}_{\mathrm{V},i} = \frac{\mathbf{j}_{\mathrm{V},i}}{4\pi V_i} \quad , \quad \mathbf{i}_{\mathrm{S},i} = \frac{\mathbf{j}_{\mathrm{S},i}}{\pi A_i} $$ Finally the temperature of each element can be determined from its emissive power using $\varepsilon$ and $\kappa$ and the Stefan-Boltzmann law: $$ T_{\mathrm{s},i} = \left( \frac{\mathbf{e}_{\mathrm{s},i} }{ \varepsilon_{\mathrm{s},i}\sigma A_{\mathrm{s},i} }\right)^{1/4} \quad , \quad T_{\mathrm{g},i} = \left( \frac{\mathbf{e}_{\mathrm{g},i} }{4 \kappa\sigma V_{\mathrm{g},i} n_i^2 }\right)^{1/4} $$ {#eq-temperature} where $\sigma$ is the Stefan-Boltzmann constant and $n_i$ is the refractive index of the medium of element $i$. As mentioned in the introduction, for a medium with varying refractive index, $\mathbf{F}$ should be obtained by solving differential equations for each ray trajectory. While differential equations require a continuously changing refractive index, the proposed transformation uses a discrete mesh of the domain. For this reason, a mean value of the refractive index in each volume element should be used in equation (@eq-temperature) to calculate the temperature. The ray trajectories in a medium with variable refractive index are found by solving the following vector differential equation for each ray [@BornWolf1980]: $$ \frac{d}{dS}\left(n \frac{d \mathbf{x}}{dS}\right) = \nabla n $$ where $\mathbf{x}$ is the position and $S$ is the distance along a ray path and $\nabla $ is the gradient operator. Since the proposed method uses only the points of emission and first interaction, everything regarding the geometric ray trajectories is handled during the ray tracing step.