# Radiosity The radiosity integral equations describe the energy transport via radiation between surfaces. This formulation not only describes direct energy transfer through absorption, but also incorporates transfer through scattered light and emission and absorption of infrared thermal energy (self heating). The most important aspect within the radiosity equations is the viewfactors. Viewfactors describe for two facets, how much of a fraction of emitted light one facet receives from the other and they are purely based on the geometry between two facets. In our code we use these equations to calculate the absolute energy value that each facet receives. We have followed the basic implementation of [Potter et al. (2023)](../literature) in our code, but their Fast Radiosity Algorithm is not yet implemented. Furthermore, we have adapted the equations to work with a incidence angle dependent albedo. To verify our implementation of the viewfactor calculation and energy calculation, we compared our results for a spherical crater with no sub-surface heat transport to the analytical solution described by [Ingersoll et al. (1992)](../literature) as can be seen in the image below. Due to the finite number of facets, the fluctuations of the values in the shadowed region in the crater are to be expected. The crater has been generated with the publicly available code by [Potter et al. 2023](../literature).
drawing
## Integral equations **Introduction** The radiosity integral equations are the mathematical description of the energy transfer via radiation based only on geometric assumptions. The basic form of the equations is $$\mathbf{B} = \mathbf{E} + \text{diag}(\mathbf{\rho}) \mathbf{F} \mathbf{B}$$ where $\mathbf{B}$, $\mathbf{E}$ and $\mathbf{\rho}$ are vectors and $\mathbf{F}$ is the viewfactor matrix [(Potter et al. 2023)](../literature). **Physics - Lambertian albedo approach** The [surface energy balance](boundary_conditions) is used to describe the temperature of (planetary) surfaces. When incorporating the radiosity equations, the incoming energy/power is not described by a single term $Q_{\text{sol}}$, but by three that include not only direct absorption of incoming light, but also reradiated heat $Q_{\text{ir}}$ and scattered light $Q_{\text{refl}}$ from surfaces within the field of view of the current facet. With these changes, the surface energy balance reads $$(1 - \alpha) (Q_{\text{sol}} + Q_{\text{refl}}) + \lambda \mathbf{n} \cdot \nabla T + \epsilon Q_{\text{ir}} = \epsilon \sigma T^4 .$$ The direct irradiation from the sun is calculated by the distance to the sun, the angle to the surface and whether it is shadowed by other parts of the body or not, but the other two contributions are described by the radiosity integral equations via $$Q_{\text{refl}}(\mathbf{x}) = \int_S F(\mathbf{x}, \mathbf{y}) \alpha(\mathbf{y}) \left[Q_{\text{sol}}(\mathbf{y}) + Q_{\text{refl}}(\mathbf{y})\right] dA(\mathbf{y})$$ $$Q_{\text{ir}}(\mathbf{x}) = \int_S F(\mathbf{x}, \mathbf{y}) \left[\epsilon (\mathbf{y}) \sigma T(\mathbf{y})^4 + (1 - \epsilon (\mathbf{y})) Q_{\text{ir}}(\mathbf{y})\right] dA(\mathbf{y})$$ and it can be shown that both of these equations are of the form of $$\mathbf{B} = \mathbf{E} + \text{diag}(\mathbf{\rho}) \int_S \mathbf{F} \mathbf{B} dA$$ [(Potter et al. 2023)](../literature). Discretizing these integrals into sums and using some simple transformations, [Potter et al. (2023)](../literature) show that the absolute absorbed energy/power can be calculated with an update at each time step ($n \rightarrow n+1$) of the thermophysical model through $$Q_{\text{rad}}^{n}(\mathbf{x}_i) = \epsilon (\mathbf{x}_i) \sigma T(\mathbf{x}_i)^4, \quad 1 \leq i \leq N$$ $$\mathbf{Q}_{\text{refl}}^{n+1} = \mathbf{F} \text{diag}(\mathbf{\alpha}) \left[\mathbf{Q}_{\text{sol}}^{n+1} + \mathbf{Q}_{\text{refl}}^{n}\right]$$ $$\mathbf{Q}_{\text{ir}}^{n+1} = \mathbf{F} \left[\mathbf{Q}_{\text{rad}}^{n} + \left(\mathbf{I} - \text{diag}(\mathbf{\epsilon})\right) \mathbf{Q}_{\text{ir}}^{n} \right]$$ $$\mathbf{Q}_{\text{abs}}^{n+1} = \left(\mathbf{I} - \text{diag}(\mathbf{\alpha})\right) \left(\mathbf{Q}_{\text{sol}}^{n+1} + \mathbf{Q}_{\text{refl}}^{n+1} \right) + \text{diag}(\mathbf{\epsilon}) \mathbf{Q}_{\text{ir}}^{n+1} .$$ In this approach, the scattering of light lags behind in the thermal model, as each multiplication of the viewfactor matrix corresponds to one order of light scattering. This means that the model only "sees" the effects of the second scattering from timestep $n$ at timestep $n+1$, the third scattering at timestep $n+2$ and so on. **Incidence angle dependent albedo** The placement of the albedo vector $\mathbf{\alpha}$ in these equations already implies a lambertian treatment of the albedo with no dependency on the incidence angle. This can be remedied while keeping the equations mostly the same through a slight simplification. If we assume that the variation of angles between the sun and a facet and between any two facets is small, we can only use one characteristic angle and do not need to treat all possible angles between the two areas. This allows forgoing any integration over the angles and simply using a fixed albedo-matrix $\alpha(\mathbf{x}_i, \mathbf{y}_j)$ for the angles between the facets and an albedo-vector $\alpha_{\text{sol}}(\mathbf{x}_i)$ for the angles between the sun and the facets. With these changes, any of the equations from the previous section that use the albedo have to be slightly altered, to incorporate the new albedo-matrix: $$Q_{\text{refl, total}}^{n+1}(\mathbf{x}_i) = \sum_{j=1}^{N} \left(1 - \alpha(\mathbf{x}_i, \mathbf{y}_j) \right) \cdot \mathbf{F}(\mathbf{x}_i, \mathbf{y}_j) \cdot \left[\alpha_{\text{sol}}(\mathbf{y}_j) Q_{\text{sol}}^{n+1}(\mathbf{y}_j) + \alpha(\mathbf{x}_i, \mathbf{y}_j) Q_{\text{refl}}^{n}(\mathbf{y}_j)\right]$$ $$Q_{\text{abs}}^{n+1}(\mathbf{x}_i) = \left(1 - \alpha(\mathbf{x}_i)\right) Q_{\text{sol}}^{n+1} + Q_{\text{refl}}^{n+1}(\mathbf{x}_i) + \epsilon(\mathbf{x}_i) Q_{\text{ir}}^{n+1}(\mathbf{x}_i) .$$ ## Viewfactors **Introduction** The viewfactor $F(\mathbf{x}_i, \mathbf{x}_j)$ is a geometrical value that is a measure for how much two facets "see" each other. In other words, if one facet emits light isotropically into one hemisphere, the viewfactor gives the fraction of that light that reaches the other facet. In a closed space with an inwards emitting surface (ideal black body), all viewfactors would add up to exactly 1. This is not to be confused with the viewfactor matrix $\mathbf{F}$, which is often employed in the actual calculations of the radiosity integral equations. For $n$ facets, this matrix is of dimension $n\times n$ and contains the viewfactors between each facet scaled with their areas. So for facet $i$ and $j$ the two entries are $$ \mathbf{F}_{i,j} = F(\mathbf{x}_i, \mathbf{x}_j) \cdot A_j$$ $$ \mathbf{F}_{j,i} = F(\mathbf{x}_i, \mathbf{x}_j) \cdot A_i$$ with the areas of the facets $A$. So while the viewfactor is symmetrical between any two facets, the respective entries in the viewfactor matrix are usually asymmetrical. **Calculation** In the most general formulation, the viewfactor between to surfaces can be expressed through a double area integral $$F(\mathbf{x}_i, \mathbf{x}_j) = \frac{1}{A_i A_j} \int_{A_j} \int_{A_i} \frac{\cos(\theta_i) \cos(\theta_j)}{\pi d^2} dA_i dA_j$$ with the angle between the surface normal and the vector connecting both surfaces $\theta$ and the distance between the two surfaces $d$ [(Siegel and Howell 1971)](../literature). If the surfaces are obstructed through a third surface, then the viewfactor would be zero. Viewfactor calculation in science and computer graphics has been done for a long time and many analytical solutions for different cases have been found and can be viewed in multiple published catalogs. But when assuming flat polygonal surfaces with arbitrary angles between each other, as it is needed for applications in planetary sciences, close form analytical solutions to these integrals only exist for a few specific cases. There are many different ways to calculate the viewfactors given an arbitrary planar geometry. The most simple is to just take the angle between the centroids (geometric centers) of each facet and approximate the integral via the midpoint rule to receive $$ F(\mathbf{x}_i, \mathbf{x}_j) = \frac{\left[\mathbf{n}_i \cdot (\mathbf{x}_j - \mathbf{x}_i) \right] \cdot \left[\mathbf{n}_j \cdot (\mathbf{x}_i - \mathbf{x}_j) \right]}{\pi \left|\mathbf{x}_i - \mathbf{x}_j \right|^4}$$ where the cosine vector equality $\cos(\theta_{ab}) = (\mathbf{a} \cdot \mathbf{b})/(|\mathbf{a}| |\mathbf{b}|)$ was used [(Potter et al. 2023, Lagerros 1997)](../literature). This is a very simplified calculation which works well for surfaces that are far away, but diverges more strongly from the expected solution for surfaces that are close. This can be seen when using this calculation on a cube with inward facing normals, as the viewfactors do not add up to one, but are instead greater than one. But for the usual application in planetary science, this approach is often adequately good. A major strength of this approach is that it is the most computationally inexpensive, as only simple arithmetic and two dot products have to be evaluated. Thus it sees far spread application in planetary science [(Potter et al. 2023, King et al. 2020, Gutiérrez et al. 2001, Lagerros 1997)](../literature). A more sophisticated, but already much more computationally expensive approach, is the one described in [Mazumder and Ravishankar (2012)](../literature). They transform the double area integral into a double line integral around the boundary of the surface with Stoke's theorem. Under the assumptions of two planar polygons with an arbitrary number of vertices and their contours $\Gamma$ they receive $$ \mathbf{F}_{ij} = \oint_{\Gamma_j} \oint_{\Gamma_i} \ln(S) d\mathbf{s}_i d\mathbf{s}_j$$ where $d\mathbf{s}$ is the differential line element on an edge of the polygon and $S$ is the magnitude of the vector connecting the line elements of both polygons, so the distance of the line elements. This distance can be expressed in the unitary coordinates $\lambda \in [0, 1]$ via $$\mathbf{S} = - \mathbf{s}_i + \overrightarrow{I_m J_n} + \mathbf{s}_j = - \lambda_i \mathbf{I_{m,m+1}} + \overrightarrow{I_m J_n} + \lambda_j \mathbf{J_{n,n+1}}$$ where $I_n$ and $J_m$ are specific vertices of their respective polygon and $\mathbf{I_{m,m+1}} = \overrightarrow{I_m I_{m+1}}$. This allows for fixed unitary boundaries on the line integrals, which are evaluated piecewise. The resulting integrals are given by $$ \mathbf{F}_{ij} = \frac{1}{4 \pi A_i} \sum_{n=1}^N \sum_{m=1}^M \int_0^1 \int_0^1 \ln\left( \lambda_i^2 |\mathbf{I}_{m,m+1}|^2 + \lambda_j^2 |\mathbf{J}_{n,n+1}|^2 - 2 \lambda_i \overrightarrow{I_m J_n} \cdot \mathbf{I}_{m,m+1} + 2 \lambda_j \overrightarrow{I_m J_n} \cdot \mathbf{J}_{n,n+1} - 2 \lambda_j \lambda_i \mathbf{J}_{n,n+1} \cdot \mathbf{I}_{m,m+1} + |\overrightarrow{I_m J_n}|^2 \right) \mathbf{J}_{n,n+1} \cdot \mathbf{I}_{m,m+1} d\lambda_i d\lambda_j$$ These integrals have no closed form analytical solution and need to be evaluated numerically. Lastly, the authors give the analytical contribution of a shared edge between polygons, as they would produce singular results if treated numerically as $$\Delta \mathbf{F}_{ij}|_{\text{shared}} = |\mathbf{I}_{\text{shared}}|^2 \left( \frac{3}{2} - \frac{1}{2} \ln \left\{ |\mathbf{I}_{\text{shared}}|^2 \right\} \right) .$$ This formulation gives reproduces analytical viewfactors up to six significant digits when using a 5-point Gaussian-Legendre quadrature scheme to evaluate the integrals [Mazumder and Ravishankar (2012)](../literature). Lastly, we want to briefly discuss a third method of viewfactor calculation which uses Monte-Carlo method. This is the most widely applicable one, as it can be used in many forms and complexities from simply calculating the viewfactors [(Emery et al. 1988)](../literature) up to completely solving the radiative energy balance equations and accounting for light scattering and absorption [(Ertürk and Howell 2017)](Litertaure). It inherently accounts for partial shadowing effects, where a surface is only partially occluded by other pieces of geometry. The problem is, that this method is prohibitatively computationally expensive. Not only are millions of rays necessary to reduce the statistical noise, but also further smoothing can often times be necessary [(Mazumder and Ravishankar 2012)](../literature). Thus we have decided not to use this method.