The Evaluation of a Novel Asymptotic Solution to the Sommerfeld Radiation Problem using an Efficient Method for the Calculation of Sommerfeld Integrals in the Spectral Domain

In this work, a recently developed novel solution of the famous ”Sommerfeld Radiation Problem” is revisited. The solution is based on an analysis performed entirely in the spectral domain, through which a compact asymptotic formula describes the behavior of the EM field, which emanates from a vertical Hertzian radiating dipole, located above flat, lossy ground. The paper is divided into two parts. First, we demonstrate an efficient technique for the accurate numeric calculation of the well– known Sommerfeld integrals, required for the evaluation of the field. The results are compared against alternative calculation approaches and validated with the corresponding Norton figures for the Surface Wave. Then, in the second part, we briefly introduce the asymptotic solution of interest and investigate its performance; we contrast the solution versus the accurate numerical evaluation for the total received EM field and also with a more basic asymptotic solution to the given problem, obtained via the application of the Stationary Phase Method (SPM). Simulations for various frequencies, distances, altitudes and ground characteristics are illustrated and inferences for the applicability of the solution are made. Finally, special cases, leading to analytic field expressions, close as well as far from the interface, are examined.

tegral representations for the received EM field, by means of a generalized solution to the respective Maxwell equations boundary value problem. Working in this domain has the advantage that no Hertz potentials and their subsequent differentiation are required for the evaluation of the fields.
In [18], the Stationary Phase Method (SPM) [19]- [21] was applied to the general integral expressions for the EM field and the well-known analytic formulas for the Space Wave, defined as the complex interference of the Line of Sight (LOS) field and a portion of the field emanating from the dipole's image point (also called the Reflected Field), were obtained as the high frequency asymptotic solution to the complete problem. In [22]- [24] we focused on the numerical evaluation of the field's integral expressions and how they compare with the respective high freq. approximation ones. It was revealed that accurately evaluating the Sommerfeld integrals in the spectral domain is also not a trivial task. The result is sensitive on the position of the singular points in relation to the integration path, an issue that has been also a major problem and a matter of debate in various related research works [6].
Then in [25], the mathematical formulation of the problem in the spectral domain was redefined for the usual case where σ ωε 0 , ie for a highly conductive interface, which is the case for most practical frequencies of interest in terrestrial communications. As shown there, a special contour integral, called "Etalon Integral", was used to deform the original contour of integration, associated with the Sommerfeld Integrals in the spectral domain, through the application of the Saddle Point Method (SDP). The above mentioned, "Etalon Integral" can be expressed in terms of Fresnel Integrals and has interesting properties, which can reduce the problem related to the vicinity of the saddle point to the pole point [26]- [31]. The result was a compact asymptotic solution that better expresses the variation of the field in the high frequency regime. Moreover, using the small and large argument approximations, associated with the Fresnel Integrals [32], pure analytic expressions were extracted, describing the behavior of the EM field close, as well as far away from the ground interface.

B. Scope of this Research
The analysis in [25] was constrained to the pure mathematical formulation of the problem. Simulations and related Fig. 1. Hertzian Dipole above infinite, planar interface. Point A is the image of the source A with respect to the ground (yz-plane), r 1 is the distance between the source and the observation point, r 2 is the distance between the image of the source and the observation point, θ 2 is the "angle of incidence" at the so-called "specular point", which is the point of intersection of the ground (yz-plane) with the line connecting the image point and the observation point, and finally, ϕ = π/2 − θ 2 is the so-called "grazing angle".
figures are still pending for validating the method. Hence, in this work extensive simulations are demonstrated and the method's efficiency is examined, with respect to the numerical evaluation of the respective integral formulas for the EM field. Moreover, we compare the method's performance against the more basic asymptotic solution of [18], which is based on the application of the SPM method. However, as mentioned above [23], [24], the accurate evaluation of the Sommerfeld integral expressions is not a straightforward task and this is due to both the presence of singularities along the integration path as well as the particular complex nature of the integrands. For that reason various specialized commercial software have been used for obtaining adequate results, for example the AWAS tool used by Sarkar et al. in [6].
In this paper we show that using an appropriate variable transformation, it is possible to convert the generalized integrals of [18] into fast converging formulas, which are rather suitable for numerical calculation, using standard numerical integration techniques. Particularly, the integral expression describing the received EM field, is broken down into two terms; one relatively easily computed definite integral, of finite integation range and another integral of semi-infinite range. However, the latter integrand proves to be a fast decaying exponential function, resulting in very fast convergence times. Comparisons against the numerical results published in [23], [24] demonstrate the advantage of the method. Also a validation against Norton's figures, associated with the well-known surface wave [14], [15], is exhibited. Preliminary results have already been demonstrated in [33].
Then, using the numerical integration results as the baseline, we juxtapose the novel asymptotic solution of [25] with the SPM based method of [18]. Moreover, an investigation of the above mentioned analytic formulas, which according to the analysis in [25], should reflect the field behavior close to as well as far away from the ground level, is made. As dictated by our simulations, the near ground-level predictions, for a sliding angle of incidence, are not validated.

C. Problem Geometry
The geometry of the problem, as given in [25] and also briefly described here for ease of reference, is shown in Fig. 1. A vertical small (Hertzian) dipole, characterized by dipole momentṗ = p ·ê x , p=const, is directed to the positive x axis, at altitude x 0 above infinite, flat and lossy ground. The dipole radiates time-harmonic electromagnetic (EM) waves at angular frequency ω = 2πf (e −iωt time dependence is assumed). The relative complex permittivity of the ground is: ε r = ε /ε 0 = ε r + iσ/ωε 0 , where σ is the ground conductivity, f the carrier frequency and ε 0 = 8.854 × 10 −12 F/m is the permitivity in vacuum or air. The goal is to evaluate the received EM field at an arbitrary observation point above the ground level, namely at point (x,y,z), shown in Fig. 1.

D. Structure of the Article
In what follows, Section II recaps the fundamentals expressions for the EM field in the spectral domain, the issues associated with their numerical calculation and demonstrates how a simple variable transformation may lead to fast converging integral formulas, suitable for evaluation in the computer. Through various simulation results, we illustrate the advantages and validate the accuracy of the redefined expressions. Then in Section III, we give a brief overview of the asymptotic solution of [25] and through an extended set of simulations -comparisons we demonstrate its efficiency. Moreover, a discussion regarding the applicability of the closed-form formulas, also predicted by [25], is given. Finally, in Section IV we summarize on the major findings and we make a brief discussion on potential extensions. The whole analysis is given for the electric field. Expressions for the magnetic field are derived similarly or by suitable use of the duality principle.

II. EFFICIENT FORMULATION FOR THE EM FIELD INTEGRAL EXPRESSIONS IN THE SPECTRAL DOMAIN A. Original Integral Expressions
According to the analysis of [18], performed in the spectral domain, the electric field at the observation point of Fig. 1 is given by the following integral expression, where E LOS denotes the direct or LOS field, E R is for the field scattered by the flat and lossy ground and the vector functions f 1 (k ρ ) and f 2 (k ρ ) are given by with H (1) 0 being the Hankel function of zero order and first kind and k 01 , k 02 , the wavenumbers of propagation in the air and lossy medium (ground) respectively.
Expressions (1) -(4) expose the following difficulties when coming to the evaluation of the respective integral through common Numerical Integration (NI) techniques: -The range of integration extends from −∞ to +∞, resulting in potential computational errors for large evaluation arguments. -The Hankel function, H 0 , exhibits a singularity at k ρ = 0 and although it is proved that this is a logarithmic singularity [34] and does not break the integrals convergence 1 , it can affect the accuracy of the numerical integration results, when implemented in the computer.
-In addition, it is obvious that k ρ = ±k 01 are also isolated singularities of (2), (3) and despite they are still integrable singularities [34] 2 , a sufficient small range around those points must be excluded, when numerically evaluating (1) in the computer. As argued in [24], doing so may severely affect the accuracy of the results. Of course, the above mentioned accuracy issues, are of practical importance, only as far as the Scattered Field, E R , is concerned, for which no analytic formula exists. For the LOS field, a closed-form expression does exist as following, which is found by solving the problem of an isolated hertzian dipole source in free space [35]. However, for verification purposes, in the sections that follow, we will also briefly examine the integral representation for the LOS field as well. Note that (5) reflects the exact solution of the problem, encompassing both the near field and far field components and is expressed in the cylindrical coordinate system, as is the case for the whole analysis herein.

B. Reformulated Integral Expressions for the EM Field
We now focus on the scattered field, i.e. the second integral expression of (1), which may be written as 0 (kρρ) = 0 [23] 2 it's a square root integrable singularity that applies to Rule1 of [34] Starting with (7a), we perform a simple variable transform, k ρ = k 01 sin ξ, which apparently maps the [−k 01 , +k 01 ] range to [−π/2, +π/2]. With this transform, (4) is translated to Ultimately and if we also take into consideration the definition for f 1 , as given by (2), the expression for I 1 becomes or equivalently, it may be written as We may further elaborate on (11), if we make use of the following properties for the Hankel function [36], [37], namely, (with the latter implying an analytic continuation of H 0 in the upper half plane) and also observe from (10) that the reflection coefficient R (ξ) is an even function, with respect to ξ. Overall, we get where J 0 denotes the zero order Bessel function. For the integrals I 2 and I 3 we follow a similar approach. Particularly, in (7b) we apply the variable transform k ρ = k 01 cosh ξ, while in (7c) we set k ρ = −k 01 cosh ξ. In both cases, the original ranges of integration, [−∞, −k 01 ] and [k 01 , +∞], in the k ρ domain, are mapped to [0, +∞] in the domain of ξ. Moreover, (4) becomes Performing the necessary calculations and also using (12), (13), we may combine the results for I 2 and I 3 as where I 23 = I 2 + I 3 and the reflection coef. R is given by Substituting (14) and (16) to (6), we reach to an integral formula for the scattered field, E R , suitable for numerical calculations. With a similar process for the first integral of (1), we get the equivalent expression for the LOS field. Overall, the redefined integral expressions for the direct and scattered fields are given by An inspection of (18) and (19) might yield useful insights, which are mentioned here, since they have not been clarified in [25]. Both formulas, express the field as a complex superimposition of plane waves. Equation (19) expresses the direct field as an integral expression over the dummy variable ξ, which is an auxiliary, transformed variable of the spectral domain coordinate k ρ . As required by the problem's geometry, the field is cylindrically symmetrical (no φ -component) and it is expressed as a complex summation of contributions, originating from the dipole's location, hence the dependence of the field on the horizontal distance, ρ and the relative height difference x − x 0 . Moreover, it is easy to identify that the x -component of the field is symmetrical, while the ρ -component is antisymetrical above and below the dipole's position, in accordance to the conventional solution of the dipole's problem [35]. The expression for the scattered field has a similar form and can be considered as the integral generalization of Fresnel's theory, due to the existence of R (ξ) and R (ξ) in (19), acting as reflection coefficients, whose values depend on the ground characteristics. Also, the field depends on the cummulative distance x + x 0 , as if the source is located at the image point A of Fig. 1.
Equations (18), (19) remedy the accuracy issues, mentioned in section II-A, above: -They utilize the zero order Bessel function, J 0 , instead of H 0 used in (1), which is a smooth, finite special function with no singularity, whatsoever.
-The singularities at points k ρ = ±k 01 have also been removed. Hence, no need to exclude any range around them is required, when using any kind of numerical integration technique, in order to calculate (18) and (19). -The result is expressed as the sum of two integrals, one bound definite integral, in the range [0, π/2] and a second improper integral, for which the range of integration extends from 0 to ∞. However, due to the presence of e −k01(x+x0) sinh ξ , the second integrand is a fast decaying function, practically making the integral a bound limits one that is fast converging, easily evaluated in the computer. The above findings are also visible in the simulations that follow.

C. Simulations Results and Comparisons
The parameters for the various simulations (i.e. transmitter -receiver heights, ground parameters, operating freq. etc) are indicated within the figures and were selected such that a comparison with preceding, referenced results of [18], [23], [24] is possible, if applicable. Fig. 2 exhibits the numerical evaluation (NI) for the scattered electric field, E R , using the redefined integral expression (19). It is compared to the equivalent values, obtained using the initial integral formulas for the electric field 3 , introduced in [16], [18] that is by using (1) -(3), above. Along with the NI results, the high frequency approximation values, are also superimposed. These values were obtained as in [23], [24], i.e. through the application of the SPM method on (1). SPM is a useful asymptotic technique for the evaluation of complex integrals, particularly when the integrands expose rapidly changing phase components 4 .
As deduced in [24], SPM results are expected to be accurate in the far field, i.e. at least at distances over 10 -15 wavelengths, or above 100 -150m, for the 30MHz case, shown in Fig. 2. Therefore, using the SPM data as the baseline, it is obvious that only the numerical evaluation of (19) achieves the required accuracy and this is noticeably evident for distances larger than the characteristic distance of the so-called Pseudo -Brewster angle, defined as the angle of incidence, θ B , where the reflected field is minimized [35] 5 . On the contrary, the numerical computation of (1) -(3) fails to describe the electric field behavior, which may be attributed to the reasons analyzed in Section (II-A), above.
In Fig. 3 we demonstrate various field types and components for the case of a hertzian dipole radiating at 300 KHz, which is regarded as the frontier between the Low Frequency (LF) and Medium Frequency (MF) bands [35]. For the LOS field we used (5), while the space wave was evaluated as in [6], i.e. by using the concept of the Fresnel Reflection Coefficient for the reflected field. The scattered field 6 was numerically computed via (19).
Due to the small antenna heights and the long distances involved (10 -20 km), the space wave is expected to diminish [35]. Therefore, the link is established primarily by means of the Surface Wave, defined as the remaining field, after subtracting the geometrical optics field (or space wave) from the complete or Total Field [38]. This is actually verified in the top plot of Fig. 3, with the total field curve being very close to the surface wave results. As a confirmation of the validity of the results, our surface wave calculations are compared with Norton formulas [14]. The respective curves are almost identical.
The bottom half of Fig. 3 displays the behavior of the integrand associated to the second integral of (19), i.e. the generalized integral over the [0, ∞) range. Actually, we are dealing only with the x-component 7 of this integrand, de- 6 The terms scattered field and reflected field are not equivalent [6]. 7 This is the major field component for the considered problem [6]. noted as function g ex (ξ) in Fig. 3 (the behavior for the ρcomponent is similar). The integrand is confined in a small window of the integration variable, ξ, outside of which and especially for large values of ξ, it is practically equal to zero. This is an outcome of the fact that the exponential function e −k01(x+x0) sinh ξ decreases much faster than the increase rate of cosh 3 ξ 8 . The bottom line is that the second integral of (19) essentially becomes a bound limits definite integral, easy to accurately evaluate in the computer, using common numerical integration techniques, which gives our formulation a computational advantage.
In Fig. 3 it is also interesting to notice the fluctuating behavior of g ex (ξ). This is an outcome of the oscillating nature of the Bessel function J 0 . Its effect on g ex (ξ) is apparent by observing the bold line of Fig. 3 (g' ex (ξ) in the figure), which demonstrates how the integrand would behave, if it hadn't been for J 0 . Again, the confinement of the integrant within a "narrow-band" of the variable ξ is apparent. It also seems that g' ex (ξ) acts like a slightly-shifted envelope function of g ex (ξ). However, notice that this is a normalized illustration of g' ex (ξ), to the respective magnitude of g ex (ξ), since the order of magnitude between the two is totally different.
The simulations of Fig. 3 are now repeated for a high frequency scenario in the VHF / UHF band. The source and observation points are located even closer to the ground level, in an attempt to detect meaningful surface wave values, if possible, in this higher frequency case. Nevertheless, as illustrated in Fig. 4, this is a situation where the space wave almost completely dictates the field behavior. The pursued surface wave becomes very quickly negligible and this is actually in accordance with Norton's predictions, where the large values for the so-called Arithmetic Distance, results in very small values for the attenuation coefficient, hence small surface wave figures in the high frequency regime [35]. These results are also a validation of the SPM method, which as mentioned in Section I, it emerges as the asymptotic solution for the complete problem, in the high frequency case.
Finally, notice in the bottom graph of Fig. 4, how quickly g ex (ξ) vanishes (in this case the real part is shown), making thus the convergence of (19) very fast. Moreover, due to the alternating positive and negative values, it is expected that the effect of g ex (ξ) on the overall result will be insignificant. The same arguments hold for the ρ-component of (19), justifying the small observed values, as far as the surface wave field is concerned. Put it differently, for the case shown in Fig.  4, the major contribution in (19) comes from a narrow area around the Stationary Point, which in this problem lies within the [−π/2, +π/2] range [24]. This contribution yields the reflected field in an asymptotic sense, as first shown in [18] with the application of the SPM method. In the rest of the integration range, the integrand is related with the surface wave [38] and exposes a behavior similar to Fig. 4, thus having minimum impact to the final result. This was a major assumption for the application of the SPM method in [18], which is now numerically validated in this high freq. scenario.
As a last validation, in Fig. 5 we demonstrate various field components for the exact scenario, illustrated in Fig. 4 of [18]. The simulation parameters are as those of Fig. 2, except for the horizontal distance range. In [18] only the Norton's surface wave was evaluated, whereas here we also compare with the NI results. Moreover, we perform a comparison between the analytic expression for the LOS field and its equivalent integral form ("LOS field NI" in Fig. 5), as both given in Section II by (5) and (18) respectively. Again, our numerical evaluation for the surface wave is more or less identical with Norton's values. No needless to say that we also achieve a perfect match between (5) and (18), essentially meaning that our redefined integral formulation for the EM field, described in Section II-B, is effective and accurate. In other words, the drawbacks, associated with the original integral expressions of Section II-A, do seem to have been mitigated.  Convergence times in milliseconds (ms). The horizontal distance was set to ρ = 1 km. The rest of the problem parameters: X, X 0 , σ, εr, I, 2h, were set as in Fig. 5. Simulations performed on a 64bit, Quad Core CPU@2.60 GHz, 16MB RAM platform, using MatLab.
We close this section with a few comments regarding the method's efficiency. The convergence time of the method depends on four key aspects; a) the utilized HW and SW platform, b) the selected NI algorithm for the calculation of (18), (19), c) the required error tolerance and d) the problem parameters; especially the frequency of operation, for a given Transmitter -Receiver (T-R) distance and altitude, or on the electric distance k 01 r, when their combined effect is accommodated. The first three factors seem quite reasonable. Regarding the fourth one, which in first glance may seem less relevant, Fig. 6 provides a good reasoning. It displays the behavior of the first integrand of (19) 9 , shown as h eρ (ξ) in the figure. It is obvious that higher frequencies contribute to additional oscillations and this is not surprising if one observes that h eρ (ξ) includes a Bessel and a phase function (cosine or sine function when the real or imaginary part is considered), which are increasingly fluctuating for larger arguments, as occurs in this case, when the frequency f = k01·c 2π increases. Hence, one might expect that more steps or intervals are required, for the NI algorithm to achieve a given error threshold. Table I demonstrates the measured performance of our method, at various frequencies, utilizing two widely used NI techniques for the evaluation of (19), namely the Adaptive Simpsons and the Trapezoid method [39]. We are able to calculate the fields at almost an arbitrary accuracy level and at very reasonable computational times. 10 Table I, also exposes the effectiveness of adaptive quadrature NI techniques for the evaluation of such ill-behaved, rapidly fluctuating functions, such as g ex (ξ) and h eρ (ξ) of Figs. 3, 4 and 6 [40], [41]. Finally, the effect of the frequency on the convergence times is apparent. Depending on the required error allowance, it seems that above a certain frequency level, the selection of an adaptive quadrature technique, such as the Adaptive Simpsons, used in our case, might be necessary for getting timely results.

III. EVALUATING A NOVEL ASYMPTOTIC SOLUTION TO
THE SOMMERFELD'S PROBLEM Now that we have a solid method for the numerical calculation of Sommerfeld Integrals, we may use it to examine a newly introduced asymptotic solution to the well-known Sommerfeld Radiation Problem. The method was first presented in [25] and briefly discussed below for ease of reference.

A. Synopsis of the Asymptotic Method
Using the rigorous mathematical analysis of [25], the field scattered by a planar interface, can be expressed as where, with respect to Fig. 1,ê θ2 =ê ρ cos θ 2 −ê x sin θ 2 refers to the unit vector, along the θ 2 -direction of a spherical coordinate system, whose origin is the dipole's image (A ) and R (θ 2 ) is given by (10), for ξ = θ 2 . Moreover, in (20) ζ p = ξ p −θ 2 , where ξ p is the pole of R (θ 2 ). Also, notice that (20) is derived under the usual case scenario, where σ ωε 0 , in which case ξ p may be approximated by The most interesting part in (20) is special function X, the so-called 'Etalon Integral' [26]- [31]. For parameters k, α, it is defined as the contour integral along path S of Fig. 7 11 . The 'Etalon Integral' has useful properties and as shown in (22), it can be expressed in terms of Fresnel Integrals, which enable its easy evaluation via the complementary error function. Keep in mind that to reach (20), the Saddle Point Method was used, in order to deform the original Sommerfeld contour of integration, S z , into S, so as the expression for the Etalon Integral, (22), could be used. Therefore, the method is still another another "high freq." asymptotic method. The procedure is described in detail in [25], of which Fig. 3 is replicated as Fig. 7, in this manuscript, such that the relevant contours and the mapping process are briefly clarified. Finally, notice that the pole ξ p of R (θ 2 ) does not influence the result, since it is kept outside the contour of integration. This is why the condition σ ωε 0 that ensures the above argument, is important in our case.
It is also possible to further elaborate on (22), if one applies the large and small argument approximations for the complementary error function [32]. As a result the following asymptotic formulas are obtained, which when applied to (20), i.e. for k = k 01 r 2 and α = −ζ p , they yield the following analytic expressions, , ϕ → 0 (26) 11 Regarding the notation in (22): Expression (25) indicates the geometric optics reflected field, emanating from A , the dipole's image point (Fig. 1). It should be accurate, for a long electric distance, k 01 ·r 2 , i.e. at the far field region, provided that at the same time the grazing angle ϕ = π/2−θ 2 is not very small. In [18], we also reached (25), using the SPM method. However, as stated in [24], the SPM required only the fulfillment of a large electric distance. The effect of the grazing angle was essentially overlooked and hence the propagation mechanism for the case of low height transmission link (where the angle of incidence is small) could not be highlighted. Pay attention to the fact that if (25) was absolutely accurate, even for sliding angles of incidence, just because of a high frequency transmitting source, the field to be received would essentially be imperceptible, since, in this case, the reflection coefficient, R , approaches to −1 and E R would simply cancel E LOS .
Regarding (26), we are given with an expression, describing the behavior of the scattered field for sliding angles of incidence. It appears to have surface wave characteristics, due to the existence of the exponentially decaying factor, with respect to the altitude, e −δk01(x+x0) . However, as mentioned in [25], for the derivation of (26) various assumptions and approximations were made, relating the electrical and geometrical characteristics for the problem (see also Section IV). The validity of these assumptions remain to be validated.
In the simulations that follow, we compare the closed-form asymptotic solution of [25], i.e. (20), against the SPM-based solution of [18], which essentially leads to the geometrical optics field expressions, in the high frequency regime. The reference for our comparisons will be the numerical integration results that we obtain for the EM field, using the methodology of Section II, above that is, the evaluation of (18), (19) and the respective formulas for the magnetic field. In addition, we also examine and comment on the accuracy of the analytic expressions (25), (26).

B. Simulation Results
We exhibit two sets of simulations in Figs. 8 and 9, below. Fig. 8 demonstrates the effect of the frequency on the total received electric field, for a number of scenarios, regarding the Transmitter -Receiver (T-R) horizontal distance, denoted with "d" in the respective plots. With the exception of Fig.  8(f), the basic simulation parameters are shown in Table II. The ground parameters, ε r , µ, σ, are indicative for the case of sea water and do fulfill the basic requirement, σ ωε 0 , mentioned in Section III-A, for all the examined cases. The altitudes X 0 and X are kept constant, at 60 m and 15 m respectively, however by increasing the horizontal distance, d (up to 30 km in Fig. 8(e)), we essentially simulate sliding angles of incidence as well. Only in the case of Fig. 8(f), where the frequencies involved are significantly lower, did we further lower the antennas' heights and this was done to examine the degree to which the methods are able to detect the so called surface wave, which in this case should be more significant [35]. We also focus on the far-field behavior and for this reason we don't exhibit sub-wavelength scenarios, although our simulations have revealed that some of our findings could Adaptive Simpsons relative error tolerance 10 −6 a much smaller than the wavelength λ = c/f b pertains for the case of see water be extended to near-field region as well 12 . Finally, for the complementary error function in (22), which due to (20) now includes a complex argument (−ζ p ), the algorithms described in [42], [43] were utilized, which very accurately evaluate such special functions in the complex plane. The case shown in Fig. 8(a) is indicative of non-near ground level terrestrial communication. The T-R relative position is such that the angle of incidence is ϕ 15 • . It is evident that there is an almost perfect match between the results obtained numerically, labeled as "NI" in the plots and what is predicted by the newly introduced asymptotic solution (20), depicted via the "Etalon" indicated lines in Figs. 8 and 9. It is also equally interesting that the older asymptotic SPM method, also yields similar results, which for frequencies around 20 MHz and above are consistently almost identical with what is numerically computed. Keep in mind that the SPM solution is essentially the expression given by (25), which is derived as a special case of (20), as already mentioned in Section III-A. However, the restriction for (25), namely √ k 01 r 2 · sin ϕ 2 1, is not strictly fulfilled in our case. For the scenario of Fig. 8(a) it goes from 0.31@1MHz to 3.1@100MHz. At 20MHz this quantity is about 1.4. Therefore, it seems that (25) is an accurate analytic expression, to be used for non-sliding angle of incidence reception, whose validity could be practically extended beyond the strict restrictions imposed for its derivation.
In Fig. 8(b) the T-R distance is increased to 3 km and as a result the angle of incidence is radically reduced to ϕ 1.43 • . In this scenario, we do observe a discrepancy between the two asymptotic solutions and of both of them with the reference numerical integration (NI) results for (19), which pertain to the complete solution for the Sommerfeld's Radiation problem. Of course, this discrepancy appears to be relatively small and if examined in a broader frequency range, as in Fig.8(c), it may be practically regarded negligible. Nevertheless, it is important to note the tendency of (20) to better follow (19), something that is even more apparent in the diagrams (d) and (e) of Fig.8. In these cases, the T-R distance is further increased to 10 km and 30 km, with the incidence angles now being as sliding as ϕ 0.43 • and ϕ 0.14 • respectively. Overall, compared with the asymptotic solution of [18] (SPM   Figs. 9(e), (f) the Scattered Field is illustrated. The term Etalon "Surf", refers to the evaluation of (26). Fig. 9(a), exhibits the magnetic field. The rest of the labeling convention of Fig. 8 applies. based solution), the recently introduced asymptotic solution in [25] (Etalon based), is a better estimate to the total solution of the Sommerfeld's radiation problem. It is also apparent that both methods smoothly converge to (19) in the high frequency regime, however the solution of [25] converges faster. On the contrary, Fig. 8(f) verifies a somehow expected behavior. In lower frequencies, both methods fail to describe the propagation mechanism, for being unable to capture the effect of the so-called surface wave, which in this scenario should be rather significant 13 . Indeed, in this case, (20) behaves only marginally better than the respective asymptotic formula of [18], which as already stated essentially yields the space wave component, ignoring the contribution of the surface wave. We mentioned before that the results of Fig. 8(f) were somehow expected, since if one follows the derivation process of [25], he/she will identify the use of the Saddle Point method for reaching the final formulas, which may therefore yield accurate results only in the high frequency regime [44].
To confirm and further solidify the above arguments, in Fig. 9 we exhibit the field behavior from the perspective of a varying T-R distance. Starting from the low frequencies scenario, in Fig. 9(a) 14 it is apparent that both asymptotic methods fail to produce accurate results. Actually, according to our detailed simulations, this situation holds true almost up to approximately 1 MHz. Moving, towards the HF frequency zone, Fig. 9(b), the advantages of the newly introduced asymptotic solution show up. The difference between (20) and the previous SPM-based solution of [18], is more evident for large distances, where the effect of the scattered field to the total field is more significant; hence, the improvement that the "Etalon" function, X, yields in (20) becomes more visible. Finally, if we further proceed to the VHF zone of Fig. 9(c), we realize that both asymptotic methods begin to converge and ultimately they coincide with the complete solution at even larger frequency bands, as indicatively shown in Fig. 9(d). At those frequencies and in accordance with what is known in the literature, the surface wave is almost negligible. Therefore, there is almost nothing extra left for special function X to expose and (20) simply yields the reflected field, exactly as the asymptotic solution of [18] does.
The last two diagrams of Fig. 9 are devoted to the investigation of (26), an interesting expression, which attributes surface wave characteristics to the near-ground-level scattered field, never encountered before, in this analytic form, in the literature. For such purpose, extended simulations were run, the outcome of which may be summarized as follows: up to the HF frequency band, (26) does converge to (20), from which it was derived when ϕ → 0. In the case of Fig. 9(e), the convergence occurs approximately at 7 km, which for the selected T-R altitudes, is equivalent to ϕ 0.61 • . However, for different scenarios, regarding the T-R heights, the required value for ϕ may change to up to 2 • approximately. Of course, as already pinpointed in Fig. 9(b), in this frequency band, (20) and therefore (26) as well, are not accurate approximations of the complete solution. They are simply better estimates compared to the SPM approximation. On the contrary, at higher frequencies, we experience a total failure of (26) to follow the behavior of (20). This situation is illustrated in Fig.  9(f). Regardless of the selected T-R heights and their horizontal distance 15 , (26) always yields significantly underestimated values. Overall, we were unable to find a set up for which (26) can provide meaningful results and definitely a reconsideration of it is required.

IV. CONCLUSION AND FUTURE RESEARCH
We demonstrated an efficient method for the numerical evaluation of Sommerfeld Integrals in the spectral domain. The method proves fast and accurate and when applied to the evaluation of the EM field of a radiating vertical dipole above flat lossy ground, it fits with existing asymptotic solutions and Norton's results.
With a reference numerical method to accurately evaluate the integral representation to the Sommerfeld's ratiation problem, i.e. (18), (19), we then focused on the evaluation of a recently developed asymptotic solution. The solution uses the saddle point method and utilizes the properties of the socalled 'Etalon Integral', as a means to increase the accuracy of the results. Through extensive simulations, we verified that for the usual case, where σ ωε 0 , the method does succeed to provide better estimates to the complete problem, as compared with a more basic asymptotic approach, which is based on the application of the stationary phase method and essentially yields the well-known geometric optics field or Space Wave. Moreover, further asymptotic properties for the 'Etalon Integral' allowed us to reach analytic formulas for the scattered field. Unfortunately, only the analytic expression that pertains to the non-sliding angle of incidence case is validated.
Based on the findings of this work, we intend to further investigate (26) and identify the reason for its mismatch vs (20). Possibly, this is related with an assumption made in [25], through which the infinitesimal quantity, δ = ωε 0 ε 1 /2σ of (26), was essentially related with angle θ 2 . This correlation may be arbitrary since a parameter associated with the electric characteristics of the problem is related with the geometry set up and this relation affects the phase of (20), as described in Appendix E of [25]. It is known that approximations made to the phase component of rapidly oscillating complex functions can be very sensitive, with respect to the accuracy of the final outcome. Probably this is why (26) does ultimately follow (20) in Fig. 9(e) but completely fails to do so in Fig. 9(f). In the latter case, the frequency is ten times bigger and hence (20) fluctuates too fast, for making the assumptions in (26) invalid.
As mentioned in [25], the ultimate goal is to provide useful asymptotics, applicable for every possible scenario, not just only for the usual σ ωε 0 case, considered here. For that purpose, we will insist on the investigation of special function X (k, α)and its properties, as well as other special functions that could be used to describe the behavior of Sommerfeld's integral expressions. Finally, a similar analysis for the case of a horizontal radiating dipole above flat lossy ground is also to be considered.

ACKNOWLEDGMENT
The Authors would like to thank Prof. George J. Fikioris, of the National Technical University of Athens, for making useful discussions and constructive comments towards the preparation of the present work.