K расчету распространения электромагнитных импульсов 
Electromagnetic Pulse Propagation over Nonuniform 
Earth Surface: Numerical Simulation 
Alexei V. Popov and Vladimir V. Kopeikin 
Pushkov Institute of Terrestrial Magnetism,  
Ionosphere and Radio Wave Propagation 
IZMIRAN, Troitsk, Moscow region, 142190 Russia 
popov@izmiran.ru, kopeikin@izmiran.ru  
 1. Introduction 
Parabolic equation method proposed by Leontowich and Fock [1,2] is an efficient 
simulation approach to VHF propagation over the earth surface. Deep physical analysis and 
advanced mathematical methods [3,4] turned Leontovich’s PE into a universal tool of 
diffraction theory. Its applications go far beyond the initial problem circle – e. g. [5-8]. The 
key role in this development played the decisive turn to straightforward numerical techniques 
pioneered by Malyuzhinets and Tappert [9,10]. 
 In radio wave propagation, PE was used first to derive explicit analytical formulae for 
the EM field strength in model environments. A simplification has been reached by 
introducing the impedance boundary condition (BC) [11]. Taking into account tropospheric 
refraction ducts required the use of sophisticated asymptotic methods [12]. Further 
development (almost exclusively towards numerical implementation) was aimed at refined PE 
modifications [13-15], account for irregular terrain [16], introducing artificial transparent 
boundaries [17,18] and nonlocal BC to describe rough interfaces [19]. A non-stationary PE 
counterpart and a finite-difference (FD) scheme for its solution have been proposed by 
Claerbout and applied to seismic problems [13]. Afterwards, this “time-domain parabolic 
equation” (TDPE) was used to calculate acoustic propagation in ocean [20]. At the same time, 
little attempts of using TDPE to simulate EM pulse propagation in realistic environments are 
known. 
 In this paper, we consider computational aspects of EM pulse propagation along the 
nonuniform earth surface. For ultrawide-band pulses without carrier, TDPE results directly 
from the exact wave equation written in a narrow vicinity of the wave front. To solve it by 
finite differences, we introduce a time-domain analog of the impedance BC and a nonlocal 
BC of transparency reducing the open computational domain to a strip of finite width. 
Numerical examples demonstrate the influence of soil conductivity on the received pulse 
waveform which can be used in remote sensing.  
mailto:popov@izmiran.ru
mailto:kopeikin@izmiran.ru
 For a high-frequency modulated EM pulse, TDPE arises as a convolution of PE 
solutions with the pulse envelope spectrum. In order to overcome computational difficulties, 
we develop an asymptotic approach based on the ray structure of the monochromatic wave 
field calculated at the carrier frequency. To accommodate complex-valued asymptotic 
solutions to the real initial condition we use the “analytic signal” approach introduced by 
Vainstein, Heyman and Felsen [21, 22]. An explicit solution of the time-domain transport 
equation reduces the computational procedure to numerical integration of standard PE at the 
carrier frequency and evaluation of a given 1D function in time domain. This diminishes 
computational expenses by 2-3 orders of magnitude and allows for pulsed wave field 
calculation in vast domains measured by tens of thousands wavelengths. As an example, we 
consider a problem of target altitude determination from overland radar data [23].  
This work has been done in collaboration with the Institute for high-frequency 
technique (IHF), Stuttgart University. Preliminary results appeared as short papers [24,25], a 
Russian version has been published in [26]. We dedicate this publication to the memory of 
Leopold Benno Felsen. 
 2. Monochromatic wave propagation 
 Omitting technical details and method refinements – see [12, 19], recall PE based 
scheme of monochromatic wave propagation over a smoothly nonuniform earth surface 
 - Fig.1.  )(xhz =
0 z                                             0 1 00 , = = σ ε 
                                      x 
ε, σ 
)(xhz =
Fig. 1. Elevated source illuminating smoothly rolling terrain (sketch). 
Horizontal magnetic component ),( zxHH y =  satisfies Helmholtz equation 
 0~22
ε           (1) 
with complex permittivity ωσπεε i4~ += , where σ  is soil conductivity in the Gaussian 
units set. In the upper medium 0,1 == σε  and the contact conditions at  are )(xhz =
, .         (2) 
where xzn ∂∂+∂∂=∂∂ αα sincos , )('arctan xh=α . At large distances from the source the 
wave field is sought as a plane wave with slowly varying complex amplitude: 
 [ )(exp),,(),,( tkxikzxutzxH ]ω−≈          (3) 
Here, λπω 2≡= ck  is the wave number, and the complex “attenuation function”  
satisfies the Leontovich PE 
),,( kzxu
 )(,02 2
ik >=
         (4) 
In 3D, divergence factor x1  must be added in (3). In this paper we use Gaussian initial 
condition 
( ) (
−−−⎟⎟
exp),,0( zzikzz
kzu β
)       (5) 
corresponding to an exact solution of PE (4) 
= xzz
kzxui 22
exp),,(
     (6) 
- skewly propagating Gaussian beam with initial width  and wave front radius 0w 0ρ , 
determined by complex parameter ( ) 12000 21
+= wkix ρ ;  β  being a small elevation angle.  
 Impedance approximation is based on wave beam contraction when entering a denser 
dielectric medium. Standard Leontovich BC [11] 
 0,~ =−=∂
                     (7) 
results from the contact conditions (2) under the assumption of almost vertical propagation in 
the lower medium: ( )ε~exp),( ikzTzxH −≈− . For grazing angles (Fig. 2) this assumption 
breaks and a plane incident wave ( )[ ]ββ coscosexp),( zxikzxH −=+  with small 1<<β   
Fig. 2. Derivation of impedance BC for grazing angles. 
enters the half-space  close to the total internal reflection angle 0<z ( )210 ~arccos εγ = : 
0 sin~2 γε
γγ +≈ . Hence ensue 
( )[ ]
( )[ ] (
( )[ ] ( )000
exp)(
expexp)(
sincos
~exp)(),(
tgzxVtgzxikT
ctgzxk
tgzxikT
zxikTtxH
      (8) 
In virtue of the superposition principle, Eq. (8) holds for an arbitrary paraxial wave packet 
with the corresponding slowly varying function . Eliminating the latter by differentiation 
and making use of (2), we obtain 
 0,0~~
1~ ==+
ε        (9) 
This modified impedance BC provides a more accurate approximation of the reflection 
coefficient, especially, in a vicinity of the Brewster angle ( ) 210 1~arcsin
+= εβ .  
Fig. 3 allows one to compare the exact Fresnel reflection coefficient  
cos~sin~
cos~sin~
=FR                 (10а) 
with the Leontovich approximation 
1sin~
1sin~
βLR                   (10b) 
and that resulting from the modified impedance BC (9)  
Fig.3. Comparison between exact Fresnel reflection coefficient 
                          and impedance   approximations. 
( ) ( )
( ) 1~cos~sin~
1~cos~sin~
εβεβε
εβεβε
βMR                (10c) 
 Taking into account the boundary tilts )(xh′  and using "parabolic" approximation (3), 
we derive a modified BC for the attenuation function : ),( zxu
 )(,0)('~
xhzuxhik
      (11) 
Contrary to the standard Leontovich BC (9), here it is not necessary to assume 1~ >>ε  - Eq. 
(11) breaks only for 11~ <<−ε  when nonlocal effects of wave interaction in both media are 
to be taken into account [19]. The impedance BC grants uniqueness of the solution of PE (4). 
In fact, calculating the energy flow through a vertical cross section one obtains 
),()(
dzzxuxI , 0),(~
−= hxu
   (12) 
which proves stability and uniqueness of the boundary value problem solution. 
 Finite-difference methods of PE solution have been studied in early works by 
Malyuzhinets and coauthors [9, 27]. Further method development is described in monographs 
[5, 19]. We employ a six-point implicit FD scheme supplemented with the impedance BC 
(11) at  and a discrete approximation of the nonlocal transparency BC [17,28] 
imposed at the artificial computational boundary 
)(zhz =
maxzz = : 
 ( ) ( )
maxmax ,
,       (13) 
An example of simulated VHF propagation over irregular terrain is illustrated by Fig. 4. 
Fig. 4. VHF attenuation function over irregular earth surface. 
 3. Radio pulse propagation: Fourier synthesis 
 A straightforward way to describe EM transients is to convolve monochromatic wave 
fields with the signal spectrum. In a 1D case, the propagating pulse is a superposition of plane 
waves  
 [ ]∫
−= dkekFtxH tkkxi )()(~
),( ω
       (14) 
In free space, kc=ω , and formula (14) yields a dispersion-less traveling wave 
−=−= dkekFsFxctFtxH iks)(~
)(),(),(
     (15) 
In a 2D environment, a natural generalization of the 1D solution (14) is a paraxial wave 
packet  
 dkekzxukFtzxH ctxik )(),,()(~
),,( −
      (16) 
where  is a solution of the PE (4) at a fixed frequency ),,( kzxu kc=ω .  
The superposition (16) will approximate an exact solution of the wave equation if the 
spectrum )(~ xF  is confined near a certain positive  satisfying the PE applicability 
conditions: 
LDk ππ 220 >>>>  where  and  are lateral and longitudinal characteristic 
scales of the problem. Consider a quasi-monochromatic pulse 
tctfctF 0cos)()( ω=  with 
duration 02 ωπ>>T . Its Fourier transform 
 ( ) ([ 00 ~~2
)(~ kkfkkfkF ++−= )],       (17) 
where  is the envelope  spectrum, contains negative frequencies 
not described by PE (4). Introducing complex signal 
= dtectfckf tiω)()(
)(ctf
)exp()()( 0tictfctFc ω−=  eliminates the 
second term in (17). Still, the remaining "positive" component ( )0
kkf −  centered at  may 
spread onto negative semi-axis. In order to avoid nonphysical effects of negative frequencies 
propagation, the "analytic signal" [21] can be used, defined as a one-side inverse Fourier 
transform of the truncated spectrum 
>++−=
kkkfkkfkF
kF      (18) 
Thus, by definition, the analytic signal is a Cauchy type integral 
FdkekFsF  ,     (19) 
regular in the lower half-plane 0Im <s . For real s, its real part coincides with  whereas 
the imaginary part is given by Hilbert transform  
 ( )∫
FsF V.P.
)(Im        (20) 
Introduction of the analytic signal violates the causality principle: the real signal  is 
zero before the moment of switching on the transmitter while  for . However, 
for a high-frequency radio pulse this discrepancy is small. So the analytic signal envelope 
defined as 
)(ctF
0)( ≠+ ctF 0<t
)(sF +  is close to  but, contrary to the "naïve" complex signal , admits 
analytic continuation into complex domain, compatible with asymptotic propagation laws 
[22].  As an example, consider a modulated high-frequency pulse  with  
)(sf )(sFc
sksfsF 0cos)()( =
the envelope 
      Fig. 5. Modulated radio pulse waveform (21).      Fig. 6. Analytic signal spectrum (18). 
               
0),exp(sin
sbsas
sf , 
ibkkba
222 −−+
=    (21) 
see Fig. 5. For  its length is ab ≈ acT π~=Λ . The envelope spectrum )(kf has a peak 
at  with 0=k Λ≈Δ π2k  and tends to zero for ∞→k  as ( )2kaO . Spectra )(~ kF  and 
)(~ kF +  are shown in Fig. 6; the analytic signal envelope )(~ sF +  is plotted in Fig. 7a,b for 
real and complex arguments.  
           (a)                (b)  
Fig. 7. Analytic signal envelope of real (a) and complex (b) arguments. 
For a wave packet 
−++ =
)(),,()(~
),,( dkekzxukFtzxH ctxik
,    (22) )(0 0)(),,0(
zikezAkzu Φ=
the “initial” condition 
 [ ] [ )()()(~
),,0( 00
)(0 0 zctFzAekF
tzH ctzik Φ−== +
]    (23) 
describes an analytic signal  with amplitude  and initial delay )(ctF + )(0 zA czzt )()( 00 Φ= . 
  Fig. 8. Example of a Fourier-synthesized propagating EM pulse. 
An example of modulated pulse propagation over smoothly rolling interface is 
depicted in Fig. 8. A sequence of snap-shots traces the evolution of the initial pulse envelope 
 defined by (5), (21) due to the incident Gaussian beam divergence and 
reflection from the curved interface
),()( 00 kzuctf
)(xhz = . It should be noted that Fourier synthesis is 
computationally efficient only for rather narrow-band pulses. In fact, for a good 
approximation of the convolution integral (22) one has to solve PE (4) for a set of wave 
numbers covering the spectral band kkkkk Δ+<<Δ− 00 , Λ≈Δ π2k  with a small frequency 
step ( kk Δ<<δ ) and, to avoid phantom solutions in the given range xΔ , even more restrictive 
condition must be posed: xk Δ<< πδ 2 . Adequate simulation methods for wide-band EM 
pulse propagation are discussed in the following sections. 
 4. Time-domain PE and boundary conditions 
Straightforward derivation shows that if ( )kzxu ,,  is a solution of PE (4), the transient 
wave packet (16) , as a function of variables ),,(),,( szxtzx Π≡Η xctszx −=,, , satisfies the 
Claerbout equation  
zsx ∂
        (24) 
Equation (24), usually called "time-domain parabolic equation" (TDPE), has been obtained in 
[13] by formal substitution sik ∂∂=  as well as by the reduction of the time-dependent wave 
equation 
      (25) 
in a narrow vicinity of the paraxial wave front ctx = . Introduction of scaled variables 
Λ−=== )(,, xctDzLx ηςξ , where L, D are computational domain length and width, 
 is spatial pulse length, yields Λ
ςηξ ∂
,  1<<Λ== DLDν     (26) 
Neglecting the small term  results in TDPE (24).This derivation clarifies the nature of 
the “time-domain parabolic equation”: 
)( 2νO
1) It is a hyperbolic equation written in a traveling coordinate frame ; ),,( szx
2) TDPE does not describe the backward moving waves; 
3) TDPE is a paraxial (narrow-angle) approximation valid in a narrow strip 1)( <<= νOLD ; 
4) TDPE describes short pulses 1)( <<=Λ νOD  whose length Λ  is comparable with the 
wave front deviation from the plane ctx =  (Fig. 9); 
5) TDPE solutions are not necessary modulated high-frequency signals – they can represent 
short ultrawide-band pulses  without carrier, e. g. a damped sinusoid (21).  )(ctf
Fig. 9. Derivation of TDPE (24). 
Here, a seeming contradiction may arise, as the spectral maximum of  can lie in 
the vicinity of zero frequency, not described by PE (4). As a matter of fact, at small distances 
from the wave front 
)(Λ= Os  the main part of the pulse energy is determined by the high-
frequency edge of its spectrum Dak π2~ >>  satisfying PE applicability conditions. 
 To solve TDPE (24), an FD scheme of the 2nd order of accuracy has been proposed in 
[13]: 
( ) ( )lllllllll
rrrrrrrr
,,11,1,1
2,,11,1, )(4
nnnnznnnn zsx
Π+Π+Π+Π∇
=Π+Π−Π−Π
ΔΔ ++++++++
 (27) 
Here, ,  { }mnn ll
,, Π=Π ,xnxn Δ= zmzm Δ= , ss Δ= ll ; ( ) 1,,1,,2 2 −+ Π+Π−Π=Π∇ mnmnmn
nz llll . This 
equation is solved by zigzag marching in  plane between boundary values ),( sx
( ) [ ] 0)0,,(,)()(,,0 00 =ΠΦ−=Π zxzctfzAsz      (28) 
(given source and causality condition). At each marching step , a three-diagonal linear 
equation set arises for the unknown vector 
),( mn
1,1 ++Π l
n . In order to complete the boundary value 
problem, we have to add a correct BC taking into account soil properties and to find a way of 
the domain truncation without creating spurious reflections. Both problems are resolved by 
applying Fourier transform to the frequency-domain BC (11), (13). Consider a paraxial wave 
packet 
dkekzxukfszx iks−
∫=Π ),,()(
,      (29) 
satisfying the causality condition 0),,( =Π szx  for 0<s . We rewrite the impedance BC (11) 
emphasizing the dependence of complex permittivity kciσπεε 4~ +=  on the wave number 
ck ω= : 
 0)('
,      (30) 
where )1(2,4 −== επσεπσ cqcr . Multiplying Eq. (30) by  and applying Fourier 
transform (29), we get 
dkkekzxukf
xhshx
iks),,()(
),,()('),,(
.  (31) 
Substituting here the inverse Fourier transform 
        (32) ∫
),,(),,()(
ηη ηdezxkzxukf ik
we obtain, by standard calculations, the following expression for the RHS of (31): 
dsNhxshx
)(),,(),,(
;   ⎥
−+= ∫
tqrrs qr
qtIeqesN
)( )()(    (33) 
Thus, we have derived a nonlocal 2D boundary condition 
dsNhxshx
xhshx
)(),,(),,(
),,()('),,( ηηη
 (34) 
being an exact time-domain counterpart of the impedance BC (11). Its nonlocality is a 
consequence of interaction between two waves propagating along the earth surface with 
different phase velocities. The integral term kernel )( η−sN  can be easily calculated for 
different ε  and σ . For 15.3>ε , function  monotonously tends to zero with increasing s 
– see Fig. 10.  
Fig. 10. Kernel of nonlocal impedance BC (34). 
It is interesting to note that in both limiting cases: perfectly conducting boundary )( ∞=σ  
and zero soil conductivity )0( =σ  Eq. (34) reduces to a local BC. In the former case  
tends to a delta function, for . The integral in (34) limits to 
=∞→−=
1)(,)0( dssNqrN
),,( shx
, the RHS vanishes, and we get a Neumann BC. In the opposite case )0( →σ  the 
spatial scale of  (length of the pulse “dispersion tail”) is growing but its absolute value is 
tending to zero, so only a local term 
 remains. 
In a similar way, the time-domain generalization of the transparency BC (13) is 
derived which grants the absence of spurious reflections a from the artificial computational 
boundary . Applying to Eq. (13) the Fourier transform (29), (32) and denoting maxzz = ipk = , 
we get  
( ) ( )
( ) ∫∫∫
maxmax
  (35) 
Finally, evaluating the inner integral we obtain an elegant 2D boundary condition 
 ( ) ( )
( )( )ηξ
ηξπ −−∂∂
∫ ∫ sx
,, max
max     (36) 
symmetric with respect to the variables x, z , which could be expected from the symmetry of 
the TDPE (24). 
 A simulated example of ultrawide-band EM pulse propagation over a nonuniform 
earth surface with soil parameters 10=ε , )01.0(109 17 мSs −⋅=σ  is depicted in Fig. 11а. 
Evolution of the spatial amplitude distribution for a pulsed signal generated by a Gaussian 
source: ( )[ ]20200 exp)( wzzzA −−=  with a skew curved wave front: 
)(2)()( 00
00 zzzzz −+−=Φ βρ  is shown in a grey color scale. The initial pulse waveform 
 is a damped sinusoid (21) with )(ctf ba =  and spatial length ma 30==Λ π .  
   (a)      (b) 
Fig. 11. Propagation of ultrawide-band pulse (21) over nonuniform earth surface. 
           Initial Gaussian beam parameters: ,80,3000 mmz == 0w 1.0,3000 −== βρ m . 
           Soil conductivity: mS01.0=σ  (a), 0=σ  (b). 
Fig. 12. Received pulse waveform depending on soil conductivity: 
                mS01.0=σ  (solid line), mS001.0=σ  (dashed line); initial pulse (dots). 
The snapshots clearly show the reflected pulse generation at the earth surface. The 
transparency BC (36) imposed at the height mz 500max =  assures unimpeded radiation exit 
from the computational domain. Finite soil conductivity causes signal dispersion appearing in 
a certain delay of the reflected pulse. It is obvious from the comparison with a similar plot 
calculated for a model non-conducting soil 0=σ  - Fig. 11b. A quantitative estimation of the 
effect can be made by means of Fig. 12 revealing a considerable dependence of the pulse 
waveform on the soil conductivity. This effect caused by the pulse penetration into the ground 
can be used for ecological monitoring (water pollution, earthquake precursors, etc.). 
5. Hybrid TDPE and short high-frequency pulse propagation 
 An important practical issue is overland propagation of short EM pulses with high- 
frequency carrier. Basically, having absolute stability, TDPE (24) is capable to describe wide-
band radio pulse propagation. However, the computational expense is drastically growing 
with increasing carrier frequency. If 0ω  considerably exceeds the spectral band of the signal it 
is useful to factor out the carrier: 
 [ ])(exp),,(),,( 00 txkiszxUszx ω−=Π      (37) 
and to consider the transient signal envelope  satisfying a hybrid equation
),,( szxU
ik ,        (38) 
combining the features of standard Leontovich PE (4) with Claerbout TDPE (24). 
Unfortunately, despite a relatively slow variations of  in space-time, straightforward ),,( szxU
numerical solution of this hybrid TDPE (HPE) entails considerable difficulties, as the large 
coefficient  by 0k xU ∂∂  demands a dense computational grid. On the other hand, the 
presence of a big parameter allows us to construct an asymptotic solution of HPE (38) 
radically reducing the computational burden. 
 In order to find a proper asymptotic Ansatz, consider monochromatic PE (4) at the 
carrier frequency ck00 =ω with the initial condition  . We admit 
complex values of the eikonal  to describe relatively narrow wave packets, like a 
Gaussian beam (6). For , we obtain an asymptotic solution 
00)(),0( zikezAzu Φ=
)(0 zΦ
         (39) ),(0 0),(),,(
zxikezxAkzxu Φ=
where  satisfies a “parabolic” eikonal equation ),( zxΦ
1 2 =Φ+Φ zx          (40) 
while slowly varying amplitude  is governed by the paraxial transport equation ),( zxA
=Φ+Φ+ AAA zzzzx         (41) 
Eqs. (40)-(41), being an approximate form of the well-known laws of geometric optics (GO) 
[4], can be easily solved by the method of characteristics. Consider a particular solution of the 
eikonal equation (40) corresponding to a bundle of rays spreading from a central point 
: 0,0 zzx ==
zx −−+≈
=Φ 20
),(       (42) 
At a characteristic line xzz γ+= 0  we have xxx
),( γγ =Φ . The envelope of the family 
(42) solves the boundary value problem with an arbitrary initial condition . 
Define 
)(),0( 0 zz Φ=Φ
 [ ] xzzxzzx )(
)()(, 0
0000 γγ +Φ=+Φ       (43) 
By differentiating Eq. (43) with respect to x  and  we obtain 0z
 xx zzx ')'1(,2
2 γγγγγ +Φ′=Φ+=Φ+Φ       (44) 
Function (43) will satisfy the eikonal equation (40) if the ray direction is matched with the 
local wave front tilt: )()( 000 zz Φ′=γ . Having constructed the eikonal ),( zxΦ  we reduce the 
transport equation (41) to an ODE  
( ) ( ) 0,
0 =+′+
++ xzxA
γ      (45) 
with an evident integral 
xzzxA
=+      (46) 
In a similar way, an asymptotic solution of the modified Claerbout equation (38) can 
be found. Substituting the Ansatz [ ]),(exp),,(),,( 0 zxikszxIszxU Φ=  into Eq. (38) we obtain 
0 =−+⎟
⎛ Φ+Φ+Φ−+⎟
⎛ +Φ+Φ− xszzzzzzsxxx IIIIIIikikk   (47) 
The leading term ( )20kO  disappears in virtue of the eikonal equation (40). Thus, to the 
accuracy ( )10−kO , a space-time transport equation arises for the slowly varying amplitude 
: ),,( szxI
=Φ+Φ+Φ− IIII zzzzsxx      (48) 
As the coefficients of Eq. (48) do not depend on s, it has a solution of the following form  
        (49) [ ),(),(),,( zxsgzxAszxI Ψ−= ]
Here,  is a solution of the stationary transport equation (41) while  is an arbitrary 
function of 
),( zxA )(sg
xcts −= , and  satisfies a linear PDE ),( zxΨ
 xzzx Φ−=ΨΦ+Ψ          (50) 
Solving Eq. (50) by characteristics one easily gets  
 [ ] [ ] )()(,)(, 00000 zxzzxxzzx θγγ ++Φ=+Ψ       (51) 
where )( 0zθ is an arbitrary function. So, the solution of the HPE (38) has asymptotic 
representation 
 [ ] [ ] [ ] ),,()(),(),(exp),(),(~),,( 000
kzxuzzxsgzxikzxsgzxAszxU
θ−Φ−=ΦΨ−
 (52) 
Here,  is a solution of the standard Leontovich PE (4),  and  are its 
amplitude and eikonal, respectively;  and 
),,( 0kzxu ),( zxA ),( zxΦ
)(sg )( 0zθ  are arbitrary functions, and  is to 
be found from the transcendental equation 
),(0 zxz
zxzz =Φ′+ )( 000 . Asymptotic solution (52) is a 
paraxial version of the space-time GO [29], the rays and wave fronts being defined 
numerically via parabolic equation. Functions  and ),( zxA ),( zxΦ  are generally complex-
valued, so distinction between wave amplitude and complex “phase” is made solely on the 
basis of their different dependence on frequency. In particular, complex eikonal ,  
defined as  
),( zxΦ
),,(log),(
0 kkukki
izx     (53) 
is calculated from PE numerical solutions at two close frequencies ck 2,12,1 =ω . 
 Physically, complex eikonal in Eq. (52) appears due to diffraction effects described by 
PE (4). An important consequence is the absence of singularities in the constructed 
asymptotic solution, as the "parabolic" rays do not produce caustics in the real space. Another 
effect caused by diffraction – pulse envelope distortion also is taken into account via complex 
values of the signal delay czx ),(ψ . Physical meaning of complex xcts −=  is provided by 
the theory of analytic signal [22]. 
 Arbitrary functions in Eq. (52) are determined by the initial conditions. In the simplest 
case the constructed transient (37) has the form 
 [ ]),(),(),,(),,( zxxctFzxAxctzxtzxH Φ−−≈−Π≡ +     (54) 
where  and  are complex amplitude and eikonal evolved from the initial 
,  given by Eq. (23) and  is the analytic signal (19) 
corresponding to the real signal . Physical meaning has the real part of the complex 
solution (54) or, from the practical point of view, its normalized envelope 
),( zxA ),( zxΦ
)(0 zA )(0 zΦ )(exp)()( 0 tictftcF ω−≈
)(ctf
2),,( szxH . 
 In virtue of the superposition principle, a more general asymptotic solution can be 
constructed as a number of terms (54). That is a direct analogy with ordinary GO where the 
incident and reflected waves correspond to different ray families. An important practical 
example is radar pulse propagation over the earth surface when the direct and reflected from 
the ground pulses can be distinguished and used for target location [23].  
 Consider first a model example: a short pulse with carrier frequency  
and damped sinusoidal envelope (21), propagating over a slowly rolling boundary . 
Initial pulse parameters are: 
МHz2000 =f
)(xhz =
maba 9, ≈=Λ= π ; ,150 m=w  ,01.0−=β  m2000 =ρ . 
Stationary field distribution calculated by numerical integration of PE  (4) at cfk 00 2π=  
produces a regular interference pattern – Fig. 13. It can be represented as a superposition of  
  Fig. 13. Monochromatic attenuation function.    Fig. 14. Reflected wave phase distribution. 
the incident Gaussian beam (6)  with the reflected wave ),,( 0kzxui ir uukzxu −≡),,( 0  
determined by the terrain  and the impedance BC (11). Functions  
are given by the asymptotic solution (39)-(41), and eikonal 
)(xhz = ),(),,( zxzxA ii Φ
),( zxrΦ  is reconstructed from the 
spatial phase distribution of the reflected wave – see Fig. 14. In accordance with such a 
monochromatic framework, a two-term asymptotic formula arises for the pulsed transient:  
[ ] [ ]
[ ] [ ),(),(),(),(
),(),(),(),(~),,(
),(),( zxxctFzxeuzxxctFezxu
zxxctFzxAzxxctFzxAtzxH
ri Φ−−+Φ−−=
Φ−−+Φ−−
+Φ−+Φ−
 (55) 
Note that to find the amplitude and complex delay of the incident and reflected signals we 
need just to solve the standard PE in frequency domain at two close frequencies 02,1 ωω ≈  - 
see Eq. (53). In time domain, calculation reduces to the evaluation of an analytic function 
 for the given argument values of interest. That radically diminishes the required 
computational resources compared with direct numerical integration of the TDPE (24). 
)(sF +
Fig. 15. Received pulse envelope as function of receiver altitude. 
         (a)       (b) 
Fig. 16. Comparison between numerical solution of TDPE (24) (a) and asymptotic  
             solution (55)  (b). 
The envelope of the received analytic signal (55), as a function of xcts −=  and the 
receiver height z, for  is shown in Fig. 15. One can see profound interference 
minima near the earth surface and a good separation of the direct and reflected pulse to 
heights above 400 m. Figures 16a,b compare the asymptotic solution with direct numerical 
integration of the Claerbout TDPE (24). Qualitatively, they are almost identical. Some hardly 
seen discrepancy is due to a limited accuracy of the asymptotic solution and FD scheme (27). 
This comparison demonstrates the efficiency of the developed approach. Substantial 
acceleration of the numerical procedure (by around 200 times in this example) makes a good 
reason to use it in realistic conditions. 
kmx 7=
As an example, we simulate an experimental situation [23]: radar pulse propagation 
between two aircrafts flying by parallel routs over an irregular terrain. The experiment [23] 
was aimed at simultaneous determination of the target range and altitude from the measured 
return times of the direct radar pulse and the echo signal from the earth surface. Our goal is to 
develop an efficient method of EM field calculations under conditions of multipath and signal 
distortion. At such large ranges )100( kmX =  the Earth sphericity must be taken into 
account. For this purpose, a parabolic hump *2)( earthRxXx −  has been added to the real terrain 
profile plotted in .9 of [23]. Atmospheric refraction has been considered by using the 
equivalent Earth radius earthearth RR 34
* =  [12]. Global field strength distribution produced by 
the incident carrier wave at  is depicted in Fig. 17a. Despite evident multipath 
character of the reflected wave (Fig. 17b) its eikonal 
МHz1410 =f
),( zxrΦ  has a rather regular structure. 
Therefore, our PE based version of complex GO can be applied to simulate the averaged  
       (a)       (b) 
Fig. 17. Simulation of experiment [23]: global field strength (attenuation function) at 
, reflected wave (b). )(1410 aMHzf =
parameters of the received radar pulse (the actually observed signal is a stochastic quantity 
with normal distribution [23]). Its envelope, as a function of the distance from the paraxial 
wave front xcts −=  and the receiver height z, is depicted in Fig. 18.  
Fig. 18. Gaussian pulse envelope as function Fig. 19. Target height from reflected  
            of relative delay and receiver altitude.   pulse delay [23]. 
A Gaussian pulse waveform is chosen with the parameters corresponding to the experimental 
situation [23]:  ,141,3.5 00 МHz== fkmz m75≈Λ . The direct and reflected pulses are 
distinctly separated for  which allows one to reliably solve the triangle  
for target altitude determination - see Fig. 19 borrowed from [23]. Variability and statistics of 
kmz 5.4> 21,, RRR
the simulated reflected pulse resemble the experimental plots presented in [23], and the 
calculated received pulsed signal envelope for a fixed receiver height kmz 3.5=  (Fig. 20)  
  Fig. 20. Received pulse envelope: asymptotic HPE solution (solid line), 
                                     experiment [23] (dash), stochastic model calculation [23] (dots). 
agrees well both with the experimental data and the results of thorough statistical modeling 
[23]. The quantitative discrepancy in the reflected pulse amplitude does not exceed the 
inherent uncertainty due to the errors in terrain description.  
Acknowledgements 
This work was supported in part by a joint RFBR-DFG grant No 01-02-04003. The 
authors are grateful to Friedrich Landstorfer and Ningyan Zhu who initiated this research. 
References 
1. M.A. Leontovich, A new method to solve problems of EM wave propagation over the 
earth surface (in Russian). USSR Academy of Sciences Trans., Physics Series, vol. 8, 
No. 1, pp. 16-22, 1944. 
2. M.A. Leontovich, V.A. Fock, Solution of the problem of electromagnetic wave 
propagation along the Earth’s surface by the method of parabolic equation, J. Phus. 
USSR, vol. 10, pp. 13-23, 1946. 
3. G.D. Malyuzhinets, Progress in understanding diffraction phenomena, Soviet. Phys. 
Uspekhi., vol. 69, pp.321-334, 1959. 
4. V.M. Babič, V.S Buldyrev. Short-Wavelength Diffraction Theory: Asymptotic 
Methods. Springer-Verlag, Berlin, 1990 
5. D. Lee, A.D. Pierce, E.-C. Shang, Parabolic equation development in the twentieth 
century, J. Comput. Acoustics, vol. 8, No. 4, 527-637, 2000. 
6. L. A. Vainstein, Open resonators and open waveguides (in Russian). Soviet Radio, 
Moscow, 1966. 
7. M.D. Feit, J.A. Fleck, Jr. Light propagation in graded-index fibers. Appl. Optics, vol. 
17, pp. 3990-3998, 1978. 
8. Yu.V. Kopylov, A.V. Popov, A.V. Vinogradov. Application of the parabolic wave 
equation to X-ray diffraction optics, Optics Communications, v. 118, pp. 619-636, 
1995. 
9.  G.D. Malyuzhinets, A.V. Popov, Yu.N. Cherkashin. On the development of a 
computational method in diffraction theory. 3rd All-Union Symposium on Diffraction 
of Waves, Academy of Sciences, Moscow, 1964. 
10. F.D. Tappert. The parabolic approximation method. Lecture Notes in Physics, v.70: 
Wave propagation and Underwater Acoustics, pp.224-287, Springer-Verlag, New 
York, 1977. 
11. M.A. Leontovich. On approximate boundary conditions on the surface of well-
conducting bodies (in Russian). Investigations on Radiowave Propagation, part 2, pp. 
5-12. Academy of Sciences, Moscow, 1948 
12. V.A. Fock, Electromagnetic Diffraction and Propagation Problems, Pergamon Press, 
1965. 
13. J. F. Claerbout. Fundamentals of Geophysical Data Processing with Applications to 
Petroleum Prospecting, McGraw-Hill, New York, 1976. 
14. A.V.Popov, S.A. Hosiosky. On a generalized parabolic equation of diffraction theory 
(in Russian). J. Comp. Math. and Math. Phys., vol.17, No. 2, pp. 527-533, 1977. 
15. E.A. Polyansky, On the relation between solutions of Helmholtz and SchrÖdinger type 
equations (in Russian), J. Comp. Math. Math. Phys., vol.2, No.1, pp. 241-249, 1972. 
16. M.F. Levy. Parabolic equation modelling of propagation over irregular terrain. 
Electronics Letters, vol. 26, pp. 1153-1155, 1990. 
17. V.A. Baskakov, A.V. Popov. Implementation of transparent boundaries for numerical 
solution of the Schroedinger equation, Wave Motion, vol.14, No.1, pp.123-128, 1991. 
18. S.W. Marcus. A generalized impedance method for application of the parabolic 
approximation to underwater acoustics, J. Acoust. Soc. Am., vol. 89, pp.391-398, 
1991. 
19. M.F. Levy. Parabolic Equation Method for Electromagnetic Wave Propagation. IEE 
Electromagnetic Wave Series vol. 45, 2000. 
20. M.D. Collins. The time-domain solution of the wide-angle parabolic equation 
including the effect of sediment dispersion. J. Acoust. Soc. Am., vol. 84, No. 6, pp. 
2114-2125, 1988. 
21. L.A. Vainstein, D.E. Vakman. Frequency Discrimination in Oscillation and Wave 
Theory (in Russian). Nauka, Moscow, 1983. 
22. E. Heyman, L.B. Felsen. Gaussian beam and pulsed-beam dynamics: complex-source 
and complex-spectrum formulations within and beyond paraxial asymptotics. J. Opt. 
Soc. Am. A, vol. 18, No. 7, pp. 1588-1611, 2001. 
23. L.M. Zurk. Experimental observation and statistics of multipath from terrain with 
application to overland height finding. IEEE Trans. Antennas Propag., vol. 47, No. 1, 
pp. 121-131, 1999. 
24. A.V. Popov, V.V. Kopeikin, Ning Yan Zhu, F.M. Landstorfer. Modelling EM 
transient propagation over irregular dispersive boundary. Electronics Letters, vol. 38, 
No. 14, pp. 691-692, 2002. 
25. A.V. Popov, V.V. Kopeikin, F.M. Landstorfer. Full-wave simulation of overland radar 
pulse propagation. Electronics Letters, vol. 39, No. 6, pp. 550-552, 2003. 
26. A.V. Popov, V.V. Kopeikin. Prediction of electromagnetic pulse propagation along 
the earth surface (in Russian). Progress of Modern Radioelectronics, No. 1’2005, pp. 
20-35. Radiotechnika, Moscow, 2005. 
27. A.V. Popov. Solution of parabolic equation of diffraction theory by finite difference 
method. J. Comp. Math. and Math. Phys., vol.8, No. 5, pp. 1140-1144, 1968. 
28. A.V. Popov. Accurate modeling of transparent boundaries in quasi-optics. Radio 
Science, vol. 31, No. 6, pp. 1781-1790, 1996. 
29. V.M. Babič, V.S. Buldyrev, I.A. Molotkov. Space-Time Ray Method. Linear and 
Nonlinear Waves (in Russian). SPB University Press, St. Petersburg, 1995. 

  We simulate EM pulse propagation along the nonuniform earth surface using so
called time-domain parabolic equation. To solve it by finite differences, we
introduce a time-domain analog of the impedance boundary condition and a
nonlocal BC of transparency reducing open computational domain to a strip of
finite width. Numerical examples demonstrate influence of soil conductivity on
the wide-band pulse waveform. For a high-frequency modulated EM pulse, we
develop an asymptotic approach based on the ray structure of the monochromatic
wave field at carrier frequency. This radically diminishes the computation
costs and allows for pulsed wave field calculation in vast domains measured by
tens of thousands wavelengths.

Introduction 
Parabolic equation method proposed by Leontowich and Fock [1,2] is an efficient 
simulation approach to VHF propagation over the earth surface. Deep physical analysis and 
advanced mathematical methods [3,4] turned Leontovich’s PE into a universal tool of 
diffraction theory. Its applications go far beyond the initial problem circle – e. g. [5-8]. The 
key role in this development played the decisive turn to straightforward numerical techniques 
pioneered by Malyuzhinets and Tappert [9,10]. 
 In radio wave propagation, PE was used first to derive explicit analytical formulae for 
the EM field strength in model environments. A simplification has been reached by 
introducing the impedance boundary condition (BC) [11]. Taking into account tropospheric 
refraction ducts required the use of sophisticated asymptotic methods [12]. Further 
development (almost exclusively towards numerical implementation) was aimed at refined PE 
modifications [13-15], account for irregular terrain [16], introducing artificial transparent 
boundaries [17,18] and nonlocal BC to describe rough interfaces [19]. A non-stationary PE 
counterpart and a finite-difference (FD) scheme for its solution have been proposed by 
Claerbout and applied to seismic problems [13]. Afterwards, this “time-domain parabolic 
equation” (TDPE) was used to calculate acoustic propagation in ocean [20]. At the same time, 
little attempts of using TDPE to simulate EM pulse propagation in realistic environments are 
known. 
 In this paper, we consider computational aspects of EM pulse propagation along the 
nonuniform earth surface. For ultrawide-band pulses without carrier, TDPE results directly 
from the exact wave equation written in a narrow vicinity of the wave front. To solve it by 
finite differences, we introduce a time-domain analog of the impedance BC and a nonlocal 
BC of transparency reducing the open computational domain to a strip of finite width. 
Numerical examples demonstrate the influence of soil conductivity on the received pulse 
waveform which can be used in remote sensing.  
mailto:popov@izmiran.ru
mailto:kopeikin@izmiran.ru
 For a high-frequency modulated EM pulse, TDPE arises as a convolution of PE 
solutions with the pulse envelope spectrum. In order to overcome computational difficulties, 
we develop an asymptotic approach based on the ray structure of the monochromatic wave 
field calculated at the carrier frequency. To accommodate complex-valued asymptotic 
solutions to the real initial condition we use the “analytic signal” approach introduced by 
Vainstein, Heyman and Felsen [21, 22]. An explicit solution of the time-domain transport 
equation reduces the computational procedure to numerical integration of standard PE at the 
carrier frequency and evaluation of a given 1D function in time domain. This diminishes 
computational expenses by 2-3 orders of magnitude and allows for pulsed wave field 
calculation in vast domains measured by tens of thousands wavelengths. As an example, we 
consider a problem of target altitude determination from overland radar data [23].  
This work has been done in collaboration with the Institute for high-frequency 
technique (IHF), Stuttgart University. Preliminary results appeared as short papers [24,25], a 
Russian version has been published in [26]. We dedicate this publication to the memory of 
Leopold Benno Felsen. 
 2. Monochromatic wave propagation 
 Omitting technical details and method refinements – see [12, 19], recall PE based 
scheme of monochromatic wave propagation over a smoothly nonuniform earth surface 
 - Fig.1.  )(xhz =
0 z                                             0 1 00 , = = σ ε 
                                      x 
ε, σ 
)(xhz =
Fig. 1. Elevated source illuminating smoothly rolling terrain (sketch). 
Horizontal magnetic component ),( zxHH y =  satisfies Helmholtz equation 
 0~22
ε           (1) 
with complex permittivity ωσπεε i4~ += , where σ  is soil conductivity in the Gaussian 
units set. In the upper medium 0,1 == σε  and the contact conditions at  are )(xhz =
, .         (2) 
where xzn ∂∂+∂∂=∂∂ αα sincos , )('arctan xh=α . At large distances from the source the 
wave field is sought as a plane wave with slowly varying complex amplitude: 
 [ )(exp),,(),,( tkxikzxutzxH ]ω−≈          (3) 
Here, λπω 2≡= ck  is the wave number, and the complex “attenuation function”  
satisfies the Leontovich PE 
),,( kzxu
 )(,02 2
ik >=
         (4) 
In 3D, divergence factor x1  must be added in (3). In this paper we use Gaussian initial 
condition 
( ) (
−−−⎟⎟
exp),,0( zzikzz
kzu β
)       (5) 
corresponding to an exact solution of PE (4) 
= xzz
kzxui 22
exp),,(
     (6) 
- skewly propagating Gaussian beam with initial width  and wave front radius 0w 0ρ , 
determined by complex parameter ( ) 12000 21
+= wkix ρ ;  β  being a small elevation angle.  
 Impedance approximation is based on wave beam contraction when entering a denser 
dielectric medium. Standard Leontovich BC [11] 
 0,~ =−=∂
                     (7) 
results from the contact conditions (2) under the assumption of almost vertical propagation in 
the lower medium: ( )ε~exp),( ikzTzxH −≈− . For grazing angles (Fig. 2) this assumption 
breaks and a plane incident wave ( )[ ]ββ coscosexp),( zxikzxH −=+  with small 1<<β   
Fig. 2. Derivation of impedance BC for grazing angles. 
enters the half-space  close to the total internal reflection angle 0<z ( )210 ~arccos εγ = : 
0 sin~2 γε
γγ +≈ . Hence ensue 
( )[ ]
( )[ ] (
( )[ ] ( )000
exp)(
expexp)(
sincos
~exp)(),(
tgzxVtgzxikT
ctgzxk
tgzxikT
zxikTtxH
      (8) 
In virtue of the superposition principle, Eq. (8) holds for an arbitrary paraxial wave packet 
with the corresponding slowly varying function . Eliminating the latter by differentiation 
and making use of (2), we obtain 
 0,0~~
1~ ==+
ε        (9) 
This modified impedance BC provides a more accurate approximation of the reflection 
coefficient, especially, in a vicinity of the Brewster angle ( ) 210 1~arcsin
+= εβ .  
Fig. 3 allows one to compare the exact Fresnel reflection coefficient  
cos~sin~
cos~sin~
=FR                 (10а) 
with the Leontovich approximation 
1sin~
1sin~
βLR                   (10b) 
and that resulting from the modified impedance BC (9)  
Fig.3. Comparison between exact Fresnel reflection coefficient 
                          and impedance   approximations. 
( ) ( )
( ) 1~cos~sin~
1~cos~sin~
εβεβε
εβεβε
βMR                (10c) 
 Taking into account the boundary tilts )(xh′  and using "parabolic" approximation (3), 
we derive a modified BC for the attenuation function : ),( zxu
 )(,0)('~
xhzuxhik
      (11) 
Contrary to the standard Leontovich BC (9), here it is not necessary to assume 1~ >>ε  - Eq. 
(11) breaks only for 11~ <<−ε  when nonlocal effects of wave interaction in both media are 
to be taken into account [19]. The impedance BC grants uniqueness of the solution of PE (4). 
In fact, calculating the energy flow through a vertical cross section one obtains 
),()(
dzzxuxI , 0),(~
−= hxu
   (12) 
which proves stability and uniqueness of the boundary value problem solution. 
 Finite-difference methods of PE solution have been studied in early works by 
Malyuzhinets and coauthors [9, 27]. Further method development is described in monographs 
[5, 19]. We employ a six-point implicit FD scheme supplemented with the impedance BC 
(11) at  and a discrete approximation of the nonlocal transparency BC [17,28] 
imposed at the artificial computational boundary 
)(zhz =
maxzz = : 
 ( ) ( )
maxmax ,
,       (13) 
An example of simulated VHF propagation over irregular terrain is illustrated by Fig. 4. 
Fig. 4. VHF attenuation function over irregular earth surface. 
 3. Radio pulse propagation: Fourier synthesis 
 A straightforward way to describe EM transients is to convolve monochromatic wave 
fields with the signal spectrum. In a 1D case, the propagating pulse is a superposition of plane 
waves  
 [ ]∫
−= dkekFtxH tkkxi )()(~
),( ω
       (14) 
In free space, kc=ω , and formula (14) yields a dispersion-less traveling wave 
−=−= dkekFsFxctFtxH iks)(~
)(),(),(
     (15) 
In a 2D environment, a natural generalization of the 1D solution (14) is a paraxial wave 
packet  
 dkekzxukFtzxH ctxik )(),,()(~
),,( −
      (16) 
where  is a solution of the PE (4) at a fixed frequency ),,( kzxu kc=ω .  
The superposition (16) will approximate an exact solution of the wave equation if the 
spectrum )(~ xF  is confined near a certain positive  satisfying the PE applicability 
conditions: 
LDk ππ 220 >>>>  where  and  are lateral and longitudinal characteristic 
scales of the problem. Consider a quasi-monochromatic pulse 
tctfctF 0cos)()( ω=  with 
duration 02 ωπ>>T . Its Fourier transform 
 ( ) ([ 00 ~~2
)(~ kkfkkfkF ++−= )],       (17) 
where  is the envelope  spectrum, contains negative frequencies 
not described by PE (4). Introducing complex signal 
= dtectfckf tiω)()(
)(ctf
)exp()()( 0tictfctFc ω−=  eliminates the 
second term in (17). Still, the remaining "positive" component ( )0
kkf −  centered at  may 
spread onto negative semi-axis. In order to avoid nonphysical effects of negative frequencies 
propagation, the "analytic signal" [21] can be used, defined as a one-side inverse Fourier 
transform of the truncated spectrum 
>++−=
kkkfkkfkF
kF      (18) 
Thus, by definition, the analytic signal is a Cauchy type integral 
FdkekFsF  ,     (19) 
regular in the lower half-plane 0Im <s . For real s, its real part coincides with  whereas 
the imaginary part is given by Hilbert transform  
 ( )∫
FsF V.P.
)(Im        (20) 
Introduction of the analytic signal violates the causality principle: the real signal  is 
zero before the moment of switching on the transmitter while  for . However, 
for a high-frequency radio pulse this discrepancy is small. So the analytic signal envelope 
defined as 
)(ctF
0)( ≠+ ctF 0<t
)(sF +  is close to  but, contrary to the "naïve" complex signal , admits 
analytic continuation into complex domain, compatible with asymptotic propagation laws 
[22].  As an example, consider a modulated high-frequency pulse  with  
)(sf )(sFc
sksfsF 0cos)()( =
the envelope 
      Fig. 5. Modulated radio pulse waveform (21).      Fig. 6. Analytic signal spectrum (18). 
               
0),exp(sin
sbsas
sf , 
ibkkba
222 −−+
=    (21) 
see Fig. 5. For  its length is ab ≈ acT π~=Λ . The envelope spectrum )(kf has a peak 
at  with 0=k Λ≈Δ π2k  and tends to zero for ∞→k  as ( )2kaO . Spectra )(~ kF  and 
)(~ kF +  are shown in Fig. 6; the analytic signal envelope )(~ sF +  is plotted in Fig. 7a,b for 
real and complex arguments.  
           (a)                (b)  
Fig. 7. Analytic signal envelope of real (a) and complex (b) arguments. 
For a wave packet 
−++ =
)(),,()(~
),,( dkekzxukFtzxH ctxik
,    (22) )(0 0)(),,0(
zikezAkzu Φ=
the “initial” condition 
 [ ] [ )()()(~
),,0( 00
)(0 0 zctFzAekF
tzH ctzik Φ−== +
]    (23) 
describes an analytic signal  with amplitude  and initial delay )(ctF + )(0 zA czzt )()( 00 Φ= . 
  Fig. 8. Example of a Fourier-synthesized propagating EM pulse. 
An example of modulated pulse propagation over smoothly rolling interface is 
depicted in Fig. 8. A sequence of snap-shots traces the evolution of the initial pulse envelope 
 defined by (5), (21) due to the incident Gaussian beam divergence and 
reflection from the curved interface
),()( 00 kzuctf
)(xhz = . It should be noted that Fourier synthesis is 
computationally efficient only for rather narrow-band pulses. In fact, for a good 
approximation of the convolution integral (22) one has to solve PE (4) for a set of wave 
numbers covering the spectral band kkkkk Δ+<<Δ− 00 , Λ≈Δ π2k  with a small frequency 
step ( kk Δ<<δ ) and, to avoid phantom solutions in the given range xΔ , even more restrictive 
condition must be posed: xk Δ<< πδ 2 . Adequate simulation methods for wide-band EM 
pulse propagation are discussed in the following sections. 
 4. Time-domain PE and boundary conditions 
Straightforward derivation shows that if ( )kzxu ,,  is a solution of PE (4), the transient 
wave packet (16) , as a function of variables ),,(),,( szxtzx Π≡Η xctszx −=,, , satisfies the 
Claerbout equation  
zsx ∂
        (24) 
Equation (24), usually called "time-domain parabolic equation" (TDPE), has been obtained in 
[13] by formal substitution sik ∂∂=  as well as by the reduction of the time-dependent wave 
equation 
      (25) 
in a narrow vicinity of the paraxial wave front ctx = . Introduction of scaled variables 
Λ−=== )(,, xctDzLx ηςξ , where L, D are computational domain length and width, 
 is spatial pulse length, yields Λ
ςηξ ∂
,  1<<Λ== DLDν     (26) 
Neglecting the small term  results in TDPE (24).This derivation clarifies the nature of 
the “time-domain parabolic equation”: 
)( 2νO
1) It is a hyperbolic equation written in a traveling coordinate frame ; ),,( szx
2) TDPE does not describe the backward moving waves; 
3) TDPE is a paraxial (narrow-angle) approximation valid in a narrow strip 1)( <<= νOLD ; 
4) TDPE describes short pulses 1)( <<=Λ νOD  whose length Λ  is comparable with the 
wave front deviation from the plane ctx =  (Fig. 9); 
5) TDPE solutions are not necessary modulated high-frequency signals – they can represent 
short ultrawide-band pulses  without carrier, e. g. a damped sinusoid (21).  )(ctf
Fig. 9. Derivation of TDPE (24). 
Here, a seeming contradiction may arise, as the spectral maximum of  can lie in 
the vicinity of zero frequency, not described by PE (4). As a matter of fact, at small distances 
from the wave front 
)(Λ= Os  the main part of the pulse energy is determined by the high-
frequency edge of its spectrum Dak π2~ >>  satisfying PE applicability conditions. 
 To solve TDPE (24), an FD scheme of the 2nd order of accuracy has been proposed in 
[13]: 
( ) ( )lllllllll
rrrrrrrr
,,11,1,1
2,,11,1, )(4
nnnnznnnn zsx
Π+Π+Π+Π∇
=Π+Π−Π−Π
ΔΔ ++++++++
 (27) 
Here, ,  { }mnn ll
,, Π=Π ,xnxn Δ= zmzm Δ= , ss Δ= ll ; ( ) 1,,1,,2 2 −+ Π+Π−Π=Π∇ mnmnmn
nz llll . This 
equation is solved by zigzag marching in  plane between boundary values ),( sx
( ) [ ] 0)0,,(,)()(,,0 00 =ΠΦ−=Π zxzctfzAsz      (28) 
(given source and causality condition). At each marching step , a three-diagonal linear 
equation set arises for the unknown vector 
),( mn
1,1 ++Π l
n . In order to complete the boundary value 
problem, we have to add a correct BC taking into account soil properties and to find a way of 
the domain truncation without creating spurious reflections. Both problems are resolved by 
applying Fourier transform to the frequency-domain BC (11), (13). Consider a paraxial wave 
packet 
dkekzxukfszx iks−
∫=Π ),,()(
,      (29) 
satisfying the causality condition 0),,( =Π szx  for 0<s . We rewrite the impedance BC (11) 
emphasizing the dependence of complex permittivity kciσπεε 4~ +=  on the wave number 
ck ω= : 
 0)('
,      (30) 
where )1(2,4 −== επσεπσ cqcr . Multiplying Eq. (30) by  and applying Fourier 
transform (29), we get 
dkkekzxukf
xhshx
iks),,()(
),,()('),,(
.  (31) 
Substituting here the inverse Fourier transform 
        (32) ∫
),,(),,()(
ηη ηdezxkzxukf ik
we obtain, by standard calculations, the following expression for the RHS of (31): 
dsNhxshx
)(),,(),,(
;   ⎥
−+= ∫
tqrrs qr
qtIeqesN
)( )()(    (33) 
Thus, we have derived a nonlocal 2D boundary condition 
dsNhxshx
xhshx
)(),,(),,(
),,()('),,( ηηη
 (34) 
being an exact time-domain counterpart of the impedance BC (11). Its nonlocality is a 
consequence of interaction between two waves propagating along the earth surface with 
different phase velocities. The integral term kernel )( η−sN  can be easily calculated for 
different ε  and σ . For 15.3>ε , function  monotonously tends to zero with increasing s 
– see Fig. 10.  
Fig. 10. Kernel of nonlocal impedance BC (34). 
It is interesting to note that in both limiting cases: perfectly conducting boundary )( ∞=σ  
and zero soil conductivity )0( =σ  Eq. (34) reduces to a local BC. In the former case  
tends to a delta function, for . The integral in (34) limits to 
=∞→−=
1)(,)0( dssNqrN
),,( shx
, the RHS vanishes, and we get a Neumann BC. In the opposite case )0( →σ  the 
spatial scale of  (length of the pulse “dispersion tail”) is growing but its absolute value is 
tending to zero, so only a local term 
 remains. 
In a similar way, the time-domain generalization of the transparency BC (13) is 
derived which grants the absence of spurious reflections a from the artificial computational 
boundary . Applying to Eq. (13) the Fourier transform (29), (32) and denoting maxzz = ipk = , 
we get  
( ) ( )
( ) ∫∫∫
maxmax
  (35) 
Finally, evaluating the inner integral we obtain an elegant 2D boundary condition 
 ( ) ( )
( )( )ηξ
ηξπ −−∂∂
∫ ∫ sx
,, max
max     (36) 
symmetric with respect to the variables x, z , which could be expected from the symmetry of 
the TDPE (24). 
 A simulated example of ultrawide-band EM pulse propagation over a nonuniform 
earth surface with soil parameters 10=ε , )01.0(109 17 мSs −⋅=σ  is depicted in Fig. 11а. 
Evolution of the spatial amplitude distribution for a pulsed signal generated by a Gaussian 
source: ( )[ ]20200 exp)( wzzzA −−=  with a skew curved wave front: 
)(2)()( 00
00 zzzzz −+−=Φ βρ  is shown in a grey color scale. The initial pulse waveform 
 is a damped sinusoid (21) with )(ctf ba =  and spatial length ma 30==Λ π .  
   (a)      (b) 
Fig. 11. Propagation of ultrawide-band pulse (21) over nonuniform earth surface. 
           Initial Gaussian beam parameters: ,80,3000 mmz == 0w 1.0,3000 −== βρ m . 
           Soil conductivity: mS01.0=σ  (a), 0=σ  (b). 
Fig. 12. Received pulse waveform depending on soil conductivity: 
                mS01.0=σ  (solid line), mS001.0=σ  (dashed line); initial pulse (dots). 
The snapshots clearly show the reflected pulse generation at the earth surface. The 
transparency BC (36) imposed at the height mz 500max =  assures unimpeded radiation exit 
from the computational domain. Finite soil conductivity causes signal dispersion appearing in 
a certain delay of the reflected pulse. It is obvious from the comparison with a similar plot 
calculated for a model non-conducting soil 0=σ  - Fig. 11b. A quantitative estimation of the 
effect can be made by means of Fig. 12 revealing a considerable dependence of the pulse 
waveform on the soil conductivity. This effect caused by the pulse penetration into the ground 
can be used for ecological monitoring (water pollution, earthquake precursors, etc.). 
5. Hybrid TDPE and short high-frequency pulse propagation 
 An important practical issue is overland propagation of short EM pulses with high- 
frequency carrier. Basically, having absolute stability, TDPE (24) is capable to describe wide-
band radio pulse propagation. However, the computational expense is drastically growing 
with increasing carrier frequency. If 0ω  considerably exceeds the spectral band of the signal it 
is useful to factor out the carrier: 
 [ ])(exp),,(),,( 00 txkiszxUszx ω−=Π      (37) 
and to consider the transient signal envelope  satisfying a hybrid equation
),,( szxU
ik ,        (38) 
combining the features of standard Leontovich PE (4) with Claerbout TDPE (24). 
Unfortunately, despite a relatively slow variations of  in space-time, straightforward ),,( szxU
numerical solution of this hybrid TDPE (HPE) entails considerable difficulties, as the large 
coefficient  by 0k xU ∂∂  demands a dense computational grid. On the other hand, the 
presence of a big parameter allows us to construct an asymptotic solution of HPE (38) 
radically reducing the computational burden. 
 In order to find a proper asymptotic Ansatz, consider monochromatic PE (4) at the 
carrier frequency ck00 =ω with the initial condition  . We admit 
complex values of the eikonal  to describe relatively narrow wave packets, like a 
Gaussian beam (6). For , we obtain an asymptotic solution 
00)(),0( zikezAzu Φ=
)(0 zΦ
         (39) ),(0 0),(),,(
zxikezxAkzxu Φ=
where  satisfies a “parabolic” eikonal equation ),( zxΦ
1 2 =Φ+Φ zx          (40) 
while slowly varying amplitude  is governed by the paraxial transport equation ),( zxA
=Φ+Φ+ AAA zzzzx         (41) 
Eqs. (40)-(41), being an approximate form of the well-known laws of geometric optics (GO) 
[4], can be easily solved by the method of characteristics. Consider a particular solution of the 
eikonal equation (40) corresponding to a bundle of rays spreading from a central point 
: 0,0 zzx ==
zx −−+≈
=Φ 20
),(       (42) 
At a characteristic line xzz γ+= 0  we have xxx
),( γγ =Φ . The envelope of the family 
(42) solves the boundary value problem with an arbitrary initial condition . 
Define 
)(),0( 0 zz Φ=Φ
 [ ] xzzxzzx )(
)()(, 0
0000 γγ +Φ=+Φ       (43) 
By differentiating Eq. (43) with respect to x  and  we obtain 0z
 xx zzx ')'1(,2
2 γγγγγ +Φ′=Φ+=Φ+Φ       (44) 
Function (43) will satisfy the eikonal equation (40) if the ray direction is matched with the 
local wave front tilt: )()( 000 zz Φ′=γ . Having constructed the eikonal ),( zxΦ  we reduce the 
transport equation (41) to an ODE  
( ) ( ) 0,
0 =+′+
++ xzxA
γ      (45) 
with an evident integral 
xzzxA
=+      (46) 
In a similar way, an asymptotic solution of the modified Claerbout equation (38) can 
be found. Substituting the Ansatz [ ]),(exp),,(),,( 0 zxikszxIszxU Φ=  into Eq. (38) we obtain 
0 =−+⎟
⎛ Φ+Φ+Φ−+⎟
⎛ +Φ+Φ− xszzzzzzsxxx IIIIIIikikk   (47) 
The leading term ( )20kO  disappears in virtue of the eikonal equation (40). Thus, to the 
accuracy ( )10−kO , a space-time transport equation arises for the slowly varying amplitude 
: ),,( szxI
=Φ+Φ+Φ− IIII zzzzsxx      (48) 
As the coefficients of Eq. (48) do not depend on s, it has a solution of the following form  
        (49) [ ),(),(),,( zxsgzxAszxI Ψ−= ]
Here,  is a solution of the stationary transport equation (41) while  is an arbitrary 
function of 
),( zxA )(sg
xcts −= , and  satisfies a linear PDE ),( zxΨ
 xzzx Φ−=ΨΦ+Ψ          (50) 
Solving Eq. (50) by characteristics one easily gets  
 [ ] [ ] )()(,)(, 00000 zxzzxxzzx θγγ ++Φ=+Ψ       (51) 
where )( 0zθ is an arbitrary function. So, the solution of the HPE (38) has asymptotic 
representation 
 [ ] [ ] [ ] ),,()(),(),(exp),(),(~),,( 000
kzxuzzxsgzxikzxsgzxAszxU
θ−Φ−=ΦΨ−
 (52) 
Here,  is a solution of the standard Leontovich PE (4),  and  are its 
amplitude and eikonal, respectively;  and 
),,( 0kzxu ),( zxA ),( zxΦ
)(sg )( 0zθ  are arbitrary functions, and  is to 
be found from the transcendental equation 
),(0 zxz
zxzz =Φ′+ )( 000 . Asymptotic solution (52) is a 
paraxial version of the space-time GO [29], the rays and wave fronts being defined 
numerically via parabolic equation. Functions  and ),( zxA ),( zxΦ  are generally complex-
valued, so distinction between wave amplitude and complex “phase” is made solely on the 
basis of their different dependence on frequency. In particular, complex eikonal ,  
defined as  
),( zxΦ
),,(log),(
0 kkukki
izx     (53) 
is calculated from PE numerical solutions at two close frequencies ck 2,12,1 =ω . 
 Physically, complex eikonal in Eq. (52) appears due to diffraction effects described by 
PE (4). An important consequence is the absence of singularities in the constructed 
asymptotic solution, as the "parabolic" rays do not produce caustics in the real space. Another 
effect caused by diffraction – pulse envelope distortion also is taken into account via complex 
values of the signal delay czx ),(ψ . Physical meaning of complex xcts −=  is provided by 
the theory of analytic signal [22]. 
 Arbitrary functions in Eq. (52) are determined by the initial conditions. In the simplest 
case the constructed transient (37) has the form 
 [ ]),(),(),,(),,( zxxctFzxAxctzxtzxH Φ−−≈−Π≡ +     (54) 
where  and  are complex amplitude and eikonal evolved from the initial 
,  given by Eq. (23) and  is the analytic signal (19) 
corresponding to the real signal . Physical meaning has the real part of the complex 
solution (54) or, from the practical point of view, its normalized envelope 
),( zxA ),( zxΦ
)(0 zA )(0 zΦ )(exp)()( 0 tictftcF ω−≈
)(ctf
2),,( szxH . 
 In virtue of the superposition principle, a more general asymptotic solution can be 
constructed as a number of terms (54). That is a direct analogy with ordinary GO where the 
incident and reflected waves correspond to different ray families. An important practical 
example is radar pulse propagation over the earth surface when the direct and reflected from 
the ground pulses can be distinguished and used for target location [23].  
 Consider first a model example: a short pulse with carrier frequency  
and damped sinusoidal envelope (21), propagating over a slowly rolling boundary . 
Initial pulse parameters are: 
МHz2000 =f
)(xhz =
maba 9, ≈=Λ= π ; ,150 m=w  ,01.0−=β  m2000 =ρ . 
Stationary field distribution calculated by numerical integration of PE  (4) at cfk 00 2π=  
produces a regular interference pattern – Fig. 13. It can be represented as a superposition of  
  Fig. 13. Monochromatic attenuation function.    Fig. 14. Reflected wave phase distribution. 
the incident Gaussian beam (6)  with the reflected wave ),,( 0kzxui ir uukzxu −≡),,( 0  
determined by the terrain  and the impedance BC (11). Functions  
are given by the asymptotic solution (39)-(41), and eikonal 
)(xhz = ),(),,( zxzxA ii Φ
),( zxrΦ  is reconstructed from the 
spatial phase distribution of the reflected wave – see Fig. 14. In accordance with such a 
monochromatic framework, a two-term asymptotic formula arises for the pulsed transient:  
[ ] [ ]
[ ] [ ),(),(),(),(
),(),(),(),(~),,(
),(),( zxxctFzxeuzxxctFezxu
zxxctFzxAzxxctFzxAtzxH
ri Φ−−+Φ−−=
Φ−−+Φ−−
+Φ−+Φ−
 (55) 
Note that to find the amplitude and complex delay of the incident and reflected signals we 
need just to solve the standard PE in frequency domain at two close frequencies 02,1 ωω ≈  - 
see Eq. (53). In time domain, calculation reduces to the evaluation of an analytic function 
 for the given argument values of interest. That radically diminishes the required 
computational resources compared with direct numerical integration of the TDPE (24). 
)(sF +
Fig. 15. Received pulse envelope as function of receiver altitude. 
         (a)       (b) 
Fig. 16. Comparison between numerical solution of TDPE (24) (a) and asymptotic  
             solution (55)  (b). 
The envelope of the received analytic signal (55), as a function of xcts −=  and the 
receiver height z, for  is shown in Fig. 15. One can see profound interference 
minima near the earth surface and a good separation of the direct and reflected pulse to 
heights above 400 m. Figures 16a,b compare the asymptotic solution with direct numerical 
integration of the Claerbout TDPE (24). Qualitatively, they are almost identical. Some hardly 
seen discrepancy is due to a limited accuracy of the asymptotic solution and FD scheme (27). 
This comparison demonstrates the efficiency of the developed approach. Substantial 
acceleration of the numerical procedure (by around 200 times in this example) makes a good 
reason to use it in realistic conditions. 
kmx 7=
As an example, we simulate an experimental situation [23]: radar pulse propagation 
between two aircrafts flying by parallel routs over an irregular terrain. The experiment [23] 
was aimed at simultaneous determination of the target range and altitude from the measured 
return times of the direct radar pulse and the echo signal from the earth surface. Our goal is to 
develop an efficient method of EM field calculations under conditions of multipath and signal 
distortion. At such large ranges )100( kmX =  the Earth sphericity must be taken into 
account. For this purpose, a parabolic hump *2)( earthRxXx −  has been added to the real terrain 
profile plotted in .9 of [23]. Atmospheric refraction has been considered by using the 
equivalent Earth radius earthearth RR 34
* =  [12]. Global field strength distribution produced by 
the incident carrier wave at  is depicted in Fig. 17a. Despite evident multipath 
character of the reflected wave (Fig. 17b) its eikonal 
МHz1410 =f
),( zxrΦ  has a rather regular structure. 
Therefore, our PE based version of complex GO can be applied to simulate the averaged  
       (a)       (b) 
Fig. 17. Simulation of experiment [23]: global field strength (attenuation function) at 
, reflected wave (b). )(1410 aMHzf =
parameters of the received radar pulse (the actually observed signal is a stochastic quantity 
with normal distribution [23]). Its envelope, as a function of the distance from the paraxial 
wave front xcts −=  and the receiver height z, is depicted in Fig. 18.  
Fig. 18. Gaussian pulse envelope as function Fig. 19. Target height from reflected  
            of relative delay and receiver altitude.   pulse delay [23]. 
A Gaussian pulse waveform is chosen with the parameters corresponding to the experimental 
situation [23]:  ,141,3.5 00 МHz== fkmz m75≈Λ . The direct and reflected pulses are 
distinctly separated for  which allows one to reliably solve the triangle  
for target altitude determination - see Fig. 19 borrowed from [23]. Variability and statistics of 
kmz 5.4> 21,, RRR
the simulated reflected pulse resemble the experimental plots presented in [23], and the 
calculated received pulsed signal envelope for a fixed receiver height kmz 3.5=  (Fig. 20)  
  Fig. 20. Received pulse envelope: asymptotic HPE solution (solid line), 
                                     experiment [23] (dash), stochastic model calculation [23] (dots). 
agrees well both with the experimental data and the results of thorough statistical modeling 
[23]. The quantitative discrepancy in the reflected pulse amplitude does not exceed the 
inherent uncertainty due to the errors in terrain description.  
Acknowledgements 
This work was supported in part by a joint RFBR-DFG grant No 01-02-04003. The 
authors are grateful to Friedrich Landstorfer and Ningyan Zhu who initiated this research. 
References 
1. M.A. Leontovich, A new method to solve problems of EM wave propagation over the 
earth surface (in Russian). USSR Academy of Sciences Trans., Physics Series, vol. 8, 
No. 1, pp. 16-22, 1944. 
2. M.A. Leontovich, V.A. Fock, Solution of the problem of electromagnetic wave 
propagation along the Earth’s surface by the method of parabolic equation, J. Phus. 
USSR, vol. 10, pp. 13-23, 1946. 
3. G.D. Malyuzhinets, Progress in understanding diffraction phenomena, Soviet. Phys. 
Uspekhi., vol. 69, pp.321-334, 1959. 
4. V.M. Babič, V.S Buldyrev. Short-Wavelength Diffraction Theory: Asymptotic 
Methods. Springer-Verlag, Berlin, 1990 
5. D. Lee, A.D. Pierce, E.-C. Shang, Parabolic equation development in the twentieth 
century, J. Comput. Acoustics, vol. 8, No. 4, 527-637, 2000. 
6. L. A. Vainstein, Open resonators and open waveguides (in Russian). Soviet Radio, 
Moscow, 1966. 
7. M.D. Feit, J.A. Fleck, Jr. Light propagation in graded-index fibers. Appl. Optics, vol. 
17, pp. 3990-3998, 1978. 
8. Yu.V. Kopylov, A.V. Popov, A.V. Vinogradov. Application of the parabolic wave 
equation to X-ray diffraction optics, Optics Communications, v. 118, pp. 619-636, 
1995. 
9.  G.D. Malyuzhinets, A.V. Popov, Yu.N. Cherkashin. On the development of a 
computational method in diffraction theory. 3rd All-Union Symposium on Diffraction 
of Waves, Academy of Sciences, Moscow, 1964. 
10. F.D. Tappert. The parabolic approximation method. Lecture Notes in Physics, v.70: 
Wave propagation and Underwater Acoustics, pp.224-287, Springer-Verlag, New 
York, 1977. 
11. M.A. Leontovich. On approximate boundary conditions on the surface of well-
conducting bodies (in Russian). Investigations on Radiowave Propagation, part 2, pp. 
5-12. Academy of Sciences, Moscow, 1948 
12. V.A. Fock, Electromagnetic Diffraction and Propagation Problems, Pergamon Press, 
1965. 
13. J. F. Claerbout. Fundamentals of Geophysical Data Processing with Applications to 
Petroleum Prospecting, McGraw-Hill, New York, 1976. 
14. A.V.Popov, S.A. Hosiosky. On a generalized parabolic equation of diffraction theory 
(in Russian). J. Comp. Math. and Math. Phys., vol.17, No. 2, pp. 527-533, 1977. 
15. E.A. Polyansky, On the relation between solutions of Helmholtz and SchrÖdinger type 
equations (in Russian), J. Comp. Math. Math. Phys., vol.2, No.1, pp. 241-249, 1972. 
16. M.F. Levy. Parabolic equation modelling of propagation over irregular terrain. 
Electronics Letters, vol. 26, pp. 1153-1155, 1990. 
17. V.A. Baskakov, A.V. Popov. Implementation of transparent boundaries for numerical 
solution of the Schroedinger equation, Wave Motion, vol.14, No.1, pp.123-128, 1991. 
18. S.W. Marcus. A generalized impedance method for application of the parabolic 
approximation to underwater acoustics, J. Acoust. Soc. Am., vol. 89, pp.391-398, 
1991. 
19. M.F. Levy. Parabolic Equation Method for Electromagnetic Wave Propagation. IEE 
Electromagnetic Wave Series vol. 45, 2000. 
20. M.D. Collins. The time-domain solution of the wide-angle parabolic equation 
including the effect of sediment dispersion. J. Acoust. Soc. Am., vol. 84, No. 6, pp. 
2114-2125, 1988. 
21. L.A. Vainstein, D.E. Vakman. Frequency Discrimination in Oscillation and Wave 
Theory (in Russian). Nauka, Moscow, 1983. 
22. E. Heyman, L.B. Felsen. Gaussian beam and pulsed-beam dynamics: complex-source 
and complex-spectrum formulations within and beyond paraxial asymptotics. J. Opt. 
Soc. Am. A, vol. 18, No. 7, pp. 1588-1611, 2001. 
23. L.M. Zurk. Experimental observation and statistics of multipath from terrain with 
application to overland height finding. IEEE Trans. Antennas Propag., vol. 47, No. 1, 
pp. 121-131, 1999. 
24. A.V. Popov, V.V. Kopeikin, Ning Yan Zhu, F.M. Landstorfer. Modelling EM 
transient propagation over irregular dispersive boundary. Electronics Letters, vol. 38, 
No. 14, pp. 691-692, 2002. 
25. A.V. Popov, V.V. Kopeikin, F.M. Landstorfer. Full-wave simulation of overland radar 
pulse propagation. Electronics Letters, vol. 39, No. 6, pp. 550-552, 2003. 
26. A.V. Popov, V.V. Kopeikin. Prediction of electromagnetic pulse propagation along 
the earth surface (in Russian). Progress of Modern Radioelectronics, No. 1’2005, pp. 
20-35. Radiotechnika, Moscow, 2005. 
27. A.V. Popov. Solution of parabolic equation of diffraction theory by finite difference 
method. J. Comp. Math. and Math. Phys., vol.8, No. 5, pp. 1140-1144, 1968. 
28. A.V. Popov. Accurate modeling of transparent boundaries in quasi-optics. Radio 
Science, vol. 31, No. 6, pp. 1781-1790, 1996. 
29. V.M. Babič, V.S. Buldyrev, I.A. Molotkov. Space-Time Ray Method. Linear and 
Nonlinear Waves (in Russian). SPB University Press, St. Petersburg, 1995. 
