A Modified Radiopropagation Multipath Model for Constant Refractivity Gradient Profiles

Abstract This paper presents a radiopropagation algorithm based on a Ray Tracing (RT) technique that combines a modified multipath model for constant refractivity gradient profiles and the Uniform Theory of Diffraction (UTD). A novel formulation is proposed by the authors for calculation and ground-reflection analysis of ray paths depending on atmospheric refractivity. The algorithm introduced herein was evaluated in a mixed scenario and in two more realistic case studies, under conditions of constant refractivity gradient and lossy terrain profiles. Pathloss results are obtained and compared with Parabolic Equation (PE) numerical solution results at 2.0 GHz, 3.5 GHz and 5.4 GHz. In such conditions, the modified radiopropagation multipath algorithm with atmospheric refractivity introduced herein showed satisfactory results.


I. INTRODUCTION
Latest technologies of Fifth-Generation (5G) wireless networks are aimed to the meeting of the objectives of high-speed transmission and large increase in channel capacity [1], [2].It is awaited that 5G systems configure heterogeneous networks of overlapping cell, also composed of macrocells for wide coverage.Several sub-6 GHz frequencies are being tested for pioneering 5G services, specifically, 750 MHz, 2.5 GHz, 3.5 GHz and 5.4 GHz bands are proposed in Latin American countries to be used in 5G applications [3].Within this framework, one of the main challenges of these technologies is to guarantee connectivity and to provide broadband services in rural and sub-urban areas.Thus, improved, quick and refined propagation models must be proposed to predict coverage and characterize the radio channel in this kind of scenarios, and at the mentioned frequency bands.This paper is organized as follows: Section II presents the formulation about ray paths considering atmospheric refractivity, Section III describes a proposed algorithm for a Modified RT technique, the case studies and results are shown in Section IV, and the conclusions are presented in Section V.

II. RAY PATHS UNDER CONSTANT REFRACTIVITY GRADIENT
Using GO assumptions, the following equation describes the ray paths depending on a refractive index n [16], where r represents the position vector and l is the coordinate associated with the path.We assume a 2-D propagation in Cartesian coordinates system (x is distance and z is height), a predominantly horizontal propagation (dl ≈ dx) and also that n(r) is a function of height, n(r) = n(z).Thus, Eq. ( 1) is rewritten as: where n(z) = n 0 + δz.In this expression, n 0 and δ are the surface value and the gradient of the refractive index of the medium, respectively [10].When n(z) is very close to 1 for all heights, it is possible to represent Eq. ( 2) as: The solution to Eq. ( 3) is given by [15]: Eq. ( 4) denotes the dependency of the ray path on the x coordinate, on the transmitter height Z T and on the launch angle of the ray α, as shown in Fig. 1.

A. Direct ray analysis
Knowing the location of a receiver point, the launch angle of the direct ray can be obtained from: where ∆Z R = δR 2 /2, Z R is the receiver height, and R is the distance from the transmitter to the receiver point (see Fig. 1).

B. Incident and reflected rays on flat terrain
The deduced formulation for ground-reflection analysis has been partially presented by the authors of this work in international conference papers.This approach is proposed when the location of a receiver point is known on flat terrain.Therefore, we deduce a formulation to calculate the location of a reflection point.Fig. 2 shows the incident and reflected ray paths and it can also be seen that incident ray and reflected ray at x = X make an angle α 1 and α 2 with the parallel to the earth plane, respectively.From Snell's law  and using the reciprocity theorem (in order to analyze a given antenna in transmitting or receiving mode), the variables exhibited in Fig. 2 and Eq. ( 5), we can write α 1 and α 2 as follows: Substituting the expressions ( 7) and (8) into Eq.( 6), we obtain a cubic equation [15]: The unknown X is the distance from the transmitter to the reflection point, and its solution is given by a real root X = x 1 , minimum value and x 1 > 0. So the reflection point is at the coordinate (X, 0), as can be seen in Fig. 2.

C. Analysis of incident and reflected rays over sloping terrain
In the same way as in flat terrain case, our purpose is to offer a formulation that calculates a reflection point over sloping terrain.The terrain is modeled by z = mx, where m is the slope of the terrain.According to Fig. 3, the following mathematical relationships are deduced: where, as depicted in Fig. 3, ∆z R is the height of a receiver point over terrain, ϕ corresponds to the angle associated with the slope of the terrain, H T refers to the transmitter height for this situation, H m is the maximum height of the wedge, X 1 represents the distance from the transmitter to the initial location of the sloping terrain, X C is the distance from the initial location of sloping terrain to the reflection point and X 2 indicates the distance from the transmitter to the final location of the sloping terrain.
Based on Fig. 3 and Snell's law, it is known that α 3 = α 4 , which can be rewritten as: where α 1 and α 2 correspond to the incident and reflected angles, respectively, with respect to a horizontal reference line (Fig. 3).Similarly, α 3 and α 4 refer to the incident angle and reflected angle with respect to sloping terrain.Manipulating mathematically, we have: Therefore, we use some trigonometric properties to reach at the following expression: Now, substituting the expressions ( 7), ( 8), ( 10), ( 11), ( 12) and ( 13) into Eq.( 16), it is possible to obtain a fourth degree equation: where the coefficients are defined in terms of δ and the variables exhibited in Fig. 3, as shown from expression ( 18) to ( 22) [17]: In Eq. ( 17) the solution is given by a real root X = x 1 , maximum value and x 1 > 0. Therefore, the reflection point is at the coordinate ( X, m(X − X 1 ) ).
The UTD approach used in this multipath algorithm proposal is based on the Lubbers coefficients.More details about the used UTD formulation are presented in Section III-C.

III. PROPOSED ALGORITHM FOR A MODIFIED RADIOPROPAGATION MULTIPATH MODEL
The algorithm is shown in Fig. 4 as a flowchart.The direct ray, reflection and diffraction analysis, are implemented by three blocks (see Fig. 4) which are executed in parallel.To run those three blocks it is necessary to know the following input parameters: electromagnetic properties, refractivity gradient of the medium, and locations of the transmitter and receiver.

Locate receiver
Obtain the direct ray path-Eq (4)

Calculate launch angle of the direct ray -Eq (5)
Direct ray hit receiver?

A. Direct ray analysis
Knowing the input parameters, the launch angle of the direct ray is calculated according to Eq. ( 5), then the direct ray path is obtained from Eq. ( 4).In the next step, it is investigated if a direct ray hits the receiver.If the algorithm identifies that the receiver is reached, it proceeds to save the path parameters, and then, continues to post-processing.Otherwise, the algorithm ignores the direct ray path and continues to post-processing.
Two situations can occur so that a ray path does not hit on a receiver point.The first situation is that the algorithm identifies that there is an obstacle that interrupts the ray path.The second one is that due to the curvature of the ray path, it cannot reach the receiver.These above criteria are applied to direct, incident and reflected ray paths.

B. Reflection analysis
The location of the reflection point is calculated by using the formulation proposed herein, Eq. ( 9) or Eq. ( 17) is used depending on the situation, either flat or sloping terrain.Once known the location where the ground-reflection occurs, it is possible to calculate the launch angle of the incident ray from Eq. ( 7), then the incident ray path is obtained.In the next step, it is investigated if the obtained incident ray path hits the location of the point reflection.If the algorithm identifies that the reflection point is reached, it proceeds to save the path parameters.Otherwise, the algorithm ignores the obtained incident ray path and continues to post-processing.
If the incident ray path hits the reflection point, it is possible to calculate the launch angle of the reflected ray from Eq. ( 8), then, the reflected ray path is obtained as indicated by Fig. 4. In the next step, it is learned if the reflected ray path hits the location of the receiver point.If the algorithm identifies that the receiver point is reached, it proceeds to save the path parameters, and then, continues to post processing.Otherwise, the algorithm ignores the obtained reflected ray path and the path parameters for incident ray, consequently, continues to post-processing.

C. Diffraction analysis
The algorithm has the ability to identify edges that generate diffraction.An analysis of direct, incident and reflected rays associated with the diffraction phenomenon is carried out.
The UTD formulation considers the 2-D problem of diffraction by a semi-infinite wedge with straight edges.The geometry and wedge diffraction variables involved in the formulation can be seen in [18].The distances between the source and the wedge, and between the wedge and the observation point, are calculated with the optical length equation [11].The general UTD solution for the electric field at the observation point is: where, E i (W ) is the incident electric field at the edge of the wedge, A(s d ) is the amplitude factor, s d is the distance between wedge and observer, and D is the dyadic diffraction coefficient.Assuming the classical notation of [19], the dyadic soft and hard is given as follows: where D i , for i = 1, . . ., 4, are the UTD diffraction coefficients, G 0 and G n , are grazing incidence factor, R 0 and R n are Fresnel reflection coefficients, for the analyzed faces of the wedge.The UTD formulation is based on the Fresnel reflection coefficients, specifically, the Lubbers Coefficients [20] defining the incidence and reflection angles of the incident and diffracted rays.

D. Post-processing
Everytime that a ray hits a receiver point, the following path parameters are computed and saved by our algorithm: τ i (t): time delay of arrival (TDA) of path, T i (t): polarimetric transmission matrix of path, Ω T,i (t): direction of departure (DoD) of path and Ω R,i (t): direction of arrival (DoA) of path.In order to calculate the total EM field at each receiver point, a coherent sum is computed, taking into account the individual contribution of each multipath component.
In the case studies carried out in Section IV, the following propagation mechanisms are processed by the algorithm: a direct ray path, an incident-reflected ray path, diffraction ray paths and diffractionreflection combinations arriving at the receiver point, accumulating in maximum two consecutive propagation mechanisms (see Fig. 8, in Section IV-B).

IV. CASE STUDIES AND RESULTS
The first case study consists of a mixed scenario to validate our proposal on a test that combines several canonical elements.Then, two realistic cases studies are carried out to verify the applicability of our radiopropagation modeling in environments that combine lossy terrain profiles and refractivity variation conditions.The frequencies used in the simulations correspond to the bands projected to offer 5G services.In all case studies, comparisons are made between a modified multipath algorithm and the DMFT-SSPE numerical solution for pathloss results.In addition, to demonstrate the ability to characterize the radio channel in realistic cases, Power Delay Profiles (PDP) are obtained through the proposed multipath algorithm.

A. Mixed scenario
The environment is a mixed scenario, with a total extension of 40 km.The terrain profile is a lossy wedge of 80 m height and centered at 20 km, and whose initial and final positions are located at 12 km and 28 km, respectively.In this case study, the ground region is characterized as a standard ground surface (ϵ r =15 and σ = 0.012 S/m) and a saltwater region (ϵ r =81 and σ = 2 S/m) is placed behind the wedge, from 28 km to 32 km.The refractivity variation is modeled as: N (z) = 304 − 100z, a linearly decreasing refractivity altitude profile, which corresponds to a surface duct [21].In this work, we use a gradient of modified refractive index dm/dz = dn/dz + 1.57 × 10 −7 [12].The transmitter is located at 100 m height, the simulation uses a frequency of 5.4 GHz and vertical polarization.
The source is modeled as a Gaussian antenna beam pattern.In most 2-D long-range radiopropagation models, this type of source modeling is chosen because it can adjust beamwidth and beam tilt and provide an approximate representation for paraboloid dish antennas [22].
In order to calculate the pathloss horizontal profile using the Modified RT algorithm, we define the horizontal step size (∆x = 10 m) and receiver points located at 10 m height along the terrain.Likewise, the vertical step size (∆z = 0.16 m) is used to calculate the pathloss vertical profile at x = 32 km.Fig. 5 shows the scenario scheme and depicts the conformation of ray paths to calculate the pathloss horizontal using the Modified RT algorithm.In Fig. 6 the 2-D received power distribution achieved with DMF-SSPE is shown.In this graph, the very low levels of received power that are perceived at the  closest distances to the transmitter are not expected, since it does not correspond to a natural behavior of the phenomenon.However, this response is due to the shape of the directional Gaussian antenna pattern used in the simulations.To estimate a pathloss horizontal profile, a correct range of receiver points is chosen, from 1.5 km to 40 km.In Fig. 7 the comparison between the proposed Modified RT and DMFT-SSPE numerical algorithm is presented.The horizontal and vertical profiles for pathloss results are illustrated, respectively, in Figs.7(a The statistics are summarized in Table I.It is possible to observe that a good agreement is obtained in statistic comparison between our proposal of a Modified RT and DMFT-SSPE results.

B. Cúcuta city case
This first realistic scenario is located in Cúcuta city (Colombia).It consists of a terrain profile with an approximate extension of 15 km, the surface is characterized as a standard ground and the refractivity variation is modeled as: N (z) = 305.66− 60z, according to Recommendation ITU-R P.453-14 for the environmental conditions of the mentioned city.The terrain profile information is collected using the NASA Shuttle Radar Topography Mission (SRTM) Version 3.0 dataset.
The transmitter is located at 100 m height and the simulation uses a frequency of 2.0 GHz and vertical polarization.To calculate the pathloss horizontal profile, both techniques use a horizontal step size (∆x = 50 m) and the receiver points are located at 10 m height along the terrain.
Fig. 8 shows the conformation of ray paths at x = 1000 m and 10 m height to calculate Pathloss and the PDP at this receiver point through the Modified RT algorithm.In Fig. 9 the PDP is depicted, in this way, the ability of the proposed algorithm to characterize the radio channel is demonstrated.It is possible to obtain a 2-D received power distribution with DMFT-SSPE, since this numerical approach allows to discretize and calculate the EM field solution of the entire analyzed scenario.Hence, we generate Fig. 10 using a DMFT-SSPE algorithm, which is implemented computationally by the authors of this work.In Fig. 10, it can be seen that at shorter distances there are Non-Line-of-Sight (NLOS) regions, due to the particularities of the terrain profile of the analyzed scenario, and taking into account that the height of the receiver points is 10 m.Furthermore, it is also possible to note that at longer distances there are line-of-sight (LOS) regions, and consequently, this propagation mechanism implies a greater contribution to the final value of received signal intensity at those receiver points.The previously described can be observed in the results of the pathloss curves obtained along the terrain profile (see Fig. 11).
In this realistic case, we also obtain results from a RT technique that does not consider atmospheric effects, denominated Non-Modified RT.The purpose is to observe if the inclusion of atmospheric refractivity in the RT technique leads to promising results, when comparing the results obtained with respect to reliable and consolidated techniques, such as the DMFT-SSPE numerical approach.Fig. 11(a In Fig. 11 it can be observed that in some positions no results are obtained, this happens because our RT algorithms use basic combinations of propagation mechanisms and therefore the ray paths do not hit on certain observation points.The proposed alternative is the combination of many more propagation mechanisms, however, this approach generates a high computational cost.The simulations show that the non-hit points of the Non-Modified and the Modified RT are different, which is explained by the slightly different path that the inclusion of atmospheric refractivity effect produces. Table II summarizes the statistics.It can be perceived that a Modified RT algorithm shows a encouraging result in statistical comparison.However, in order to consolidate a conclusion about the results obtained with the proposed algorithm, it is necessary to analyze its performance and behavior in another scenario and at a different frequency band.both situations are illustrated, Figs.14(a) and 14(b), respectively.It can be noted that unlike the first realistic case, in this environment a complete pathloss curve is obtained along the terrain, this is mainly due to the choice of higher receiver points.This means that for all the observation points considered over the Pacific region terrain profile, the proposed multipath approach achieves a combination of ray paths that hit a known receiver point.
Table III displays the statistics.Similarly to the realistic case above, it can be observed that a Modified RT algorithm shows a slightly better result in the statistical comparison.
The results obtained in both realistic cases allow to conclude that the proposal of atmospheric refractivity inclusion into a multipath approach, based on RT solutions, is a promising propagation model to be applied in this kind of environments.Additionally, the capacity of the proposed algorithm to estimate coverage and characterize radio channels at sub-6 GHz frequency bands (projected for 5G systems operation) is also demonstrated.
Many two-dimensional (2-D) space implementations can be found in the literature to solve radiopropagation problems, where numerical techniques such as Parabolic Equation (PE) model, Finite Difference (FD), Method of Moments (MoM), Finite Element Method (FEM) and Finite-Difference Time-Domain (FDTD) method have been successfully implemented in closed and open environments Journal of Microwaves, Optoelectronics and Electromagnetic Applications, Vol.22, No. 2, June 2023 DOI: http://dx.doi.org/10.1590/2179-10742023v22i2272846300

Fig. 3 .
Fig. 3. Incident and reflected ray paths over a sloping terrain. tan

Fig. 4 .
Fig. 4. Flowchart for a modified multipath model based on RT techniques.

Fig. 8 .
Fig. 8. Ray paths to calculate Pathloss results and PDP by using a Modified RT algorithm.
) illustrates the comparison of the pathloss horizontal profile results obtained with the Non-Modified RT and DMFT-SSPE tecniques.Similarly, Fig. 11(b) shows the pathloss curves results obtained with the proposed Modified RT algorithm and DMFT-SSPE solution.

Fig. 12 .
Fig. 12. Ray paths conformation at the final position of the radio link for Colombian Pacific region case.

Fig. 13 .
Fig. 13.Power Delay Profile at the final position of the radio link for Colombian Pacific region case.

TABLE I .
STATISTIC COMPARISON FOR PATHLOSS RESULTS IN MIXED SCENARIO CASE

TABLE III .
STATISTIC COMPARISON FOR PATHLOSS RESULTS IN PACIFIC REGION CASE RT Technique Mean Absolute Difference (dB) Standard Deviation (dB)