Propagation Attenuation Maps Based on Parabolic Equation Method

Modern wireless communication systems use various technological solutions to increase the efficiency of created radio networks. This efficiency also applies to radio resources. Currently, the utilization of a radio environment map (REM) is one of the directions allowing to improve radio resource management. The REM is increasingly used in emerging mobile ad-hoc networks (MANETs), in particular military tactical networks. In this case, the use of new technologies such as software-defined radio and network, cognitive radio, radio sensing, and building electromagnetic situational awareness made it possible to implement REM in tactical MANETs. Propagation attenuation maps (PAMs) are crucial REM elements that allow for determining the ranges of radio network nodes. In this paper, we present a novel algorithm for PAM based on a parabolic equation method (PEM). The PEM allows determining the signal attenuation along the assumed propagation direction. In this case, we consider terrain topography to obtain a more realistic analysis. Then, we average the adjacent attenuation profiles defined for the selected directions in places where attenuation has not been calculated. To this aim, linear regression is applied. Finally, we define several metrics that allow for the accuracy assessment of determining the PAM as a function of its dimensions.


Introduction
Currently, we are observing the development and implementation of the fifth-generation (5G) mobile network, i.e., 5G New Radio (NR) [1,2]. At the same time, 5G standards are still being developed for other types of communication systems than terrestrial-wireless, such as satellite [3], wired, optical, hybrid [4], etc. This direction of development is so attractive that attempts are also made to replicate it in military communication systems [5][6][7][8].
A few years ago, the first works and analyses began on the legitimacy of using 5G technologies for defense. 5G technologies are a set of radio and network technologies, i.e., in relation to the physical and higher layers of the ISO open systems interconnection (OSI) reference model, respectively. These technologies often also include technologies that were developed earlier, e.g., with the Long Term Evolution (LTE), LTE-Advance (LTE-A), or LTE-A Pro standards, and which are also implemented in the 5G NR. The most important 5G technologies include massive multiple-input-multiple-output (massive-MIMO), beamforming, millimeter waves, full-duplex communications, software-defined radio (SDR), radio resource management (RRM), dynamic spectrum access (DSA), interference mitigation, energy harvesting communications, ultra-dense network, self-organizing network (SON) [9][10][11][12][13].
• determining the attenuation for given terrain profiles on specific azimuths outgoing radially from the transmitter (TX); • using the PEM [34][35][36] for attenuation calculation along with the profiles; • using linear interpolation to determine the attenuation values between the profiles to determine the entire PAM. According to the best knowledge of the authors, the proposed approach can be considered original. The choice of PEM and linear interpolation resulted from the limited computing power of the devices on which REM will be implemented. In the general case, these methods may also be replaced by others, which provide similar abilities.
The basic PEM approach allows considering the terrain topography (with or without building heights), the electrical ground, and air-refraction parameters for determining propagation attenuation (i.e., path loss). In this paper, we present the novel PAM creation algorithm. To the best of the authors' knowledge, it is an innovative and original solution. Typically, the PEM is used to determine the attenuation along a given propagation direction related to a specific terrain profile as a function of height above ground level. We show how to use the PEM at a given height of the receiving antenna to determine the attenuation on the surface.
The remainder of the paper is organized as follows. Section 2 provides a brief overview of the REM. The idea, short description of PEM, and algorithm of the proposed PAM approach with exemplary results are shown in Section 3. In Section 4, we analyze the efficiency of creating the PAM as a function of its dimensions. Next, in Section 5, an evaluation of the PAM accuracy versus angular resolution of terrain profiles is shown. A summary of the paper is presented in Section 6.

Radio Environment Maps
The REMs are a way of representing complex and multi-domain information for cognitive radio. By sharing and disseminating information contained in the REM, it is possible to provide cognitive functionalities such as situational awareness, deducing, learning, planning, and decision-making support. Consequently, the REM can be a database distributed among secondary users, centralized, or in a hierarchical topology where the central database interacts with local databases. REM are multi-domain databases containing terrain information, sensor data, simulation results of propagation models, regulations, and policies. PAMs are used in the planning of MANETs (also as cognitive radio networks) to increase the efficiency of such a network [37]. Using the measurements anywhere in the network area is impractical. Therefore, as part of the REM, spectrum measurements from available radio network nodes are collected. Then, data fusion is performed, i.e., available measurements are converted to estimate the level of interference in other, unexplored locations. The REM contributes to the development of cognitive mechanisms and the building of long-term knowledge. The REM introduces environmental awareness that is more difficult to achieve by the individual capabilities of individual network nodes. Therefore, the REM enables network nodes to be transformed into smart ones [22,23]. REM facilitates the adaptation of the radio network to a new environment. Therefore, the REM is a promising concept for the efficient operation of a cognitive network without overloading the cognitive node and can be connected to a network-based cognitive engine. Due to the REM dissemination, an individual secondary user can get to know the radio environment much more than is the case with the limited sensing possibilities. In this case, it is possible to use cooperative techniques to interact with primary users, avoid the hidden nodes problem, and increase the entire system's efficiency. Due to its architecture, the REM can be broadly divided into global (G-REM) and local (L-REM). The REM complexity depends on the number of nodes and channels, the network coverage area, and the information granularity. REM synchronization is crucial for the accuracy and reliability of the information provided. Data must be effectively disseminated to avoid information that is obsolete or out-of-date.
The basic functionality of the REM is to build a dynamic interference map for selected frequencies in every interesting place. For this purpose, maps are used that represent the coverage of a given area with a radio signal, called radio-frequency (RF) layers of the REM, or shortly RF-REM [19]. For example, it is achieved by spectrum measurements obtained from sensors (i.e., network nodes) at specific locations. Since it is impractical and often impossible to perform measurements anywhere in the area of network operation, the available measurements are used to estimate the level of disturbances in other locations. Another approach is to use propagation damping prediction methods. It is a universal solution that does not depend on the current location of the sensor and allows you to analyze the area in which the network did not work before. Hence, the use of different propagation models for this purpose is common practice. The methods of creating the RF-REM can be divided into three main categories [19]: • direct methods based on the interpolation approach, using spatial statistics and measurements of signal levels made in specific locations, make a direct estimation of the missing data (e.g., nearest neighbor, inverse distance weighting (IDW), and kriging); • indirect methods based on the use of the TX location and its parameters as well as propagation models; • hybrid methods combining both of the above approaches.
The above groups of methods have their advantages and disadvantages and are characterized by different computational complexity. The selection of an appropriate method is also dependent on the availability of specific data, e.g., the TX location and parameters, measurement results of signal levels in selected locations, variability of parameters over time, etc. Propagation methods allow for the creation of more accurate RF-REMs [38]. The developed solution based on the PEM can be used in indirect methods of creating the RF-REM.

Propagation Attenuation Map
PAMs show the path loss values relative to the defined point, usually located in the center. This part presents the proposed algorithm for PAM creation. It includes the general idea, algorithm view, and exemplary results.

Idea of PAM Algorithm
The main idea of the proposed PAM is based on the determination of a sparse matrix for given terrain profiles on specific azimuths outgoing radially from the TX. In the next step, the full PAM matrix is determined based on the interpolation method. We suggest using PAM in the first stage and linear interpolation-in the second stage. However, it should be emphasized that both the PEM and linear interpolation can be replaced by other methods that will allow for determining the full PAM matrix. Our choice resulted from the low computational expenditure of these methods.
The main idea of the PAM algorithm is presented in Figures 1 and 2. In this case, the PAM is identified with the path loss matrix, PAM, whose center cell corresponds to the TX location. In the first stage (see Figure 1), the so-called sparse PAM matrix is created based on the PEM for the terrain profile along a specific azimuth. In the second stage (see Figure 2), the sparse PAM matrix is completed based on interpolation of the non-zero values of the matrix cells. In this way, the final PAM is obtained, represented by the so-called full PAM matrix.

Parabolic Equation Method
For modeling tropospheric ducting propagation (TDP), the mode theory, RTM, and PEM are usually used. The first two methods require greater computational effort and have significant problems with modeling radio wave propagation in the troposphere, analysis of range-dependent environments, and higher frequency ranges. Moreover, the RTM does not allow for precise estimating of the field strength [36]. Hence, the PEM is widely used as an effective tool for modeling the TDP considering the terrain topography. Therefore, we used the PEM in the proposed PAM algorithm.
The development of the PEM and its widespread use took place in the early 1980s. However, the origins of the PEM as applied to the TDP modeling date back to 1946. Then, Leontovich and Fock [39] presented an analysis of TDP based on a solution to the parabolic equation. The first solutions were based on the first-order approximation of Taylor expansion, considering only the first two terms. This resulted in a low accuracy that allowed only analyzing a narrow beam in the elevation plane. This method was called narrow-angle PEM [36].
In 1973, Hardin and Tappert [40] proposed using the split-step Fourier transform (SSFT) algorithm to sequentially solve a parabolic equation along a given radius. This approach is one of the commonly used methods alongside two other numerical methods for solving a parabolic equation, i.e., a finite-difference method (FDM) and finite-element method (FEM).
From the TDP viewpoint, the PEM allows considering the terrain irregularity, the earth conductivity, and air refractivity for determining distributions of the electric field strength or signal attenuation. The phenomenon of air refraction is considered by using different refractivity profiles, which depend on changes in air parameters (i.e., temperature, pressure, and humidity) at different altitudes. The earth's surface is also an important factor influencing the shape of these profiles, i.e., different nature of changes in the air refractivity index is above the surface of a sea, forest, mountains, desert, etc. [41]. This factor is also important in the analysis of the earth's electrical parameters. These parameters are considered in the lower impedance boundary conditions for the parabolic equation [42,43]. The lower boundary conditions are applied at the contact level of the air and soil layers. In practice, when analyzing the TDP at greater distances, the earth's curvature radius should be additionally taken into account [41].
Generally, the PEM is based on a numerical solution of the 2D parabolic equation in the Cartesian system: where E(x, z) is electric field strength for the assumed radius-terrain profile (generalized x coordinate) and for any height (z coordinate), i = √ −1, k 0 = 2π/λ is the wavenumber in the vacuum, λ is the wavelength, n = √ µ r ε r is the refractivity index of air, µ r and ε r are the relative magnetic and electric permeability of air, respectively.
Equation (1) is a parabolic approximation of the Helmholtz wave equation, whose full-wave solution is provided by the PEM [43]. This equation is solved numerically using one of the three mentioned methods (i.e., SSFT, FDM, or FEM) in a sequential manner from the TX along an analyzed terrain profile determined on a given azimuth direction of a propagation radius. Detailed descriptions of the algorithms are presented in [34,35]. A free FEM implementation in the MATLAB environment is also available [44,45].
We used the PEM implementation in MATLAB based on the SSFT approach, which considers the refractivity profile, terrain topography, impedance boundary conditions, and the earth curvature radius. In this case, we use the digital terrain elevation data (DTED) for determining the terrain profiles [46,47]. From the viewpoint of the proposed PAM solution, any PEM implementation can be used. Considering additional factors (e.g., the height of buildings, and vegetation) only influences the accuracy of the estimation of field strength or signal attenuation, but it does not matter for the implementation of our PAM algorithm.
Analyzing the literature, we see some similarities of the PAM algorithm to the 3D PEM approach presented in [48,49]. However, it should be highlighted that we do not determine the field strength distribution in the 3D space, but only 2D attenuation distribution (i.e., PAM) at the specific height of the receiving antenna. Such PAM as a heat map allows determining the radio coverage or the range of individual MANET nodes. In both cases, i.e., in the PAM and 3D PEM, radial rays from the TX (i.e., terrain profiles for which field distributions are determined) are applied. These rays are determined for specific azimuth directions with a given step.

PAM Algorithm
The proposed approach allows determining path loss values in all discrete points around a fixed TX position (marked in red in Figure 1). The indexes of individual pixels are correlated with geographic coordinates, and the distance between them depends on the available topography data, i.e., DTED resolution, ∆R. The analyzed area is contained in a square with a side length of 2R 0 , where R 0 means the radius of the analyzed terrain area for which the PAM is determined. The size of the created PAM matrix is K × K, where K = 2R 0 /∆R + 1. Calculations according to the PEM method are performed on defined directions around the TX depending on the angular resolution ∆α.
An example procedure for determining attenuation on a selected azimuth direction α n = n∆α is depicted in Figure 1 (marked in green). The first step is to compute the terrain profile along the α n direction. This profile contains the height of the terrain with ∆z resolution for discrete distance steps (∆x) along R n = R 0 /AF(α n ), where AF(α n ) is the following angular factor: AF(α n ) = sin α n for 45 PEM calculations based on the SSFT approach are performed in the next phase. Then, path loss vectors for a given height of the receiving antenna, h R , are recorded into the PAM matrix. The method of determining the final PAM is explained in detail by Algorithms 1 and 2 for the sparse and full PAM matrices, respectively. Algorithm 1 shows a determination way of the sparse PAM matrix based on PEM calculations. The procedure for determining the path loss values should be repeated for every defined azimuth and the corresponding terrain profile. Depending on the selected azimuth resolution, the number of loop executions (lines 4 to 10 in the algorithm) vary.
11. Output PAM as sparse matrix based on PEM for selected terrain profiles. The full PAM matrix is obtained by launching Algorithm 2. The proposed approach allows determining empty spaces (i.e., zero value cells) in the sparse PAM matrix (white ones in Figure 1). For this aim, we use linear regression. The interpolation method for creating the full PAM is explained in Figure 2. The input for this procedure is a partially completed path loss matrix, i.e., the sparse PAM matrix. Attenuation values in gray-marked points (see Figure 1 or Figure 2) are calculated according to Algorithm 1 described above. Zero cells (white ones in Figure 1) of the sparse PAM matrix are calculated based on linear regression and the two attenuation values located in the nearest non-zero cells. Depending on the coordinates of the considered point, i.e., indexes of matrix cell, there are two possible ways of linear interpolation: • in a column-for the azimuth directions:

Algorithm 2 (Creating full PAM matrix based on interpolation)
Require: sparse PAM matrix (PAM). Interpolation in columns: Find gth row (g ≥ k) corresponding first non-zero cell in kth column (see one of yellow cells in Figure 2).

3.
Find hth row corresponding next non-zero cell in kth column (see one of yellow cells in Figure 2).

4.
Calculate all zero cells in kth column between gth and hth rows, i.e., PAM(j, k) for g < j < h, based on Equation (2) (see blue cells in Figure 2). 5.
Set g = h.
Find gth row (g ≥ K-k + 1) corresponding first non-zero cell in kth column (see one of yellow cells in Figure 2).

8.
Find hth row corresponding next non-zero cell in kth column (see one of yellow cells in Figure 2). 9.
Calculate all zero cells in kth column between gth and hth rows, i.e., PAM(j, k) for g < j < h, based on Equation (2) (see blue cells in Figure 2). 10. Set g = h.
until h ≤ k.
repeat 13. Find pth column (p ≥ j + 1) corresponding first non-zero cell in jth row (see one of orange cells in Figure 2). repeat 14. Find gth column corresponding next non-zero cell in jth row (see one of orange cells in Figure 2). 15. Calculate all zero cells in jth row between pth and qth columns, i.e., PAM(j, k) for p < k < q, based on Equation (3) (see purple cells in Figure 2). 16. Set p = q.
18. Find pth column (p ≥ K-j + 2) corresponding first non-zero cell in jth row (see purple cells in Figure 2).
repeat 19. Find qth column corresponding next non-zero cell in jth row (see purple cells in Figure 2). 20. Calculate all zero cells in jth row between pth and qth columns, i.e., PAM(j, k) for p < k < q, based on Equation (3) (see purple cells in Figure 2). 21. Set p = q. until q ≤ j-1.

Set
23. Output PAM as full matrix (without zero cells) based on interpolation.

Exemplary Results
To check the developed algorithms, two positions of the TX have been proposed (see Figure 3). One of them is located in a lowland area (52.51 • N, 18.48 • E), the other-in a hilly terrain (51.22 • N, 15.55 • E). For each transmitting point located in the center of the map in Figure 3, terrain profiles are determined for the propagation directions (i.e., azimuths), which are shown by rays marked in red. In this case, the elevation maps obtained based on DTED2 [46,47], i.e., with an average resolution ∆R = 30 m, are illustrated in the background.  The path loss read from the determined field distributions at the receiving antenna height above the terrain profile is the basis of the created PAM matrix. Figure 5b presents exemplary path losses versus terrain profile radius (i.e., length) for the two analyzed profiles and field distributions from   The path loss read from the determined field distributions at the receiving antenna height above the terrain profile is the basis of the created PAM matrix. Figure 5b presents exemplary path losses versus terrain profile radius (i.e., length) for the two analyzed profiles and field distributions from Figure 4.    Figure 6. In this case, we show PAMs (i.e., full matrices) for two selected TX locations (see Figure 3) in lowland and hilly areas, respectively. These PAMs were obtained based on 360 terrain profiles with a step =°Δ 1 . α The PAM quality depends on the density of profiles on the basis of which it is generated. On the other hand, we may see the differentiation of the obtained PAMs depending on the terrain and TX position for which they are determined.  The path loss values versus terrain profile radius, L(D), for analyzed receiving antenna height, h R , and calculated field distributions (e.g., see Figure 5b) are written into the PAM matrix in the cells corresponding to the analyzed azimuth, α n . In this way, the sparse PAM matrix, PAM(j, k), is created. It is the basis for determining the full PAM matrix by applying the interpolation process described in Algorithm 2. The final result of the developed algorithms is illustrated in Figure 6. In this case, we show PAMs (i.e., full matrices) for two selected TX locations (see Figure 3) in lowland and hilly areas, respectively. These PAMs were obtained based on 360 terrain profiles with a step ∆α = 1 • . The PAM quality depends on the density of profiles on the basis of which it is generated. On the other hand, we may see the differentiation of the obtained PAMs depending on the terrain and TX position for which they are determined. is created. It is the basis for determining the full PAM matrix by applying the interpolation process described in Algorithm 2. The final result of the developed algorithms is illustrated in Figure 6. In this case, we show PAMs (i.e., full matrices) for two selected TX locations (see Figure 3) in lowland and hilly areas, respectively. These PAMs were obtained based on 360 terrain profiles with a step =°Δ 1 . α The PAM quality depends on the density of profiles on the basis of which it is generated. On the other hand, we may see the differentiation of the obtained PAMs depending on the terrain and TX position for which they are determined.

Efficiency of Creating PAM versus Its Dimensions
This part of the paper is devoted to studying the effectiveness of creating PAMs according to the previously presented algorithms. First, the metrics are defined, and then the research results are presented.

Metrics
As seen in the analyzed areas, there are points (i.e., sparse PAM matrix cells) for which calculations have not been made. Path loss values in these locations must be interpolated according to Algorithm 2. On the other hand, there are points (i.e., PAM cells), where calculations were performed repeatedly (more than once). This is due to the high density of the profiles near the TX. Let us introduce the redundancy parameter r as the number of path loss value calculations at a given point. As a consequence, it is possible to define the calculated points ratio CPR for the analyzed area where CP and TP are the number of calculated points (with r ≥ 1) and total number of points at the analyzed area, respectively.
To determine the number of all calculations performed, the redundancy values for all points in the analyzed area should be summed up where s is the point index and r s is the redundancy value for a specified point s.
The CP measure introduced earlier considers the number of points, where the calculations were made at least once (i.e., r ≥ 1). By correlating this parameter with the redundancy r, it is possible to define a new CP r measure meaning the number of points where the calculations were made exactly r times. Therefore, the previously introduced CN metric can also be expressed as where r max is the maximum value of r for the analyzed area. The normalized number of points where the calculations were made exactly r times can be defined as follows: NCP r (%) = CP r CP ·100.
Additionally, for evaluating the PAM creation efficiency, we use a computation time CT expressed in seconds. The measure considers the execution of all algorithm steps presented in Section 4.2, excluding the calculation of terrain profiles. The absolute values of this measure depend on the processing power of the used computer. Therefore, this parameter is used to compare different variants, e.g., different sizes of the analyzed area.

Result Analysis
The results presented in this part of the paper relate to studies carried out for various sizes of the analyzed area (square with a side length of 2R 0 ). The angular resolution ∆α for all tests are equal to 1 • . Figure 7 shows the effect of increasing the radius R 0 on the calculated points ratio (i.e., CPR) values (blue color) and the calculation time (i.e., CT) needed to obtain the results (orange color).
As mentioned in the previous section, the CT metric depends on the processing power and is used for comparison purposes. The greater the value of the radius R 0 , the longer the time needed to perform the calculations. For example, in the case of a tenfold increase in R 0 (from 2 to 20 km), the CT increases about six times. In the case of analyzing small areas, the vast majority of points on the PAM are calculated based on PEM (some points even multiple times, as shown in the following figures-see Figure 8). The larger the area of the analysis, the CPR value significantly decreases, and more points must be determined using the interpolation method presented in Figure 2. In the case of the radius R 0 equal to 1 km,  As mentioned in the previous section, the CT metric depends on the processing power and is used for comparison purposes. The greater the value of the radius 0 , R the longer the time needed to perform the calculations. For example, in the case of a tenfold increase in 0 R (from 2 to 20 km), the CT increases about six times. In the case of analyzing small areas, the vast majority of points on the PAM are calculated based on PEM (some points even multiple times, as shown in the following figures-see Figure 8). The larger the area of the analysis, the CPR value significantly decreases, and more points must be determined using the interpolation method presented in Figure 2. In the case of the radius 0 R equal to 1 km, about 90% of the points on the PAM are determined based on calculations. Increasing 0 R to 10 km results in a decrease in CPR to the level of only about 21%. In this case, 79% of points must be interpolated. Figure 8 depicts the effect of duplicate calculation of points in generating the final PAM. In particular, the greatest redundancy of calculation is for the area next to the TX. Therefore, the increase in the number of calculations (i.e., CN ) does not translate proportionally to the increase in .  As mentioned in the previous section, the CT metric depends on the processing power and is used for comparison purposes. The greater the value of the radius 0 , R the longer the time needed to perform the calculations. For example, in the case of a tenfold increase in 0 R (from 2 to 20 km), the CT increases about six times. In the case of analyzing small areas, the vast majority of points on the PAM are calculated based on PEM (some points even multiple times, as shown in the following figures-see Figure 8). The larger the area of the analysis, the CPR value significantly decreases, and more points must be determined using the interpolation method presented in Figure 2. In the case of the radius 0 R equal to 1 km, about 90% of the points on the PAM are determined based on calculations. Increasing 0 R to 10 km results in a decrease in CPR to the level of only about 21%. In this case, 79% of points must be interpolated. Figure 8 depicts the effect of duplicate calculation of points in generating the final PAM. In particular, the greatest redundancy of calculation is for the area next to the TX. Therefore, the increase in the number of calculations (i.e., CN ) does not translate proportionally to the increase in . CPR  Another metric, independent of the processing power of the used computer, is the total number of calculations (i.e., CN). The absolute value of this measure for three distances R 0 is presented in Figure 9a. The normalized values of the calculations for these three cases are shown in Figure 9b.
Another metric, independent of the processing power of the used computer, is the total number of calculations (i.e., CN ). The absolute value of this measure for three distances 0 R is presented in Figure 9a. The normalized values of the calculations for these three cases are shown in Figure 9b.

PAM Accuracy Evaluation versus Angular Resolution of Terrain Profiles
In this section, we present a preliminary assessment of the accuracy of the proposed algorithm versus the angular resolution of the terrain profiles. This assessment was made for the selected TX location with the coordinates ( )°°4 9.28 N, 19.96 E and = 0 1 km. R  Figure 10 shows the normalized number of points on the map for which the attenuation calculations are repetitive. In this case, the parameters NCP 2 and NCP 3 mean that 37% and 14% of points are calculated two or three times, respectively. On the other hand, NCP >5 means the sum of points on the map is calculated 6 or more times. For PAM with a larger radius, i.e., R 0 = 5 km, the number of points without redundancy (NCP 1 ) is 92% and and the number of points without redundancy (NCP 1 ) is 97%.
Another metric, independent of the processing power of the used computer, is the total number of calculations (i.e., CN ). The absolute value of this measure for three distances 0 R is presented in Figure 9a. The normalized values of the calculations for these three cases are shown in Figure 9b.

PAM Accuracy Evaluation versus Angular Resolution of Terrain Profiles
In this section, we present a preliminary assessment of the accuracy of the proposed algorithm versus the angular resolution of the terrain profiles. This assessment was made for the selected TX location with the coordinates ( )°°4 9.28 N, 19.96 E and = 0 1 km. R Figure 10. Normalized number of points where calculations were made exactly r times for TX-RX distance equals to 1 km.

PAM Accuracy Evaluation versus Angular Resolution of Terrain Profiles
In this section, we present a preliminary assessment of the accuracy of the proposed algorithm versus the angular resolution of the terrain profiles. This assessment was made for the selected TX location with the coordinates (49.28 • N, 19.96 • E) and R 0 = 1 km. The terrain around TX is shown in Figure 11a, while Figure 11b illustrates the PAM obtained for the angular resolution of the terrain profiles equal to ∆α = 1 • . In our analysis, this PAM is used as a reference to determine the algorithm errors at another angular resolution of the terrain profiles. Figure 12 shows exemplary PAMs for the analyzed area and two different angular resolutions equal to 2 • and 10 • . The terrain around TX is shown in Figure 11a, while Figure 11b illustrates the PAM obtained for the angular resolution of the terrain profiles equal to =°Δ 1 . α In our analysis, this PAM is used as a reference to determine the algorithm errors at another angular resolution of the terrain profiles. Figure 12 shows exemplary PAMs for the analyzed area and two different angular resolutions equal to 2° and 10°. In the following, for every map point, we calculate the attenuation error are the PAM values (i.e., attenuations) obtained for the analyzed Δα and reference =°Δ 1 α angular resolutions of terrain profiles. In this way, the error maps, i.e., attenuation error matrices, are determined. Figure  13 shows sample error maps that were obtained for the two PAMs illustrated in Figure  12. The terrain around TX is shown in Figure 11a, while Figure 11b illustrates the PAM obtained for the angular resolution of the terrain profiles equal to =°Δ 1 . α In our analysis, this PAM is used as a reference to determine the algorithm errors at another angular resolution of the terrain profiles. Figure 12 shows exemplary PAMs for the analyzed area and two different angular resolutions equal to 2° and 10°. In the following, for every map point, we calculate the attenuation error are the PAM values (i.e., attenuations) obtained for the analyzed Δα and reference =°Δ 1 α angular resolutions of terrain profiles. In this way, the error maps, i.e., attenuation error matrices, are determined. Figure  13 shows sample error maps that were obtained for the two PAMs illustrated in Figure  12. In the following, for every map point, PAM(j, k), we calculate the attenuation error where PAM ∆α (j, k) and PAM ∆α=1 • (j, k) are the PAM values (i.e., attenuations) obtained for the analyzed ∆α and reference ∆α = 1 • angular resolutions of terrain profiles. In this way, the error maps, i.e., attenuation error matrices, are determined. Figure 13 shows sample error maps that were obtained for the two PAMs illustrated in Figure 12. These error maps are the basis for determining the empirical cumulative distribution function (CDF), CDF(∆L), and root-mean-square error (RMSE), RMSE(∆L), of the attenuation error. Exemplary CDFs for the two analyzed error maps (see Figure 13) are presented in Figure 14a. Generally, we can judge that the obtained functions are symmetrical about the sign of ∆L. This means that the attenuation errors in the PAM with ∆α oscillate around the attenuations for the reference PAM. Hence, we also computed the CDFs for the module of the attenuation error, CDF(|∆L|), which are depicted in Figure 14b. The CDFs analysis shows that increasing Δα causes a significant increase in the attenuation errors in determined PAMs. This is due to the smaller number of terrain profiles included in the calculations. For The CDFs analysis shows that increasing Δα causes a significant increase in the attenuation errors in determined PAMs. This is due to the smaller number of terrain profiles included in the calculations. For The CDFs analysis shows that increasing ∆α causes a significant increase in the attenuation errors in determined PAMs. This is due to the smaller number of terrain profiles included in the calculations. For ∆α = 2 • , 80% of the absolute attenuation errors are less than 3 dB, while for ∆α = 10 • , it is 7 dB. For CDF(|∆L|) = 0.9, |∆L| = 5 dB and |∆L| = 10 dB for the angular resolutions equal to 2 • and 10 • , respectively. Therefore, we may assume that for small Da, the number and time of calculations can be reduced at the expense of the accuracy of the PAMs obtained. This approach can be used in the preliminary stage of PAM determination, which requires initial visualization of the results. In the next step, with the platform computing resources available, the exact PAMs can be determined for the lower angular resolution of terrain profiles.
Similar conclusions might be drawn by analyzing the RMSE scalar measure. Figure 15 shows the RMSE of attenuation error for selected angular resolution. be reduced at the expense of the accuracy of the PAMs obtained. This approach can be used in the preliminary stage of PAM determination, which requires initial visualization of the results. In the next step, with the platform computing resources available, the exact PAMs can be determined for the lower angular resolution of terrain profiles.
Similar conclusions might be drawn by analyzing the RMSE scalar measure. Figure  15 shows the RMSE of attenuation error for selected angular resolution. respectively. Moreover, the graph shows that the error difference for is small. In this case, we can obtain almost unchanged PAM accuracy with a two-fold reduction in the number of terrain profiles and a two-fold reduction in computation time.
The accuracy of PAM generation depends not only on the angular resolution of the profiles but also on their length (i.e., 0 R ), i.e., the PAM matrix dimension. A more detailed analysis of the PAM accuracy will be presented in our next work.

Conclusions
This paper focuses on a novel method of creating PAM as a crucial REM element, which allows for determining the ranges of radio network nodes. Our solution is based on the PAM determination in two stages. In the first stage, we determined the so-called sparse matrix for profiles outcoming radially from the TX. In the second stage, we used the interpolation method to calculate the full PAM matrix. The proposed approach is based on the DTED-based terrain profiles, PEM, and linear interpolation algorithm. Thereby, the developed solution considers the influence of topography on path loss calculation. The proposed algorithm was implemented in the REM for determining the radio ranges of tactical MANET nodes in the emerging military communication system. On the other hand, the proposed two-step approach to the PAM can be considered universal. This means that DTED, PEM, and linear interpolation can be replaced by others, which provide calculating the full matrix.
In this paper, we have presented a detailed description of the PAM generation algorithm based on the PEM for a given terrain topography. Metrics for evaluating the algorithm effectiveness depending on the PAM size have been defined. They were the basis for the assessment of the proposed solution for various PAM. Analyzing the results presented in the paper, it can be concluded that in the 0 R range up to about 3 km, over 50% For the angular resolutions equal to 2 • and 10 • , we obtained RMSE(∆L) = 4 dB and RMSE(∆L) = 7 dB, respectively. Moreover, the graph shows that the error difference for ∆α = 5 • and ∆α = 10 • is small. In this case, we can obtain almost unchanged PAM accuracy with a two-fold reduction in the number of terrain profiles and a two-fold reduction in computation time.
The accuracy of PAM generation depends not only on the angular resolution of the profiles but also on their length (i.e., R 0 ), i.e., the PAM matrix dimension. A more detailed analysis of the PAM accuracy will be presented in our next work.

Conclusions
This paper focuses on a novel method of creating PAM as a crucial REM element, which allows for determining the ranges of radio network nodes. Our solution is based on the PAM determination in two stages. In the first stage, we determined the so-called sparse matrix for profiles outcoming radially from the TX. In the second stage, we used the interpolation method to calculate the full PAM matrix. The proposed approach is based on the DTED-based terrain profiles, PEM, and linear interpolation algorithm. Thereby, the developed solution considers the influence of topography on path loss calculation. The proposed algorithm was implemented in the REM for determining the radio ranges of tactical MANET nodes in the emerging military communication system. On the other hand, the proposed two-step approach to the PAM can be considered universal. This means that DTED, PEM, and linear interpolation can be replaced by others, which provide calculating the full matrix.
In this paper, we have presented a detailed description of the PAM generation algorithm based on the PEM for a given terrain topography. Metrics for evaluating the algorithm effectiveness depending on the PAM size have been defined. They were the basis for the assessment of the proposed solution for various PAM. Analyzing the results presented in the paper, it can be concluded that in the R 0 range up to about 3 km, over 50% of points in the PAM are calculated according to PEM. In the case of larger areas, it should be remembered that most of the points are interpolated, which affects the accuracy of determining the path loss values. These errors will get bigger with increasing distance from the TX. On the other hand, we presented a preliminary analysis of the PAM accuracy versus the angular resolution of the terrain profiles. The obtained results show that increasing ∆α reduces the accuracy of determining the attenuation. However, with a two-fold increase in ∆α, the computation time decreases twice, and the mean estimation error is equal to 4 dB. The obtained results also showed that increasing ∆α from 5 • to 10 • only slightly increases the deviation of ∆L. The utilization of larger ∆α can be used in the initial stage of PAM determination for the preliminary visualization of the results. In the next step, the exact PAMs may be determined for the lower angular resolution of terrain profiles.
Future work on PAM development will focus on assessing the impact of profile density on the path loss estimation accuracy, PAM utilization to evaluate the radio range of nodes (i.e., network coverage), comparison of the linear interpolation method with others (e.g., IDW or kriging), and PAM algorithm modification for the TX directional antennas.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to project restrictions.