Inconsistent Rays in Propagation Prediction by Ray Launching in Rectangular Tunnels

Ray tracing is replacing the less accurate statistical and empirical approaches of radio channel modeling. Although being high-frequency asymptotic method, ray tracing has been shown to be mathematically equivalent to purely analytical modal methods in rectangular tunnels and waveguides. The equivalence applies to modeling reflections using image theory, while other variants of the ray tracing algorithm are still subject to approximation errors that can systematically add up. Failure to recognize this leads to indiscriminate use of the less appropriate algorithmic variants. In particular, the use of ray tracing by ray launching can be questionable in tunnel environments, at least when characteristic sequences are used to avoid double counting errors. In addition to the known path inaccuracies, we first identify previously untreated inconsistent rays as the most problematic, leading to a significantly overestimated signal at distances greater than 100m. We show that the problem is not equivalent to double counting of rays since the inconsistent rays represent valid wavefronts for the points in space at which they are detected. The discrepancy arises from the use of reception spheres, which allow some spatial displacement of ray paths. We quantify the errors in terms of estimated power, channel impulse response, and delay spread in a rectangular tunnel. We propose an improvement in the double-count filters to detect inconsistent rays. However, evaluation of signals deeper in the tunnel at frequencies below 1GHz must still be avoided by ray launching due to the remaining path inaccuracies.


I. INTRODUCTION
Radio frequency propagation modelling relies on increasingly sophisticated techniques. Stochastic modelling has been largely replaced by discrete channel models, with ray tracing variants leading the way. On the other hand, numerical solutions based on Maxwell's equations remain computationally infeasible and have many weaknesses that still need to be addressed to compete with ray methods for electrically large problems.
Since the initial proposals by Ikegami et al [1], numerous ray tracing algorithms have been proposed. Most of them can be divided into three computationally distinct groups. The first group effectively traces a large number of The associate editor coordinating the review of this manuscript and approving it for publication was Gerardo Di Martino . rays from the transmitting source in all directions into the scene. Algorithms in this group refer to the principle as raylaunching [2], pincushion method [3], or shooting and bouncing ray (SBR) [4], which are terms used here interchangeably. The second approach combines traced rays as ray tubes [5] or beams [6] to reduce the computational cost. Finally, the image theory [7] goes the furthest by using entire scene surfaces as ray aggregation units. Although the above methods are all based on the principles of geometrical optics, they cannot be applied indiscriminately to propagation problems. Each method has its own advantages and disadvantages, and the correct method should be chosen for a particular problem.
Ray tracing by ray launching can be used to simulate various propagation phenomena such as reflection, refraction, diffraction, and scattering. The magnitude and phase of the transmitted electric field are updated as the ray travels through the scene, taking into account attenuation due to the energy spreading in space, antenna gains and polarizations of the transmitters and receivers. Rays are launched without requiring any prior knowledge of the positions of the receivers. It is common practice to distribute rays as evenly as possible in all directions to scan the unknown environment without favoring parts of the modeled geometry.
Reception spheres are used to detect rays passing in close proximity to points of interest [8], [9]. Unfortunately, the use of spheres leads to various errors that have serious consequences for the accuracy of the simulation. Rays belonging to different wavefronts should be identified, failing to do so leads to so called double counting error [10]. Moreover, path inaccuracies are generally assumed to be negligible and are not treated, which is not the best approach for high accuracy simulations. Furthermore, as shown here, inconsistent rays occur, which can cause unacceptably high errors in regular geometries such as rectangular tunnels or waveguides.
Image theory is better suited for modelling pure reflections in scenarios with large planar surfaces [11], [12], [13]. Systematic enumeration of ray trajectories between known transmitters and receivers based on the position mirroring technique avoids all the errors mentioned above. In general, intersection tests between ray trajectories and mirror surfaces still need to be performed. Although the solution is quite elegant, it can only be applied to reflection phenomena and does not scale well with massively parallel computer architectures, on which forwarding of rays through the scene can be very efficient [14]. The latter advantage and the wide availability of SBR tools are most likely the reasons that many analyses of tunnel environments are performed using less accurate approaches, even for problems for which the image method is much better fit. Often, the brute-force method is further optimized, e.g., by shooting and bouncing polygons instead of rays [15] or by launching only a subset of rays in precomputed directions [16], [17]. It is not uncommon to introduce further simplifications just to reduce the high computational cost, resulting in even more inaccurate simulations [18]. Many publications rely on variants of the SBR tool in tunnels without providing information on how the errors attributed to the reception spheres are handled [19], [20], [21], [22], [23]. The fact that some widely used commercial software offers the possibility to choose the SBR method does not mean that it is the best choice for the considered scenario. However, as shown here, when comparing the results with field measurements, where some error is naturally to be expected, good agreement is nevertheless likely for shorter axial distances and higher simulated frequencies.
This paper analyzes several types of errors specific to the use of reception spheres in ray tracing by ray launching, which can lead to substantial deviations from the theoretically correct electric fields and possibly to incorrect conclusions of the underlying study. The analysis is limited to rectangular tunnels and reflection phenomena because the errors are significantly amplified in an environment where the superposition of ray paths is governed by specific geometric patterns.
Moreover, the rectangular dielectric tunnel is an environment for which the equivalence of image theory and modal method has been demonstrated, providing an excellent opportunity to accurately quantify errors independent of other error sources that are present in more complex environments. The case is of practical importance since many studies of rectangular tunnels choose ray launching simulations as part of their evaluation. However, the limited scope of the following analysis does not preclude the types of errors identified here from occurring in other environments.
In particular, three types of errors are examined: path inaccuracy, double counting, and inconsistent rays. While double counting has been extensively studied in the literature and efficient techniques exist to avoid it, the other two types of errors have not received much attention. Path inaccuracies have been recognized before, but rarely corrected. As far as we know, the inconsistent ray error has neither been recognized nor treated. We propose to modify the characteristic sequences of the double counting avoidance to prevent the addition of inconsistent rays.
In summary, the contributions of this work can be described as follows: • A previously unknown type of error, called inconsistent rays, is identified that occurs in ray tracing by ray launching associated with the use of reception spheres, which cannot be treated as a double counting error.
• The errors are analyzed systematically and in isolation in the environment of a dielectric rectangular tunnel or equivalent waveguide with a known solution.
• Double counting avoidance algorithm is adapted to include inconsistent rays.
• Tunnel simulation scenarios that should be avoided if the SBR is the tool of choice are identified.
The remainder of the paper is organized as follows. Section II describes the tunnel channel model based on image theory, which is later used as a reference for error evaluation. The electric field along the tunnel is given assuming a sufficiently large axial distance. After discussing the necessary requirements for the radius of the reception sphere, the sphere-induced errors are presented in Sect. III, including path inaccuracies III-A, double counting of rays III-B and newly identified inconsistent rays III-C. The SBR tool used in the simulations is described in Sect. IV, including ray launching template and the architecture of ray tracing pipeline. The error analysis presented in Sect. V addresses path inaccuracies and inconsistent rays, while skipping the already solved problem of double-counted rays. It also describes the possible avoidance of inconsistent rays and quantifies the remaining errors by the SBR tool. Section VI concludes the paper.

II. TUNNEL CHANNEL MODEL
We adopt the analytical approach from [24], where similar modelling of dielectric tunnel was used to show the mathematical equivalence of image-based ray tracing and modal methods. Although the modelling applies to VOLUME 10, 2022 three-dimensional space, the enumeration of all ray paths is here conveniently achievable in a two-dimensional plane shown in Fig. 1. The modelling plane is perpendicular to the tunnel at the point of transmission S. The origin of the plane coordinate system is set at the center of the tunnel rectangle. The right-side vertical wall intersects x-axis at (a, 0) whereas the ceiling intersects y-axis at (0, b), thus having a tunnel 2a wide and 2b high. Z -axis is directed along the tunnel length.
For given regular geometry one can analytically enumerate all source images as coordinate pairs (x m , y n ) where m and n are the number of reflections from the vertical and horizontal walls, respectively, and (x 0 , y 0 ) is the location of the source S. An image is a virtual location of a source as seen by the receiver after a signal undergoes a specific number of reflections from horizontal and vertical walls. The order of reflections is important. Negative m indicates left wall as the last vertical reflection. On the other hand, n takes negative value for the last bottom wall reflection with respect to the sequence of horizontal reflections. According to the image theory, straight lines between virtual source images and a reception point represent all signal propagation paths that need to be considered when calculating the electric field. In other words, there are neither duplicate nor missing paths for a theoretically correct signal evaluation. For the sake of completeness and to ensure the reproducibility of the results, we next review the evaluation of a signal power distribution in a rectangular tunnel with smooth surfaces. The notation used in the following is also summarized in Table 1.
Assuming a sufficiently large axial distance z, the direction of the electric field along the tunnel can be neglected, since grazing incidences do not change it to the extent that is important for the following analysis. A similar assumption is made in the proof of equivalence of ray and modal methods in [24]. The electric field at the observed point is then a scalar superposition of the electric fields of all propagation paths where the contribution of a single path is Here, E t is the magnitude of the transmitted electric field, k is the wavenumber, r is the distance between the receiver and the source image, ρ ⊥ is the perpendicular reflection coefficient, and ρ || is parallel reflection coefficient. In practice, only a limited number of paths are summed, with the shorter paths being preferred over the longer ones.
For the source images at z = 0, the distance simplifies to The reflection coefficients can be expressed as where the angles of incidence relative to a reflecting surface normal follow from the geometry as and θ || = arccos |y n − y| r m,n .
⊥ and || capture the electrical properties of the reflecting surface, i.e., relative permittivity ε and conductivity σ , as complex relative permittivitȳ and are calculated as and

III. SPHERE-INDUCED ERRORS
Reception spheres are used to detect rays passing in close vicinity of the receiver in ray-launching approach. The slight deviations of discovered propagation paths lead to some serious consequences elaborated in the following. It should be noted that determining the exact propagation path after multiple reflections and possibly other modelled interactions with matter would be too complex for a ray launching that aims for independent and scalable handling of rays on massively parallel computer architectures. The size of a reception sphere is usually hidden from the end user because it is dictated by the wavefront ray density at a given propagation distance. A concept of the wavefront and its ray density is needed to properly accumulate the contribution of the rays to the observed electric field. The initial wavefront is represented by all emitted rays. The interaction of rays with a flat surface triggers a new wavefront of rays. In general, different interaction phenomena produce different wavefronts, e.g., reflected wavefront, transmitted wavefront, edge diffracted wavefront. The size of a reception sphere should guarantee at least the interception of one ray per wavefront. Since the spatial separation of the rays increases with the propagation distance, the spheres should also increase in size, which means that spheres of different sizes must be considered at a given point in space.
Assuming an ideal distribution of rays with uniform spacing αd, where α denotes a constant angular distance or FIGURE 2. Field strength differences between two alternative ray paths in brute-force ray tracing due to path inaccuracy are attributed to differences in path lengths and reflection angles. separation between rays in radians and d denotes the propagation distance, the adjacent rays form the vertices of an equilateral triangle on the wavefront surface. To enclose at least one ray, a randomly placed reception sphere should have a minimum radius of the triangle-circumscribed circle Since uniform distribution of rays in three dimensions is not possible for more than 12 rays, (12) must account for the variability in spatial separation of rays, e.g., by using the maximum angular distance. Unfortunately, the use of reception spheres with any nonzero radius is a source of errors, which can be classified into three types: path inaccuracy, double counting and inconsistent rays.

A. PATH INACCURACY
Path inaccuracy causes errors in the electric field contribution to the total field due to the spatially displaced ray path. We illustrate the problem in Fig. 2 using our tunnel model. Two rays are shown, both of which experience a single reflection from the left tunnel wall. The first path of length d a leads directly through the reception point. It corresponds to the path found by the image theory. The second path of length d b is found by a ray launching and intersects the reception sphere at some distance from the receiver. In general, the length d a deviates slightly from the length d b . However, the main cause of inaccuracies is different angles of incidence at the interaction points, which lead to slightly wrong Fresnel coefficients and generally affect both signal attenuation and polarization. In Fig. 2, only different angles in the xy-plane are shown. The angular difference due to path inaccuracy is also present in the z-axis.
Path inaccuracies are difficult to eliminate. We could strictly maintain small radii of reception spheres by using as many rays as possible, which keeps the error insignificant VOLUME 10, 2022 for practical purposes, as will be shown in Section V. Note, however, that the efficient double-counting filters presented next allow the use of spheres with radii larger than the minimum. In this case, the error in the evaluated electric field may become unacceptable even for a large number of emitted rays.

B. DOUBLE COUNTING
The size of a reception sphere should guarantee at least one hit per wavefront in order not to underestimate the electric field. On the other hand, a radius of the reception sphere given by (12) does not exclude the possibility of several hits per wavefront. This leads to a doubling of the field strength, i.e., to a power error of +6 dB per erroneously counted ray. In the ideal case of uniformly distributed rays and planar wavefront, the probability that a randomly positioned sphere spans more than one ray is 2/ , or about 20.9% [8], which is often reported as the double-count probability. We have shown in [10] that the above value is only a lower bound for several reasons. Apart from the fact that a uniform distribution of rays is not possible and a larger spatial distance must be used in (12), modelling refracted and diffracted phenomena in general scenarios change spacing between rays in non-trivial ways, thus making the theoretical double-count probability more difficult or even impossible to achieve. In [10], we have given analytical bounds on the maximum and minimum for rays distributed in space using icosahedral grid lattice.
The double counting problem has been studied in depth and efficient solutions exist. The Zero Counting (ZC) algorithm [25] shrinks the sphere to avoid double hits, but omits a number of valid wavefronts. Some of these wavefronts are caught by introducing auxiliary rays with smaller reception spheres. The probability of missing a wavefront is estimated to be 4.97% for uniformly spaced rays. However, this completely ignores the smaller spacing between rays due to the rigorous interaction modelling. The SBR, which allows for zero counting, could still experience double hits, while the total number of emitted rays must be tripled.
In the Distributed Wavefronts (DW) algorithm, all rays incident on the sphere contribute a certain fraction of the signal power. The contribution is weighted so that the signal levels are approximately invariant with respect to the local distribution of rays and the best average fit is obtained [8].
Exact detection of double counting is based on so-called characteristic sequences (CS) [26]. A characteristic sequence is a sequence of scene interactions that rays experience as they traverse the scene. Interactions with the same planar surface or with the same diffraction edge are considered identical. The sequence is constructed in real time while the ray is moved through the scene. It consists of unique surface identifiers. In case of simulating diffraction phenomena, unique identifiers are also assigned to the edges and included in the sequence. The characteristic sequence can be viewed as a wavefront signature that is matched by the receiver against a table of previously seen signatures. Unique rays are allowed to contribute to the signal, while rays that match the table entry are simply ignored.
The high resource requirements of handling signature tables and processing characteristic sequences has led to replacing a sequence with a pair consisting of the ray incidence angle and the number of interactions, which is justified in [27]. The proposed algorithm declares two rays with the same number of interactions and similar sphere incidences as belonging to the same wavefront. This is clearly an approximate solution. Similar space and time considerations apply to the variant of the algorithm in which the number of interactions is replaced by the propagation distance [2]. The fact that even rays of the same wavefront have slightly different angles of incidence, and hence path lengths, makes this wavefront differentiation even less accurate.
In [10], we proposed Bloom filters configured with a low false positive rate to replace exact wavefront differentiation and solve the problem more efficiently. Instead of long characteristic sequences represented as a concatenation of unique interaction identifiers, a small number of independent hash values of finite length must be computed as the ray moves through the scene, typically on the order of 10. The table of characteristic sequences is replaced by a lossy hash table of a filter that requires less memory. Bloom filtering provides near-optimal avoidance of double counting, i.e., it is as efficient at detecting double hits as CS, but requires significantly less memory and processing time.

C. INCONSISTENT RAYS
The above solutions to double counting do not address the third type of error that arises from the use of reception spheres. We first became aware of this type of error when modelling propagation in a tunnel environment using the ray launching approach, because of its severe impact. However, the error also occurs when modelling general environments, but with less obvious consequences.
The problem of inconsistent rays can be illustrated as follows. Two rays shown in Fig. 2 intersect the dashed vertical line only once. The crossed line corresponds to the left wall reflection. In general, the intersections of the vertical grid lines indicate the order in which the left and right wall reflections alternate. Similarly, the intersections of the horizontal grid lines indicate the order of the lower and upper wall reflections. Thus, based on the intersection ordering, the exact sequence of the modelled reflections can be determined. This sequence corresponds to the characteristic sequence used to avoid double counting. Let us now consider the ray paths shown in Fig. 3. The inconsistent ray error occurs when the ray paths from a given image intersect grid lines in different orders. Note that the double-counting filters only exclude ray paths with the same order of interactions, which is perfectly valid for the exact paths, but may fail when the paths end in close proximity to a receiver.
In the case of source image I −2,1 in Fig. 3, for example, two rays with different characteristic sequences can be seen. The upper ray is first reflected from the right wall, followed by the upper and right wall reflections. Finally, the ray crosses the reception sphere. The lower ray, on the other hand, is reflected from the tunnel walls in a different order, starting with the upper wall reflection and followed by the right and left wall reflections. Obviously, the double-counting avoidance allows both rays, resulting in incorrect field strength. As shown in the next section, this error occurs frequently in tunnel geometry, with the calculated field deviating significantly from the correct value.
The problem of inconsistent rays cannot be solved with the existing double counting prevention, since the latter relies on characteristic sequences. Moreover, the inconsistent rays represent valid wavefronts for the points in space where they are actually detected. The question arises whether using the image location with some margin-instead of characteristic sequences-resolves the error at least for the reflected rays. In the case of rectangular tunnels, the source images are well separated for such an approach and provide nearly correct field strengths, as shown below. Note that the path inaccuracies described earlier are still present.

IV. SBR TOOL
The ray tracing tool used in the evaluation is a variant of SBR with initial rays emitted in an icosahedral grid pattern [28]. Icosahedral grids are commonly used to distribute the rays as evenly as possible around the transmission point. The choice of a particular ray launching grid is important because it determines ray spacing and thus the sizes of reception spheres. In general, there are two variants of icosahedral grid construction algorithms, recursive and non-recursive [29], with slightly different distribution of rays in space. In Fig. 4, both approaches are illustrated. We use the first variant, where each recursive step bisects the sides of the icosahedron and projects newly introduced points onto the enclosing sphere, as shown in subfigures (a) through (c). Subfigure (d) shows an alternative non-recursive construction. The number of recursive steps n is called the refinement level, which establishes the final 10 × 2 2n + 2 grid points that serve as initial ray directions. For the recursive icosahedral grids, the exact maximum angular separation α n at the refinement level n is given by where the initial angular separation of icosahedral grid points α 0 , measured from the center of the icosahedron, is The SBR tool is based on a highly optimized GPU core using the NVIDIA OptiX ray tracing engine [14] adapted for radio frequency simulations. Fig. 5 shows a call diagram of the NVIDIA ray tracing pipeline including the adaptations for radio channel modeling. The blue colored boxes are OptiX internal algorithms common to all ray tracing applications, from computer graphics to scientific modeling and visualization. The yellow boxes are user-defined functions, in our case radio frequency propagation prediction. Clearly separated modules that integrate seamlessly into the ray tracing pipeline generate the initial rays using an icosahedral grid template, perform basic intersection tests with the triangles describing the environment and with the reception spheres, as well as the channel modeling tasks in the two separate closest hit procedures. The closest hit calls are invoked by the engine itself. The triangle closest hit module computes the reflection, refraction, or diffraction of a single ray at the points of electromagnetic interaction with objects in the scene. The sphere closest hit module, on the other hand, accumulates the power delay profile at the receiving point and includes highly configurable Bloom filtering to avoid double counting of rays, as described in [10]. Both procedures can generate new rays by recursively calling rtTrace, either to simulate a new wavefront or simply to force the continuation of a ray in case a reception sphere is hit. The objects in the scene VOLUME 10, 2022 FIGURE 5. Custom NVIDIA OptiX ray tracing engine for radio channel modeling is at the core of our SBR tool. The user-defined yellow boxes in the control flow shown perform specific tasks for modeling radio propagation, while the blue boxes are internal to OptiX. Some functions of the general ray tracing pipeline are not needed (empty boxes).
are held entirely on a GPU in a bounding volume hierarchy, with rays generated in parallel threads and traced through the scene. The bounding volume hierarchy allows ray tracing to be accelerated by restricting intersection tests to only plausible parts of the geometry. For more details on ray types, their payloads, acceleration structures and scene representation, see [14] and NVIDIA's OpiX Programming Guide [30].
The assumption of a large axial distance that we made in Section II does not apply, which means that the electric fields are treated as vectors in three-dimensional space. For transmission and reception, the orientation and gain of the antennas are taken into account. However, in order to compare the simulation with the modal analysis from [24], an isotropic antenna is used in the following analysis. The number of receivers observed in a single simulation run is limited only by the available memory. With many receivers scattered in space, the brute force approach has an advantage over the tools that launch only a subset of rays in precomputed directions. Such tools should either be run multiple times or use a subset of rays that covers almost all directions.

V. ERROR ANALYSIS
We use the same tunnel geometry as in [24] where the signal along the tunnel was calculated using image theory and the modal method to prove the equivalence of both approaches at larger axial distances. The authors of [24] also validated the results with field measurements [31], making the presented scenario an excellent reference. The tunnel is 1.83 m wide and Vertically polarized omnidirectional antennas are assumed. Fig. 6 shows the scenario considered with two selected ray paths experiencing several reflections before hitting a reception sphere. Note that the size of a reception sphere is greatly exaggerated, and although a single sphere is drawn at the reception point, each ray path has its own sphere size. The simulation parameters are summarized in Table 2.

A. SIGNAL POWER
In Fig. 7, the signal power along the tunnel is plotted for the selected frequencies at a transmit power of 0 dBm. The plots are the result of a Matlab program that implements image theory using the modelling approach for a rectangular tunnel described in Section II. The order of images or modes in the modal theory used to calculate Fig. 7 was 100, which conceptually corresponds to the maximum number of reflections configured in the SBR tool.
Icosahedral grid refinement levels from 10 to 14 were used for error assessment, with the upper limit chosen by the maximum run time of several hours. The number of rays, maximum angular separation and radii of the reception spheres   at propagation distances of 1 and 600 m, determined by (12), are listed in Table 3. Note that the propagation distances after multiple reflections can be larger than 600 m, which means that the reception spheres can be even larger.
To evaluate path inaccuracy error independently of other types of errors, the Matlab program implementing the image theory and used to plot Fig. 7 was extended as follows. The straight lines between the virtual source images and the receiver points were randomly shifted at the receiver point in a range given by the size of a reception sphere. The point closest to the center of the sphere is chosen as the end point of the ray path. The randomly generated points in a sphere follow a uniform distribution. In Fig. 8, the signal including the path inaccuracy error is plotted for the four frequencies and reception sphere sizes at refinement level 13. Visible deviations from the plots in Fig. 7 are observed at the two lower frequencies. The received power of the 455 MHz signal degrades significantly around the 200 m distance, with the calculated power remaining roughly constant deeper in the tunnel, while the signal at 915 MHz begins to fluctuate at the end of the simulated length.
Quantitative assessment of path inaccuracy error is given in Table 4 for five refinement levels, including standard deviation, maximum, and average of errors. Note that the refinement level is used here only to define the sizes of the reception spheres. As expected, using smaller spheres reduces the path inaccuracy. The two higher frequencies, 2.45 and 5.8 GHz, show minimal errors in the range of a tenth of a dB for the observed tunnel length. Caution should be exercised when interpreting the results in context of the SBR tool. The grid pattern of emitted rays gives a systematic superposition of the rays at the receivers, which may not be as random as assumed in our experiment. This may lead to an increase in path inaccuracy error, as shown at the end of this section.
The double counting errors have been treated extensively in the literature. We refer the interested reader to [10] for a comprehensive survey and evaluation. The use of bruteforce ray tracing in tunnel scenarios should not be affected by double-counting errors if the SBR tool implements the proposed exact algorithms.
The evaluation of the third type of errors caused by inconsistent rays was performed using the SBR tool. In this way, the spatial dependencies of the rays are taken into account. This type of error cannot be studied completely in isolation because in the SBR tool path inaccuracy cannot be easily  eliminated. However, since we know the magnitude of the error from our first experiment, we can safely draw conclusions about the magnitude of errors caused by inconsistent rays.
The received power shown in Fig. 9 corresponds to the same scenarios as in Fig. 8. Here, the discrepancies between the ideal and calculated values are even more pronounced, with the power for the two lowest frequencies being significantly overestimated at distances greater than 100 m. As can be seen from Table 5, the error remains largely independent of the refinement level. With respect to frequency, it increases by up to two orders of magnitude for the two highest frequencies, resulting in a standard deviation of several dB, while the calculated power for the two lowest frequencies beyond 100 m has no scientific or practical value. One can improve the brute-force treatment of reflections to remove inconsistent rays while maintaining resistance to double-counting errors by redefining the characteristic sequences as follows. Instead of creating a sequence of reflection surface identifiers [26] or updating the hash of the sequence [10], the source image is associated with each ray and updated for each reflection event. At the receiver, the coordinates of the source image should be rounded to a certain precision. This step is necessary to guarantee the exact same location even if the order of reflections changes and the limited floating point accuracy affects the calculations. Of course, the simulated scenario should not have source images closer than the selected accuracy. The coordinates of such a source image then serve as a key for searching and updating the database or, in the case of avoiding double counting by Bloom filters, as a basis for computing multiple independent hash values. We implemented the latter solution in our SBR tool and rounded the coordinates to the nearest centimeter. The improved ray tracing significantly reduced the remaining error, as shown in Fig. 10 and Table 6.
The similarities to Fig. 8, where the path inaccuracy error was emulated by a Matlab code, are high with the SBR results showing greater variation, especially at the lower frequency of 455 MHz. In addition to the path inaccuracies due to regular grid patterns instead of random offsets, other implementation-dependent sources of error also contribute to the differences between Fig. 8 and Fig. 10, such as the use of single floating point precision and the accuracy of the actual routines for ray-triangle intersections. However, the standard deviations and averages of the remaining error for the two highest frequencies are less than 1 dB, suggesting that ray tracing by ray launching with efficient treatment of inconsistent and duplicate rays may be suitable for some analyses in rectangular tunnels at higher frequencies.

B. CHANNEL CHARACTERIZATION
In addition to the presented signal power calculation, rigorous error correction is important for advanced channel characterization. Therefore, in the following we analyze the effect of inconsistent rays on the estimated channel impulse response, delay spread, and angle of arrival. We have chosen a receiver located 300 m from the transmitter, since at this distance the ray launching already shows significant deviations of the signal power from the theoretical values.
Since the channel impulse response (CIR) fully describes a linear, time-invariant communication channel at all frequencies, and CIR measurements have become an essential part of pulse-based physical layers, there is a need for accurate CIR simulations. The impulse response of the channel is usually modeled as the sum of a time-varying number of multipath components in a tap-delayed linear filter. A time-invariant model is used here because the time-varying aspects of the channel modeling are not relevant to a static simulation scenario. The time-invariant CIR can be formulated as where each of L taps represents a multipath component with sign-extended (±) real amplitude A multiplied by a timedelayed Dirac delta function. Knowing (15), one can compute a received signal as a convolution of CIR with the transmitted signal, i.e., the received signal is a sum of attenuated, delayed and possibly overlapping versions of the transmitted signal, plus white Gaussian noise. The computation of h(τ ) by ray tracing is not as straightforward as one would expect, since ray paths, amplitudes, and propagation delays are functions of frequency. In ray tracing, the signal is evaluated at a single frequency, so only a narrowband CIR can be computed that applies to bandwidth-limited transmission. In [32] and [33], subband-split ray tracing was proposed to fully evaluate (15). The method involves multiple simulation runs for a set of center frequencies of complementary subbands. These narrowband CIRs are then combined in frequency space. An inverse Fourier transform yields the final CIR over a wider bandwidth. In the following, only simulations of narrowband CIRs are presented and analyzed.
Field measurements do not produce individual CIR taps with arbitrary temporal precision, as simulations do. The actual resolution is the inverse of the sampling pulse bandwidth [34]. Instead of individual taps, we obtain bin amplitudes that are functionally dependent on these taps. To analyze the inconsistent ray error in this context, power delay profiles with 5 ns resolution (200 MHz bandwidth-limited pulse) were calculated and compared. Fig. 11 compares two narrowband impulse responses at the lower frequency of 455 MHz for a selected receive point; one is from image theory, the other was computed with the SBR tool. The latter contains inconsistent ray error in addition to negligible path inaccuracies. Shown are the first bins covering a power range of 150 dB. The average error of 9.0 dB affects most of the bins. The amplitudes of the bins decrease with delay, as expected, because the latter bins correspond to rays that experience a larger number of reflections and thus longer propagation paths and higher attenuation. The first bin, which contains the line-of-sight component, is less affected. It is the sum of the shortest paths with a small number of reflections, resulting in a less likely change in the reflection order. A similar comparison for the highest frequency of 5.8 GHz is shown in Fig. 12. The average bin error is 1.1 dB, which is hardly noticeable at the scale shown.
The root mean square (RMS) delay spread used in wireless channel characterization can be calculated from the impulse response [35]. In the presented case, the delay spread of 3.1 ns  is incorrectly extended to 3.13 ns (1% increase) regardless of the simulated frequency when inconsistent rays are present.
The angle of arrival as a measure of propagation direction is even less affected. Although some reflected paths erroneously receive a power increase because of inconsistent rays, this occurs in a rectangular tunnel from symmetric directions, at least near the centerline, which cancels the angular error.
The improved SBR that removes inconsistent rays as proposed in the previous subsection, does not deviate significantly from image theory with respect to the above metrics. The average bin error decreases to 0.04 dB at 455 MHz and to 0.01 dB at 5.8 GHz. The RMS delay spread increases by 0.03% in both examples.

VI. CONCLUSION
This work investigated ray tracing errors in the brute-force method based on ray launching that are caused by the use of reception spheres. The popularity of the general-purpose technique, which covers most propagation phenomena and can be efficiently executed on parallel computer architectures, has led to its widespread use. The method is also applied to simpler scenarios, such as reflection modelling in rectangular tunnels, where a more accurate image theory is a better choice. Failure to address the errors discussed here can lead to erroneous conclusions.
In addition to the presented signal power calculation, the effect of inconsistent rays on the channel impulse response, delay spread, and angle of arrival for a receiver deeper in the rectangular tunnel was investigated. We have shown that path inaccuracy can be kept sufficiently low by densely represented wavefronts in the GHz range. Further, double-counting errors can be efficiently removed by known solutions. On the other hand, inconsistent ray errors cannot be treated by the characteristic sequences of the double counting avoidance. This type of errors leads to unusable results at tunnel axial distances above 100 m and below GHz frequencies, while several dB errors occur above the GHz threshold. The proposed treatment of the inconsistent rays by integrating the source images into the double-count filtering eliminates most of the errors, except for the tolerable path inaccuracies, which must be taken into account when interpreting the results.
The considered errors also occur in simulations of other propagation phenomena in irregular environments. Since the irregularity prevents systematic amplification of the less common inconsistent ray errors, they are less likely to have significant or even observable effects. Nonetheless, no attempt was made here to address refraction, diffraction, or scattering. Evaluation of these phenomena, as well as ways to resolve inconsistent rays in these cases, remains a future work. In summary, cautious use of ray tracing by ray launching is advisable for regular geometries.