Accuracy of view factor calculations for digital terrain models of comets and asteroids

Context. Detailed shape and topographic models coupled with sophisticated thermal physics are critical elements to proper charac- terization of surfaces of small bodies in our solar system. Calculations of self-heating e ﬀ ects are especially important in the context of thermal evolution of non-convex surfaces, including craters, cracks, or openings between “rocks”. Aims. Our aim is to provide quantitative comparisons of multiple numerical methods for computing view factors for concave geome- tries and provide a more rigorous criteria for the validity of their application. Methods. We contrasted ﬁve methods of estimating the view factors. First, we studied speciﬁc geometries, including shared-edge facets for a reduced two-facet problem. Then, we applied these methods to the shape model of 67P / Churyumov-Gerasimenko. Nev-ertheless, the presented results are general and could be extended to shape models of other bodies as well. Results. The close loop transformation of the double area integration method for evaluating view factors of nearby or shared-edge facets is the most accurate, although computationally expensive. Two methods of facet subdivision we evaluate in this work provide reasonably accurate results for modest facet subdivision numbers, however, may result in a degraded performance for speciﬁc facet geometries. Increasing the number of subdivisions improves their accuracy, but also increases their computational burden. In practical applications, a trade-o ﬀ between accuracy and computational speed has to be found, therefore, we propose a combined method based on a simple metric that incorporates a conditional application of various methods and an adaptive number of subdivisions. In our study case of a pit on 67P / CG, this method can reach average accuracy of 2-3% while being about an order of magnitude faster than the (most accurate) line integral method.


Introduction
The main energy input into a cometary nucleus is solar irradiation, which is then redistributed via energy sink processes such as sublimation of ices, conduction into the subsurface, or infrared surface emissions. Thermal re-emission from a surface element may be absorbed at another distinct surface element, becoming its secondary source of energy. This is the so-called selfheating effect, which is most pronounced for certain topographic environments such as pits, craters, cracks, and cliffs (Colwell et al. 1990;Lagerros 1997;Ivanova & Shulman 2006). In this work, we ignore the reflection of visible light as a secondary source of energy, as considered for example in Rozitis & Green (2011). The importance of the self-heating effect is well recognized and nowadays many thermophysical models incorporate these calculations in the context of high-resolution shape models of comets (e.g., 67P/Churyumov-Gerasimenko) and asteroids (e.g., Ryugu, Bennu, etc) (Gutiérrez, P. J. et al. 2001;Rozitis & Green 2011;Hu et al. 2017;Shi, X. et al. 2016;Pelivan 2018;Hamm et al. 2018). We should add that the complexity of calculations of the net radiation exchange among surface elements (facets) can be made less complex by assuming that each surface rezac@mps.mpg.de is diffusive and has homogeneous emissive power (Lambertian surface). This work is also done within this framework.
The mathematical formalism for including the self-heating effect requires calculations of the so-called view factor(s). A vast body of literature in the field of computer graphics and engineering simulations deals with the complexities of these calculations (Siegel & Howell 2002;Modest 2013). There is still ongoing research into improving these algorithms (Francisco et al. 2014;Kramer et al. 2015;Narayanaswamy 2015). By definition, the view factor represents a fraction of radiation from a surface element, i, directly reaching a facet, j. Its differential form is be expressed as In this form and under stated assumptions a view factor is purely a geometric quantity independent of surface properties and temperature. This expression relates cosines of angles (φ i j , φ ji ) between the normal direction of each facet and the ray connecting the two differential areas dA i and dA j , on facet i and j, respectively (Fig. 1), and the square of the distance between the differential areas, d i j . The binary flag, δ i j , is a convenient parameter indicating whether there is a direct visibility between Article number, page 1 of 9 Article published by EDP Sciences, to be cited as https://doi.org/10.1051/0004-6361/202038462 A&A proofs: manuscript no. selfheat-final facets. View factors possess an important reciprocity relation, dA i dF i j = dA j dF ji . Several methods can be used to solve the integration over Eq. 1 to find self-heating of a finite area. In this work, we investigate the accuracy and numerical efficiency of these methods when applied in the context of concave geometries (e.g., pits, cracks, and fractures) for a high-resolution nucleus shape model. From a numerical point of view it is challenging to properly evaluate view factors for neighboring facets sharing an edge and adjacent facets inclined at an angle, but not necessarily sharing an edge.
The paper is structured as follows. In Section 2, we provide a detailed description of the five methods and three geometry cases we study in this work. In Section 3, we show the simulation results of comparison between those methods for various geometry conditions. We also evaluate performance of these algorithm for view factor calculation in the context of the surface of 67P/CG. In Section 4, we describe and summarize the main lessons we learn from this work and provide our suggestions for view factor evaluation calculation for comets and asteroids.

Method description
When considering finite areas, A i and A j , the view factor expression involves a double integration working with an assumption that the two surfaces are isothermal, black, and diffusively emitting (Siegel & Howell 2002), that is, In certain cases, usually for simple geometries, analytical solutions exist for Eq. 2; a comprehensive list is provided in Siegel & Howell (2002). The most simple form, which we label as a method (M 1 ), can be expressed as (Ozisik 1985;Siegel & Howell 2002;Modest 2013) This straight forward approach is often applied in thermophysical studies of comets and asteroids (Davidsson & Rickman 2014;Delbo et al. 2015). It is also acknowledged that this approach is only valid when A j d i j , and that it can reach a singularity as d i j approaches zero.
In search of a simple analytical solution for small d i j cases, we found an expression in the documentation of the Abaqus® commercial package (Abaqus 2020), which aims to address this issue, as follows: We designate this method as M 2 . The limit of Eq. 4 is Eq. 3 when d i j is large. Generally, the Monte Carlo method is convenient for computing the four-dimensional integral in Eq. 2 (Modest 1978;Yarbrough & Lee 1986;Modest 2003;Howell et al. 2010;Maltby & Burns 1991;Baranoski et al. 2001;Vujičić et al. 2006;Mirhosseini & Saboonchi 2011). However, it may be very computationally expensive for high-resolution digital terrain models (DTM) of comets, which makes other methods more practical for the view factor evaluation.
In particular, the subdivision approach which splits facets into smaller ones (Rozitis & Green 2011; Davidsson & Rickman 2014;Hu et al. 2017;Pelivan 2018;Hamm et al. 2018), offers a reasonable way to fulfill the crucial condition for Eq. 3 and avoid the risk of grossly overestimating view factors. The strategy is to divide facets into smaller triangles A im or/and A jn and apply Eq. 3 to determine the view factors between the smaller facets. The expression of F i j is written as where N di and N d j are the number of subdivisions. This method, adopts a linear subdivision 1 algorithm to split the triangular facets, labeled as M 3 -Linear. This works well for nonadjacent facets, but may prove to be inaccurate when applied to touching (shared-edge) facets. In practice, the approach to the singularity may still be avoided because the contribution of the border region to F i j is scaled down by the fractional area (which acts as a weight in the sum), however, the cost in accuracy has not been yet quantified for typical applications to cometary surface models. The derivative approach of the facet subdivision method, which we call M 3 -R2, is also analyzed in this work. The goal is to avoid the explicit numerical calculation of splitting facets into smaller triangular elements. Instead, it relies on an algorithm efficiently generating an even distribution of points over an arbitrary shaped triangle using a quasi-random sequence (Roberts 2018) (see appendix A for the R2 algorithm description). Our motivation for testing this approach is twofold. First, the M 3 -Linear can subdivide triangles only in multiples of four by definition, which induces limits for optimization as to the number of subtriangles. In addition, the actual (linear) subdivision algorithm (Warren & Weimer 2001) may be computationally costly for repeated calculations as required for applications when the shape model itself is evolving with time (Zhao et al. 2020).
Another approach to solve Eq. 2, first proposed by Sparrow (1963), is to transform the area integrals into two contour integrals via the Stoke's theorem, which could be expressed as (Modest 2003) where Γ i and Γ j represent the contours bounding surface i and j, while s i and s j are differential length vectors along contours. Mazumder & Ravishankar (2012) applied the method for polygon planar surfaces in the form of discrete summations written as where N and M are vertex numbers of facet i and j; for triangular facets we have N = M = 3. The indices m and n indicate the m th and n th vertex of a facet, following anti-clockwise sequence. The expression i m i m+1 is the vector from the m th vertex to m+1 th vertex of the facet i, when m = M; m+1 represents the first vertex. The quantity λ represents the ratio of the length between the starting vertex of the edge to the tail of the differential length vector ds to the length of the edge. In the above expression we note that the forth term in the bracket has a positive sign. There seems to be a misprint in the original reference indicating a negative sign. A Gauss-Legendre quadrature schemes with ten points is employed to solve Eq. 7 (Mazumder & Ravishankar 2012).
In the case of faces sharing an edge, the contribution from the shared edges to Eq. 7 is in the form of (Ambirajan & Venkateshan 1993) where i s is the vector of the shared-edge. This method, labeled as M 4 in this work, is the most accurate in the case of nearby and/or adjacent facets, and it serves as a baseline result.
In the next section, we present results of inter-comparisons of the five different methods expressed in Eq. 3 (M 1 ); Eq. 4 (M 2 ); Eq. 5 (M 3 -Linear); the experimental M 3 -R2 method; and the two equations (Eq. 7 with 8), which we label as M 4 approach. The main goal is to learn about the numerical accuracy and efficiency of these methods for shared-edge and nearby facets for viewing factor computations.

Results
First we investigate the accuracy of methods described above for three cases of simplified geometry of a two facet system as shown in Fig. 2. These reduced complexity numerical experiments aim to study the accuracy of each method from a different point of view: as a function of distance, the angle between facets, and the number of subdivisions for a facet. In these numerical experiments we define a simple metric of "closeness", which combines the relative distance between facets as well as their areas, as This parameter can serve as a good criterion for the validity of Eq. 3, and also can be reliably used to optimize subdivision methods in terms of the number of required samples (see Appendix B for example). The number of subdivisions (N di = N d j ) are the key parameters for the two M 3 methods, and in the following cases we set this number to 256. In the last numerical test we investigate accuracy as a function N di ,and N d j to justify this choice. As already noted, the M 4 method is taken as a reference value because of its high accuracy, especially for shared-edge and close-by facets (Mazumder & Ravishankar 2012). In Fig. 3, a parallel facets case is considered. For both panels of the figure, the abscissa shows the N s metric (Eq. 9), while the ordinate indicates the view factor F 21 (top panel), and relative differences with respect to M 4 in the bottom panel. In this scenario, the N s only varies due to changes in the distance between facets, while areas of the two triangles remain constant. We can immediately see (bottom panel) that the M 1 method, as expected, is the least accurate for close distances (small N s ), however, this approach improves in accuracy as we increase the distance between the facets. The M 2 method demonstrates improved accuracy in comparison with M 1 in this setup. The subdivision techniques are superior in accuracy to both M 1 and M 2 for the entire range of N s values. Even at the smallest N s both M 3 approaches reach a sub-percent difference relative to the M 4 reference. We should point out that this very high accuracy is contingent on the number of subdivisions, which is taken in this work as N di = N d j = 256. This means that the accuracy of M 3 methods can be further improved or degraded by increasing/decreasing N di and N d j , respectively. The improvements in the view factor accuracy come at a computational cost, which are discussed later in this work in connection with a more practical example. Figure 4 also has the same two panels as the previous figure. However, these results correspond to the case in which one of the facets is inclined at 45 • as illustrated in panel B) in Fig. 2. Naturally, the absolute value of the view factor is reduced compared to the parallel facet example. Here, the relative accuracy between M 1 and M 2 does not show any difference, and they only reach the 1% mark of accuracy (with regard to M 4 ) for large N s values of about 500-600. As in the previous case study the two subdivision methods show an accuracy of < 1% for the entire range of N s values. We can also observe a discontinuity in the relative accuracy curve for the M 3 -R2 method (also visible in the previous example). This feature, manifesting as a rapid increase in accuracy over a narrow range of N s values, is purely a numerical artifact of the method. It is related to the mapping of the point via the parallelogram algorithm (see appendix A). However, its dependence on N s is not trivial. From these and other tests we see that it depends on the exact shape of the two facets, as well as the angle between them, and the number of sampling points applied. For instance, increasing the N di and N d j over 300 seems to completely eliminate this feature for this test case.
Fig 5 shows calculations of the view factor F 21 as a function of the opening angle between the two facets sharing an edge (panel C) in Fig. 2). The top panel shows the absolute value of F 21 , while the relative deviation with regards to the M 4 method is plotted in the bottom panel. This geometrical setup is challenging for all methods, but in particular for M 1 and M 2 for small opening angles. In addition, both M 1 and M 2 methods do not reach better than 10-20% accuracy even for large angles, although the M 1 method is not monotonic in its behavior. On the other hand, the two M 3 methods show a trend of increasing accuracy for angles larger than 25 • until the value of about 170 • , where there is a small spike of worsening. The curve for M 3 -R2 stays more flat with an increasing angle going from about 8% (at 25 • ) to 5% at 170 • . The M 3 -Linear method starts at ∼10% accuracy reaching a value 1-2% for large opening angles. This result is important especially for modeling cracks, pits, and other concave surface structures. Nevertheless, a physical shape model of a comet or an asteroid is not expected to have two touching facets with an opening angle smaller than perhaps 30-40 • , however, there are likely to be many cases of large opening angles (> 120 • ).  Fig. 2). The different curves correspond to different methods as labeled. (Bottom) Relative deviation from M 4 . As in the previous case, even in this case the M 3 -Linear and M 3 -R2 methods remain accurate within 1% even for very low N s values. However, the M 3 -R2 method does not asymptotically approach the reference value with increasing N s in this setup. The M 2 and M 1 agreement shown in this figure is just a coincidence due to facet parameters.
In connection with the previous figure (Fig 5), additional numerical experiments show that when we increase the number of sampling points the M 3 -R2 method reaches an accuracy of 2-3%. On the other hand, the M 3 -Linear result improves at slower rate with increasing number of subdivisions. This is demonstrated, in a more general way, in Fig. 6. This is a case of two touching (shared-edge) facets that have 60 • opening angle between them. The abscissa shows the different number of subdivisions and sampling points relevant to the two M 3 methods, and the ordinate corresponds to the relative deviation in F 21 view factor from the reference value, M 4 . M 1 and M 2 methods do not depend on subdivision, therefore they keep a constant value. The first thing to note is that the M 3 -Linear method with a subdivision of zero recovers the M 1 view factor value. Then, as the number of subdivisions increases the value approaches the M 4 reference, however, this value can reach below about 5% accuracy only for 256 sub-facets. We note that the M 3 -Linear can increase subdivisions only in steps of 4 n for n = 0, 1, ..., hence the segmented curve in Fig. 6. The M 3 -R2 method, as it depends on quasi-random sequence of points, does not behave deterministically for a small number of subdivisions and does not recover the M 1 value. In this, as well as other numerical experiments, we observe that a reliable accuracy is reached when using 50-80 L. Rezac and Y. Zhao: View factor accuracy sampling points. However, as the number of sampling points is increased, it does approach asymptotically the M 4 value and may be marginally better than the M 3 -Linear approach when number of points exceeds 300-400. The decision of how many sub-facets or samplings to use is the key aspect of both M 3 methods from the point of view of accuracy, but also from the point of view of computational speed.
Finally, we apply the M 1 , M 4 , and both of the M 3 methods for view factor calculations for a section of the 67P/CG SHAP7 (Preusker et al. 2017) shape model. The extracted surface shown in Fig. 7 belongs to the northern hemisphere (includes parts of the Seth region) and is made of 1000 facets. It is chosen because it has several interesting topographical features desirable for testing the view factor algorithms. In particular, there is the deep pits/sink hole, for which accurate view factor calculations are important for a long-term evolution models. In Fig. 8 we show results of the accuracy comparisons.
The first panel in Fig. 8 shows the percent difference between M 1 and M 4 calculated as (M 1 -M 4 )/M 4 . As expected, the M 1 method in this example is applied beyond the scope of its assumptions, and the deviation inside the pit are on average around 15%, with couple of smaller facets near the floor and wall of the pit exceeding 50%. In the rugged region surrounding the pit, with many adjacent facets with large opening angles, it also performs Fig. 6. Relative deviation of view factor F 21 with respect to M 4 for two facets sharing an edge with a 60 • opening angle between them, shown as a function of subdivisions of the two triangles. Only the two M 3 methods are dependent on the number of subdivision. Both M 3 -Linear and M 3 -R2 become more accurate with larger number of subdivisions. We find that around 100 subdivisions the deviation from the reference is about 10% for M 3 -Linear and about 5% for the M 3 -R2. However, the M 3 -R2 is unpredictable for very small number of subdivisions (see text for details). Fig. 7. Small subset of SHAP7 shape model of 67P/CG featuring a large pit to be used for performance inter-comparison of the view factor calculation method. The pit can be found in the northern hemisphere of the comet (border of Hapi and Seth region). The color scale corresponds to the total F 12 view factor for each facet evaluated with the M 4 method, which is used as a reference (true value). In this calculation the effects of the surrounding shape model are not taken into account. However, this should not carry importance for facets inside the pit, which is our main focus in the following comparisons. The bright blue colored facet indicates total view factor below the minimum value of 10 −8 .
poorly. The top right and bottom left panels show the M 3 -Linear and M 3 -R2 methods with 64 sub-sampling, respectively. We observe that for facets inside the pit the accuracy at this level of subdivision is acceptable and both methods provide a similar level of accuracy within 0-3%; only few facets have larger deviations. Outside the pits, the view factors absolute values are smaller and mutual facet geometry leads to larger discrepancies of about 20%. Increasing the number of subdivisions to 256 carries a larger computational burden than the more accurate M 4 ap- Fig. 8. Inter-comparison of the various methods of estimating the view factors for a 1000 facet segment dominated by a large pit. In each of the four panels the color scale represents absolute value percent differences in total view factor of a facet with respect to the M 4 method. The limits on the color scale are 0 to 5%, and facets above this limit are shown in yellow. This includes a couple of facets at the edges of the extracted surface that do not have the visibility of any other facets. The labels in each panel indicate a particular method to which the results belong. The numbers in parentheses indicate the number of sub-facets for the M 3 methods, for example, M 3 -Lin(64) means that M 3 -Linear with each facet split into 64 sub-facets. A comparison of the "combined" method aimed to optimize speed and accuracy is presented in the bottom right panel, which is detailed in Appendix B.
proach, hence it is more suitable to employ the combined method (see appendix B). This method also offers a good trade-off between accuracy and calculation speed as summarized in Table  1.

Conclusions
The goal of this paper is to present quantitative comparisons of several numerical algorithms for calculating view factors. These are routinely required for accurate thermophysical modeling of cometary and asteroid surfaces. In particular, we want to investigate challenging cases of shared-edge facets as a function of the angle between them. Having the accuracy of these calculations under control becomes especially important when modeling energy dissipation inside cracks, pits, and/or fractures. As a benchmark of accuracy we provide the double-line integral method, M 4 , which does not suffer shortcomings of assumptions inherent in the more typical technique of a facet subdivision. The five methods investigated in this work include the "standard" ap- Table 1. Computational time of each method for the section of 67P/CG surface presented in Fig. 8. The reported value is the output of the time command of Unix. The (64) and (256)  First, we reduced the problem to a two-facet system to clearly demonstrate the performance of each approach as a function of a single parameter, N s , that combines the distance and areas of the two facets (Eq. 9). Second, we present a view factor accuracy as a function of opening angle between the triangles with a shared edge, and we also show how the performance depends on number of subdivisions and samples for the two M 3 methods.
Finally, we also applied the main algorithms to a subset of 67P/CG surface and discussed the accuracy and computational burden for each method. The main conclusions of this work can be summarized as follows: 1. The view factor calculations transformed into closed loop integral, M 4 , is the most accurate and it is especially suitable when high accuracy calculations are required for touching or very nearby facets. This is, however, very time consuming, more than two orders of magnitude compared to the most simple solution, M 1 . 2. The applicability of the simple approach, M 1 , method, depends on the N s value and to some degree on the relative orientation of the facets in question. For example, an accuracy of 1% for parallel facets seems to require N s parameter larger than about 100. However, for facets with smaller opening angles (45 • in our case) in our setup these conditions could reach values of 500-600 (see Fig. 4). Nevertheless, if an accuracy of 10-20% is acceptable the M 1 method appears to be generally proper for conditions of N s > 100, provided that facets do not share an edge. 3. We find the accuracy of M 2 approach not significantly different from M 1 except in cases of few special geometries. We do not see an advantage of using this method for general application to comet or asteroid shape models. 4. The M 3 -Linear approach, often used in applications to comets or asteroids (e.g., (Hu et al. 2017) and references therein), is a viable option for view factor calculations of nearby and/or adjacent facets. However, the number of facet subdivisions used to improve the accuracy plays a crucial role and depends on facet configurations, as shown in Fig.  5 and 6. If the number of sub-facets needs to be increased beyond 64 we found that the combined method presented in appendix B offers a better choice in terms of accuracy and computational speed. 5. The M 3 -R2 method also turns out to perform well, being able to deal with nearby or shared-edge facets. For the same number of sampling points it is very close in accuracy to the M 3 -Linear approach and even may outperform it for selected cases of geometry. This method has the advantage that the number of sampling points can be optimized easily, it is simpler to implement, and the computational burden is about 20% less than the M 3 -Linear. However, this approach relies on a quasi-random sequence of points distributed over a facet, hence, it cannot be used for number of sampling points smaller than about 60 (an empirical finding). In addition, this method does not behave predictably for low number of sampling points. 6. In larger scale applications where the view factor accuracy is a crucial aspect we suggest combining the methods M 1 , (M 3 -Linear or, M 3 -R2) with M 4 to strike a good balance between accuracy and speed. We present an algorithm in appendix B, which we used in this work to reach an average accuracy inside the pit better than 2% (Fig. 8). The algorithm is robust in that it relies on the N s metric to decide which method to apply and how many subdivisions to use. However, there are still several tuning parameters (the constants used in the algorithm) that allow a user the freedom to adjust to trade between accuracy and speed, depending on the nature of the problem at hand.
As a final remark, we acknowledge that our numerical implementation of the discussed algorithms is not optimized for calculation speed. The absolute, perhaps even the relative times, might be improved with a more rigorous approach (taking advantage of implementation in compiled programming language, hardware acceleration, and/or better parallelization). This work made use of the open source software PyVista (Sullivan & Kaszynski 2019) to produce Fig. 8 and for efficient M 3 -Linear facet subdivision implementation. (left) One hundred sampling points distributed on a triangle via the R2 method, and (right) 100 sample points from a uniform random number generator mapped onto the same triangle. The low discrepancy (the largest possible distance among the points) of R2 method is well demonstrated, producing a homogeneous coverage that allows for a quick, implicit area subdivision. dot product between the two normal vectors of facets in question to select which method is applicable. The N s also determines the number of sub-samplings (or sub-facets) when either of the M 3 methods is used. We would like to point out that for the linear method the subdivision can be scaled only as 4 n , and in this algorithm we do not allow n smaller than two, effectively yielding N d of either 8 or 64. Also, we opted out for the M 3 -Linear instead of the slightly faster M 3 -R2 method due overall smaller scatter in accuracy as shown in In this comparison we selected facets with F i j > 0.01, which includes all the facets inside the large and small pits (Fig. 7). The circle marker represents the mean percent difference with respect to the M 4 reference and error bars correspond to the 1 sigma standard deviation. Although on average the method M 1 shows acceptable accuracy (∼2.5%) the size of the error bars indicates rather large deviations among the facets; this is also seen in Fig. 8, where many facets show significantly larger deviation. Increasing subdivision number of M 3 -Linear method from 64 to 256 marginally improves the overall accuracy and reduces the deviations among individual facets. A similar conclusion is reached for the R2 method, however, it has noticeably larger deviations compared to the M 3 -Linear for the same number of subdivisions. The combined method, as given in the algorithm in appendix B, provides excellent agreement and very small deviations (which are about the size of the marker). This plots highlights the advantage of using the combined method, which balances speed and accuracy.