I. INTRODUCTION

Modern communication and radar systems need effective tools to foresee propagation. In many cellular communication and HF and VHF radar sites, high precision propagation simulation is vital. The results are necessary for choosing proper places for ground stations and optimum height for transmitters and receivers, and also reducing the costs of implementing communication systems. One of the primary aims of propagation simulation is to forecast path loss (*L _{P}*) with high precision. Path loss is defined as:

Where *P _{T},G_{T},G_{R},P_{R},L_{P},* and

*L*denote transmitted power, transmitter antenna gain, receiver antenna gain, propagation loss, and system loss, respectively. In (1), antenna gains are measured with respect to an isotropic reference antenna, and

_{S}*L*represents all system losses, including thermal loss, signal loss in cables, mismatch loss, etc. For example, this equation is used in cellular communication system design to foresee each antenna coverage and dead zones before they are installed. Propagation loss denotes all loss caused by factors such as propagation medium refraction index, multipath interference, and reflections from natural and man-made obstacles. So the simulator must take all such parameters into account.

_{S}In the past two decades, researchers have been in favor of path loss calculation, with numerous efforts made in this direction. In [^{3}], virtual tools for ground-wave propagation modeling and numerical simulation strategies are reviewed. The propagation loss of 71 floors in 17 different buildings, in an urban area, is measured in 4 frequencies ranging from 800 MHz to 8 GHz in [^{4}]. In [^{5}], a study of path loss models in the 2 – 73.5 GHz range in an Urban-Macro Environment is performed. In [^{6}], electromagnetic (EM) wave measurement in 2.4 GHz in underwater environment is conducted. Also, it is shown that the propagation velocity and absorption coefficient in fresh water are frequency independent.

So far, several numerical methods have been introduced for the analysis of electromagnetic problems. Most of these methods are based on full-wave analysis in the time or frequency domain; and utilize one or two difference equations, plus boundary conditions (depending on the specific case) to solve the problem [^{7} - ^{10}]. Each method is used to solve a specific class of problems. Although these methods are powerful and can be used to solve a wide variety of problems, due to improper usage and intrinsic limitations, sometimes the derived answer is erroneous. Unless errors such as overflow or underflow have occurred, computers always produce answers. The case is whether this answer correctly represents the physics of the problem in the real world.

One of the most ubiquitous methods is the finite-difference time domain, or FDTD method for short. This method is performed in the time domain and has been used for decades. Although real wave propagation takes place in three-dimensional space, a two-dimensional FDTD approximation can be used, if the problem is symmetric relative to one of the coordinate axes.

In this paper, the propagation factors derived using the TDWP and FDTD methods are compared. In order to find differences in the results obtained from the two methods and their causes, spontaneous fields in different time periods might be studied. In the first part, the TDWP method and its application in solving problems are explained. In the next part, a simulation in the whole space is performed using both TWDP and FDTD, and the results are compared. Finally, the simulation process is modified to reduce these differences.

II. TDWP METHOD

The wave propagation problem may be set in city streets with numerous buildings, in a vast desert, or in a mountainous area. The classical FDTD method faces shortcomings when simulating wave propagation along long paths. Suppose we would like to study wave front propagation with spatial Gaussian distribution (along the height axis) in a 10 km path and a height of 50 m. If we divide this space into a mesh of 10 cm cubic cells, there would be 5000000 cells in this space, not considering PML blocks located at the left, right, and top of the computational space. The propagated wave consists of three components, regardless of its polarization. Every time these three components are updated, the new values should be stored in matrices. So we need three matrices with 5 milli on elements each, all with 15 million elements updated in every time step. Considering the fact that the utilized program uses some of the computer's RAM and CPU cycles, it can be concluded that processing such huge amounts of data with personal computers is impossible. Even if sufficient RAM capacity is obtained, processing the data takes considerable time.

To overcome this problem, instead of updating all components for every computational cell at every moment, a rectangular window around the pulse can be presumed, containing most of the propagated pulse energy. In this manner, the computational space can be divided into smaller windows, as shown in fig. 1 and 2.

These smaller rectangular windows must be surrounded by PML from the left, right, and top. PML blocks surrounding the computational space prevent reflection and refraction of the waves along the boundaries. The aim is to compute propagation loss along the wave's path and also along the height axis. To this end, we trigger the computational space with a short duration pulse. The trigger source is located at the beginning of the first window. When it begins to stimulate field components in the adjacent cells, the wave front begins to move along the length of the window to its border. When the wave front first reaches the last column of cells in this window, the magnitude of field components in this column are stored in some part of the computer memory, until the effective pulse energy reaches the end of the current window. After the resulting data is stored, the first window is closed and a new computational space is defined as the second window, containing the physical properties of the next compartment. In the new window, all cell field components are initially zero, which means that none of the cells are stimulated. The field components of cells in the last column of the preceding window are immediately copied as field components of cells of the first column of the current window (fig.3). Similar to the first window, the pulse travels to the second window, and the same procedure is repeated. By continuing in the same manner for numerous windows, pulse flight can be simulated for paths as long as tens of km.

In fact, in this manner, a rectangular computational space with very great length is divided into rectangular computational spaces with smaller length, but the same height as the huge rectangle, so that the simulation can be readily done with a personal computer.

III. Problem formulation

A. Wave Propagation Formulation

Although the physics of propagation is taking place in a three-dimensional space, for a symmetrical problem with respect to one of the Cartesian coordinates, we can use a two-dimensional FDTD approximation. Maxwell's equations and the constitutive parameters, in the form of vector differentials, in an isotropic space containing magnetic and dielectric materials, are as follows:

Where
*σ* And *σ ^{*}* are the electric and magnetic conductivity of the material, and

*ρ and ρ*are the electric and magnetic charge, respectively. By assuming

^{*}*TM*(

_{y}*H*,

_{y}*E*,

_{x}*E*) polarization for the propagated wave and considering a zero value for the electric and magnetic conductivity coefficient, equation 2 and 3 can be simplified as:

_{z}In Fig. 4, the configuration of the electric and magnetic field nodes are shown for *TM* polarization. In this figure, fields components within closed dash lines have the same index in the implementation of the algorithm. The nodes of the electric field, with lines to illustrate the direction of field components, are shown with solid circles; and the nodes of the magnetic field are shown with empty squares. Triangle shaped dashed lines represent a group of nodes that have identical array indices in the computer algorithm. For example, in the bottom left of the grid, all three nodes inside the triangle have (*i* = 1, *k* = 1) indices. The diagram on the bottom left of the grid shows nodes related to each field and their indices for a typical group. In the diagram on the right side of the grid, the same nodes and their array indices have been shown when used in a computer program. In this polarization, none of the nodes of the fields were located in places with the correct values of discretized space.

With the above explanations, and assuming a discrete time form of

The components of the electric and magnetic field are calculated in the Δ*t*, 2Δ*t*, 3Δ*t*, …, *n*Δ*t*, and

If a node of the magnetic field is placed within the setting of the conductor, it is assumed that the entire cell containing that node is placed inside the conductor. Thus, all electric field components of that cell (which include that node of the magnetic field) are fixed to zero.

IV. Path loss calculation

In this problem, the propagation space is a semi-confined space, limited by *z* = 0 and *z* → ∞ vertically, and by *x* → ∞ and *x* → ∞ in the plane parallel to the ground. We would like to calculate the loss caused by the propagation medium, without considering the effect of the receiver and transmitter antenna. In other words, we would like to calculate the power lost in propagation, after the wave is transmitted by the transmitter antenna, before it is absorbed by the receiving antenna. Therefore, in (10), *G _{T}, G_{R},* and

*L*should be substituted with 1. In every given point,

_{s}*P*can be calculated in terms of field magnitude [

_{R}^{2}]:

Therefore, to compute the propagation loss, we need to calculate the electric field magnitude at every given point in the propagation space.

To be able to compare the results for different media, we need to eliminate the effect of the reduction of wave power due to the distance from the source via the inverse square law, which is the only loss associated with a lossless medium without obstacles. The propagation factor at every given point is defined as |*E/E*_{0}| (for any specific frequency). Where *E* and *E*_{0} are the measured and reference electric field at the desired point, respectively [^{2}]. The reference field is the field propagated from the same source, not taking the terrain into account (fig.5). When the wave reaches the boundaries of the propagation space into free space, it must be ensured that it leaves the studied medium and will not be reflected back. To simulate the free space surrounding the propagation space, electromagnetic absorber layers must be used. These layers (fig.5) attenuate the wave, so that it doesn't return to the propagation medium after leaving it. one of the most appropriate types of such absorber layers, for simulations, is the “Perfectly Matched Layer” or PML for short [^{11}]. So in this method, we first perform the simulation for the terrain and store the derived *E* (*t*) values for both a specific height and different distances, or also a specific distance and different heights. In this step, the propagation medium is surrounded by the PML from the top, left, and right. In the next step, the simulation is repeated with the terrain removed and the propagation space covered by PML from the bottom as well. In this step, the *E*_{0}(*t*) values are stored in the same cells as the last simulation [^{12}].

Finally, after both steps are done, by performing DFT analysis on the time-domain data available (in the desired bandwidth), |*E/E*_{0}| for each desired frequency is derived. Usually, it would be favorable to calculate and plot this ratio in decibel units (2*10*log*|*E/E*_{0}|) for a specific distance or height.

A. Source Stimulation

To stimulate the fields of the computational space, we need to model a physical source. Because finite time distributions, which also have limited frequency bandwidth, provide better responses in FDTD simulations, the pulse shape (in the time domain) is usually selected as Gaussian. However, here we prefer to use the derivative of the Gaussian function, in order to prevent DC effects that can cause shifting of the wave propagation frequency in some applications. However, we should note that the desired pulse for our source should cover the range of frequencies determined by the bandwidth. The Gaussian function derivative used in our algorithm is as follows:

Where *ϕ* is the time shift of the function toward the positive direction of the time axis; which is equal to
*BW* = 200 *MHz*.

Now which *i* and *j* the simulation applies to depend on the type of stimulation. The distribution, e.g. pulse time distribution, can also be in different forms. One of the most popular spatial distributions along the vertical axis is the Gaussian distribution; where its centricity is shown with *h _ source*.

The statistical process of spatially distributing the pulse along the *z* axis is achieved by its variance.

*h*_*source* is the height of the spatial distribution center for the source from the ground and *σ _{ν}* is the spatial variance of the source. In fig.7, the spatial distribution of the source with centricity of

*h_source*= 25

*m*from the ground at a height of 50 m is shown for different variances.

Thus, the pulse function *pulse*(*y, t*) will be the multiplication of the two functions:

In order to compare the FDTD and TDWP methods, a simulation is set up as shown in fig. 8, containing two conducting triangles in free space. The ground plane is also assumed to be a perfect electrical conductor. The dimensions of the simulated space are shown in fig.8.

The fields used in this simulation are Gaussian derivative functions of time with a central frequency of 50 MHz and a bandwidth of 100 MHz. The spatial distribution is a Gaussian derivative function of height as well. The center of the spatial distribution is taken to be 15 m high. The stimulation source is assumed to consist of a number of dipoles perpendicular to the ground. Therefore, the polarization is TM along the x direction. We would like to obtain the propagation factor as a function of height for four different distances, as shown in fig. 8. We compute the desired factors by two different procedures. In the first procedure, the simulation is performed for the whole space (using FDTD method) as shown in fig. 9 and the results are recorded.

In the second procedure, the space depicted in fig.8 is cut into two halves, or in fact two windows (TDWP method), and the simulation is rerun. To achieve a better understanding of the propagation, the vertical component of the normalized electric field in the computational space is color-plotted. The instantaneous images of the wave front, yielded from the FDTD and TDWP methods, are shown in fig.9.

Propagation factor curves derived from these two methods are shown in fig.10. It can be seen that the two curves do not match perfectly. The reason can be deduced by analyzing the images in fig. 9.

As it can be seen in frames 3 and 4 of fig.9 (in the TDWP method), when the propagated wave reaches the end of the first window, it is restored in the beginning of the second window. But, as seen in frames 6 of fig. 9, reflections from the right side triangle (situated in the second window) are neither propagated back into the first window, nor absorbed by the PML located at the left boundary of the second window. Instead, the reflection is re-reflected into the second window.

The reason the reflected wave from the second triangle is not reflected back into the first window is obvious. As consequent computational spaces are discrete, i.e. not continuously connected as in the FDTD method, when the wave front passes through a window for which calculations are over, no field will be regenerated in the previous windows. It is apparent that should we apply the method we use to store fields yielded from previous window and regenerate the field in current window in reverse direction, the result would be a heavy computational burden as FDTD. The waves reflected from the second triangle are absorbed by the PML blocks located at the right and top borders of the window; however, they are reflected by the left border. The reason is that the first column in the window is regenerating fields from the last window during the whole duration of the simulation. In other words, these cells do not take part in the wave propagation that happens in the right to left direction and solely represent the values stored from the previous window. It may seem that when the reflected wave from the second triangle reaches this column, the column's cells are not regenerating any field. But, in fact, the visible portion of the wave generated by these cells is only the tangible part, and these cells continue to regenerate fields, even though the resulting magnitudes are too small to be normally noticed. It is obvious then that as these cells do not take part in wave propagation, they practically act as a thin conducting wall and cause the reflected wave front from the second triangle to be re-reflected back into the second window. This is actually also the case for the cells generating the stimulation. As in this particular example it is possible to predict when the effective portion of the wave has passed through the previous window to the current window, the latter problem can be tackled by limiting the regeneration of field components in the first column. Moreover, considering the fact that the time domain Gaussian derivative function is known, the same trick can be used to limit the duration of the stimulation in the first window. Another trick which can be used for the stimulation is to set some threshold, so that when the source field magnitude falls below this threshold, the cells generating the source fields are set to zero. Fig 10, shows the comparison between the propagation loss results of two method (TDWP vs FDTD) in several lines.

B. Limiting field regeneration in the first column of the second window cells

In fig.11, instantaneous images taken from TDWP simulations with and without regeneration limit for the first column of the second window cells are compared.

Frames 1 to 6 demonstrate snapshots taken from the 608^{th} through 928^{th} time steps of the simulation, with 64 time steps between each two consecutive frames.

If the propagation factor computed for lines *x*_{3} and *x*_{4} using the modified simulation method is compared with FDTD, it can be observed that the modified TDWP method provides results closer to the results derived from FDTD [^{5}]. Fig 12, shows the comparison between the propagation loss results of two method, Optimized TDWP and FDTD, in Lines 3 and 4.

As it can be sometimes impossible to predict the time when wave reflection from the terrain reaches the end of the previous window, field regeneration duration cannot be limited in general. This problem can be solved by applying the substitution
*F* can be any field component *H _{x}, H_{y}*, or

*H*and

_{z},*f*(

*t*) represents the stored value of the component from the last stage. It might seem that the aforementioned substitution causes new values to sum up with old values; but it must be borne in mind that field regeneration takes place when the field component values are being updated, so that every time a value is updated, it does not contain the previous value anymore.

C. Analysis of Improving Factors of Optimized TDWP Method

In following some relations are defined which can help us to make a suitable analysis regarding the efficiency of proposed Technics, in order to improve traditional TDWP method.

If we denote the propagation loss factor of Traditional and our Optimized TDWP method as *PF _{TDWP}*

_{(op)}and

*PF*

_{TDWP}_{(re)}, respectively, the amount of generated Error in terms of y (axis) will be calculated for Traditional and Optimized methods as below, respectively. It should be emphasized that the one-window FDTD method is not applicable for long pathes. It is applied to this study just for a path of 100 m.

The ratio of *δ _{op}*(

*y*) to

*δ*(

_{re}*y*) is defined as a Improving Factor for each height of y.

For investigating the improving factor of proposed method in one given line (like 0 < *y* < *H*, *x* = *x*_{o}), we should first calculate the below Integrals for each line.

The ratio of Δ* _{re}* to Δ

*is defined as a Improving Factor for each vertical line x.*

_{op}In similar way, for calculating the propagation loss factor along a horizontal line (like 0 < *x* < *L*, *y = y*_{o}) we will have :

The ratio of *δ _{op}*(

*x*) to

*δ*(

_{re}*x*) can be defined as improving factor for each length of x.

For investigating the improving factor of proposed method in one given line (like 0 < *y* < *H*, *y* = *y*_{o}), we should first calculate the below Integrals for each line.

The ratio of Δ* _{op}* to Δ

*is defined as a Improving Factor for each vertical line x.*

_{re}At last, for attaining to a total improving factor in one window, relation 16 (or 20) should be employed over a given surface like
*x*_{o} = 52.5*m* (*Line*3), the amounts of *δ _{ορ}*(

*x*) and

*δ*(

_{re}*x*) are 1.26 and 3.55 respectively in 19 m; Thus

*cf*(19) = 0.35. by calculating the improving factor for

*x*

_{o}= 52.5

*m*in all heights we will have

*CF*(

*x*= 52.5) = 0.16, While this factor will be

*CF*(

*x*= 70) = 0.05 for

*x*

_{o}= 70

*m*(

*Line*4). It shows that this technic make a better operation for

*Line*4 in comparison to

*Line*3.

V. Conclusion

In this paper, the TDWP method for wave propagation simulation and propagation loss computation was introduced. The benefits of this method in comparison to the FDTD method where discussed, and by comparing these two methods for a terrain with two simulation windows, the shortcomings of TDWP methods where demonstrated. After that, by proposing the approach of limiting field regeneration for first column cells of each window, the matching between the propagation factors curves derived from the two methods was ameliorated. The TDWP method enables wave propagation simulation with personal computers; but provides better results for plain surfaces rather than for terrains. Such as the simulation of wave propagation for islands located near coasts where surface impedance changes as a function of distance from the source. In order to achieve more accurate results for surfaces with terrains, modifications should be made, as discussed above.