Consider a point source $𝐬$. The measurement plane, denoted as $ℳ$, is partitioned into $N_M$ cells, each with an equal area $\left|C\right|$. The radio map value for the $i$-th cell, where $i\in \{1,\mathrm{\dots },N_M\}$, is defined as:

Here, $\overline{e}\left(𝐭\right)$ denotes the channel gain at the position $𝐭$ on the measurement plane, considering all paths originating from the source $𝐬$. Assuming a discrete set of paths impacts the measurement plane at $𝐭$, $\overline{e}\left(𝐭\right)$ can be expressed as

where $𝒫(𝐬,\mathrm{\ell },𝐭)$ represents the set of all paths with depth $\mathrm{\ell }$ that originate from $𝐬$ and intersect the measurement plane at $𝐭$. The term $e(p)$ is defined as:

In this context, $\lambda$ is the wavelength, $\mathrm{\ell }$ is the number of scatterers along the path $p$, ${𝐂}_{\text{T}}$ is a complex-valued vector of dimension 2 representing the transmit antenna pattern evaluated at the angles of departure of $p$, $({\theta }_{\text{T}},{\phi }_{\text{T}})$, and ${𝐓}^{({\mathrm{\ell }}^{\prime })}$ is a $2\times 2$ complex-valued matrix modeling the interaction with the ${\mathrm{\ell }}^{\prime }$-th scatterer. Unlike in (4), the receive antenna pattern is not applied here. Instead, the squared norm of the electric field is used. This is equivalent to assuming that a receiver positioned on the measurement plane uses a dual-polarized isotropic antenna, and that both components are combined non-coherently.

While the definition (33) is comprehensive, it is not directly practical for computation. A straightforward but inefficient method to compute a radio map is to calculate the paths for each measurement cell using the path solver (Section 3) by placing one or more targets within each cell. However, this approach becomes computationally prohibitive as the number of cells increases. Therefore, we reformulate the radio map definition (33) to enable efficient computation via Monte Carlo integration.

To build intuition, let’s first consider the scenario where all paths from the source $𝐬$ to the measurement cells are specular chains, meaning they contain no diffuse reflections. In this case, for a given scene and source $𝐬$, a path of depth $\mathrm{\ell }$ is fully characterized by its sequence of interaction types ${𝝌}^{(\mathrm{\ell })}=({\chi }^{(1)},\mathrm{\cdots },{\chi }^{(\mathrm{\ell })})$ and its departure direction ${\widehat{𝐤}}^{(0)}$. This is because, in specular chains, the direction of the scattered wave at each interaction point is uniquely determined by the interaction type and the direction of incidence, as illustrated in Figure 23. In this case, the definition (33) can be reformulated over the unit sphere $\mathrm{\Omega }$ centered on $𝐬$, which represents all possible ray departure directions. To achieve this, we define $p(𝐬,{\widehat{𝐤}}^{(0)},{𝝌}^{(\mathrm{\ell })},𝐭)$ as the path of depth $\mathrm{\ell }$ from $𝐬$ to $𝐭\in ℳ$, characterized by the departure direction ${\widehat{𝐤}}^{(0)}$ and the sequence of interaction types ${𝝌}^{(\mathrm{\ell })}$. The length of the corresponding ray tube is denoted by $r$, and the angle between the plane normal and the path’s incident direction is denoted by $\theta$. By applying the change of variables $d𝐭=\frac{r^2}{\cos \theta }d{\widehat{𝐤}}^{(0)}$, we can reformulate the definition (33) as:

In this context, $𝐭$ denotes the position on the measurement plane where the path intersects. The double summation is over all path depths $\mathrm{\ell }=0,1,\mathrm{\dots }$ (with $\mathrm{\ell }=0$ corresponding to LoS) and all possible sequences of interaction types ${𝝌}^{(\mathrm{\ell })}$, excluding diffuse reflections. The indicator function $𝕀\left(𝐭\in C_i\right)$ equals $1$ if $𝐭$ is within $C_i$ and $0$ otherwise. If the path does not intersect the measurement plane, then $𝕀\left(𝐭\in C_i\right)=0$. The integration is carried out over the unit sphere, where $d{\widehat{𝐤}}^{(0)}$ represents the infinitesimal solid angle element corresponding to the ray tube in the direction ${\widehat{𝐤}}^{(0)}$.

To extend the definition to include diffuse reflections, we first modify the definition (36) to explicitly include the sequence of propagation directions of the scattered rays at the interaction points, denoted as ${\overline{𝐤}}^{(\mathrm{\ell })}=({\widehat{𝐤}}^{(1)},\mathrm{\cdots },{\widehat{𝐤}}^{(\mathrm{\ell })})$. For this purpose, we introduce ${\overline{𝐧}}^{(\mathrm{\ell })}=({\widehat{𝐧}}^{(1)},\mathrm{\cdots },{\widehat{𝐧}}^{(\mathrm{\ell })})$, representing the sequence of normals at these interaction points. We then define the function $\overline{f}(\cdot )$, which ensures that the direction of the scattered rays aligns with the interaction type, as follows:

where

and, $\delta (\cdot )$ is the Dirac delta distribution. The definition in (36) can be then be reformulated as

Compared to (36), the definition in (39) integrates over all propagation directions of the scattered rays at the interaction points. The integral over ${\mathrm{\Omega }}^{\mathrm{\ell }}$ combined with the function $\overline{f}(\cdot )$ merely filters out invalid directions and, strictly speaking, does not perform any other function. However, while the definition in (39) is equivalent to (36), it facilitates the extension to paths with diffuse reflections:

where the inner sum covers all possible sequences of interaction types of depth $\mathrm{\ell }$, including diffuse reflections. The path $p(𝐬,{\widehat{𝐤}}^{(0)},{\overline{𝐤}}^{(\mathrm{\ell })},{𝝌}^{(\mathrm{\ell })},𝐭)$ represents a path of depth $\mathrm{\ell }$ from the source $𝐬$ to the target $𝐭\in ℳ$. This path is defined by the initial propagation direction ${\widehat{𝐤}}^{(0)}$, the sequence of scattered ray directions ${\overline{𝐤}}^{(\mathrm{\ell })}$, and the sequence of interaction types ${𝝌}^{(\mathrm{\ell })}$. Compared to (39), the paths depend on the directions of propagation of the scattered waves, which are not solely determined by the initial propagation directions and interaction types when diffuse reflection occurs. In (40), $f(\cdot )$ in (37) is extended to include diffuse reflections as follows:

where $𝕀(\cdot )$ is the indicator function, which equals $1$ if the condition is true and $0$ otherwise. Note that in (40), $r$ is the length of the ray tube impinging on the measurement plane, which, in the case of paths including diffuse reflections, corresponds to the distance from the intersection with the measurement plane to the last diffuse reflection point.

In practice, we estimate $C_i$ using Monte Carlo integration. Similar to the path solver, we employ random sampling of interaction types (see Section 3.1.2) and limit the integration to paths with a maximum depth of :

Assuming $N_S$ samples for Monte Carlo integration, we sample the unit sphere around the source $𝐬$ using $N_S$ samples, resulting in initial ray tubes with a solid angle of approximately $4\pi /N_S$. For diffuse reflection points, we emit a single ray on the hemisphere, leading to ray tubes with a solid angle of $2\pi$. This results in the following estimator

Here, the subscript $n$ denotes the $n$-th path, and ${\omega }_n^{({\mathrm{\ell }}^{\prime })}$ represents the solid angles of the ray tubes corresponding to sampling the ${\mathrm{\ell }}^{\prime }$-th integral:

In the case of diffuse reflections, the multiplication by the ray tube solid angle is already applied through the computation of the scattered electric field (123). It has been included here for the sake of clarity. As we will see in the next section, this estimator can be efficiently computed using a single SBR loop.

The Sionna RT built-in radio map solver calculates the estimator (43) using a single SBR loop. Unlike the path solver described in Section 3, it does not require an image method-based algorithm, as it only needs to test intersections with the measurement plane rather than connecting paths to precisely located target points. Additionally, the electric field for each ray tube is computed dynamically during the SBR loop as the paths are traced.

Algorithm 3 outlines the SBR loop utilized by the Sionna RT radio map solver. While the algorithm is presented for a single path $n$, the radio map solver processes all paths in parallel. The process begins by sampling the departure direction for each ray using a Fibonacci lattice, as detailed in Section 3.1.1. The electric field is initialized using the transmit pattern evaluated at the path’s departure direction ${\widehat{𝐤}}_n^{(0)}$. The ray tube length is initially set to zero, and its solid angle is set to $4\pi /N_S$, since the $N_S$ initial rays have departure directions on the unit sphere. The radio map is initialized to zero, and the path is marked as active. The SBR loop continues to trace the paths through the scene until a user-defined maximum depth $L\in \mathbb{N}$ is reached or the ray is deactivated. The loop’s first step involves checking for intersections with the measurement plane. If an intersection is detected, the radio map is updated using the function Add_To_Radio_Map. The ray-plane intersection-checking function is not detailed in this document but is assumed to return the position of the intersection point ${𝐯}_{\text{mp}}$. The function Add_To_Radio_Map accumulates contributions to the radio map in the vector $𝐂=[C_1,\mathrm{\cdots },C_{N_M}]$ at the intersected cell entry, as described in (43). Note that by iterating over path depth $\mathrm{\ell }=0,\mathrm{\dots },L$ and adding contributions to the radio map at every iteration, the radio map solver accounts for paths of all depths, from LoS to the maximum depth. The ray tube length is updated with the distance between the last intersection point and the measurement plane intersection point. Additionally, ${\widehat{𝐧}}_{\text{mp}}$ represents the normal vector to the measurement plane, and $\left|{\widehat{𝐤}}^{𝖳}{\widehat{𝐧}}_{\text{mp}}\right|$ is the cosine of the angle between the ray direction and the measurement plane normal. The weight ${\alpha }_n$ accounts for the transmitter array, as explained in Section 4.1.2. Next, the algorithm checks for intersections with the scene. If no intersection is found, the ray is deactivated as it has exited the scene. Otherwise, an interaction type ${\chi }_n^{(\mathrm{\ell })}$ is sampled using the probability distribution from Section 3.1.2, and the radio material is evaluated for the sampled interaction type. This results in the transfer matrix ${𝐓}_n^{(\mathrm{\ell })}$ and the scattering direction ${\widehat{𝐤}}_n^{(\mathrm{\ell })}$. Unlike Algorithm 2, the function Sample_Evaluate_Interaction samples the interaction type, scattered direction, and evaluates the radio material. The electric field is updated by applying the transfer matrix ${𝐓}_n^{(\mathrm{\ell })}$ and dividing by the probability of the sampled interaction type. If the interaction type is diffuse, the ray tube length is reset to zero, and the solid angle weight is set to $2\pi$, as diffuse reflection points act as new point sources, and as a single ray is emitted in the hemisphere containing the incident direction. Note that in the case of diffuse reflection, the matrix ${𝐓}_n^{(\mathrm{\ell })}$ applies the scaling by ${\omega }_n^{(\mathrm{\ell }-1)}$. The SBR loop concludes by updating the ray tube length with the distance between the last and current intersection points. Finally, the radio map computation is completed by scaling by ${\lambda }^2/{(4\pi )}^2$ according to (35), and by the surface area of measurement cells, following the radio map definition in (33).

Figure 25(b) illustrates the computation time needed to generate a radio map for a single source as the number of samples $N_S$ and the number of transmit antennas $N_A^{\text{tx}}$ are varied. The setup, depicted in Figure 25(a), employs the “simple street canyon” scene built into Sionna RT. The red ball indicates the source, positioned 20 meters above the ground. The radio map in Figure 25(a) is calculated using ${10}^8$ samples and a linear array of 512 transmit antennas. Notably, the radio map solver can compute radio maps with ${10}^9$ samples on a single NVIDIA RTX 4090 GPU. The computation time increases with both the number of samples and antennas, exceeding one second only for the most demanding computations. It is important to note that while the computation time increases with the number of antennas due to the precomputation of the array weighting factor $\alpha$, it does not increase the memory footprint of the radio map solver (see Section 4.1.2).