IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 68, NO.

8, AUGUST 2020 6287

Effect of the Zenith Angle on Optical Wave

Propagation in Anisotropic Non-Kolmogorov
Atmospheric Turbulence: A New
Experiment-Based Model
Elad Dakar, Ephim Golbraikh , Natan S. Kopeika, and Arkadi Zilberman

Abstract— The spectral properties of the structural function of of the index of refraction fluctuations, a basic assumption was
atmospheric turbulence can differ in the vertical and horizontal suggested that the refractivity fluctuations are flattened in the
directions. Taking into account the influence of this anisotropy on vertical direction by a factor which is constant for all scales
the signal along the direction of propagation at an angle γ , to the
horizon, can be useful in modeling and interpreting the results but is allowed to vary with altitude.
of measurements. Usually when modeling the propagation of an Then, the spectrum  on each isosurface in Fourier space
electromagnetic signal in the anisotropic turbulent atmosphere, is taken as
the effective anisotropic parameter η is introduced. In the latter
η2 · Cn2 μ
case,the dependence of the wave vector is given by the expression n (q) = ·q (1)
k = ηk2ρ +k 2z , where k z is a vertical component of k, and kρ is 4π
horizontal. However, the spectral index of the structural function, where Cn2 is the refractive index structure constant of the
α, remains constant and independent of γ . Here, we consider the spectrum, μ is the parameter that determines the power law
structural function of atmospheric turbulence with dependence of the spectrum, and q2 = η2 kρ2 + kz2 is an anisotropic
on the propagation angle γ . The existing experimental data show
that the spectral index α is a function of height and γ . The function of two wavenumbers, one vertical (kz ) and one
dependence of the α parameter on γ is introduced, and it is horizontal (kρ = kx2 + ky2 ). The same value of the exponent
determined from the experimental measurements of atmospheric 
turbulent spectra. As an example of using the proposed approach, of k ≡ kx2 + ky2 + kz2 in the isotropic spectrum.
we present the results of calculations of the log-amplitude
The fact that the horizontal wavenumber used by those
variance and propagation of a Gaussian beam through the
anisotropic turbulent atmosphere. models is kρ which comes from the fact that all of those
models are axisymmetric around the vertical axis. η is a
Index Terms— Atmospheric propagation, electromagnetic wave
propagation, optical propagation in anisotropic media.
parameter that determines the level of anisotropy of the
model and usually is selected between 2.5 and 50 in different
I. I NTRODUCTION models.
In various experiments on the study of the atmospheric

O NE of the important problems in the propagation of

electromagnetic waves in the Earth’s atmosphere is asso-
ciated with the influence of its turbulence on this process [1].
turbulence, 1-D spectra are measured in selected directions
(for example, [11]–[19] and references therein). As follows
from these experiments, the spectrum power (the structure
In this case, account must be taken of the fact that atmospheric function exponent) may depend on the zenith angle of obser-
turbulence may be anisotropic [2]–[8], which can significantly vation. Then, in contrast to the previous anisotropic models,
affect the characteristics of the received signals or images. in the proposed model, the turbulent structure function in the
To consider the anisotropy of atmospheric turbulence, it was stratified atmosphere can be described as
customary to use various versions of the model proposed 
in [9]. In all of these models (see [10], [7]), in the context Dn (ρ, z) = C̃n2 · ( ρ 2 + z2 )α(γ ) (2)
Manuscript received January 9, 2020; revised March 5, 2020; accepted where C̃n2 is the generalized structure constant, α is the power
April 1, 2020. Date of publication April 15, 2020; date of current version of the structure function which changes according to γ —the
August 4, 2020. (Corresponding author: Ephim Golbraikh.)
Elad Dakar and Natan S. Kopeika are with the Electrooptical and Photonics angle of a direction in space to the horizon (either from above
Engineering Department, School of Electrical and Computer Engineering, or below the horizon) [5], [12], [13], [15], [16], [20] which,
Ben-Gurion University of the Negev, Beer-Sheva 84105, Israel (e-mail: in the context of propagation of electromagnetic waves in–
[email protected]; [email protected]).
Ephim Golbraikh is with the Physics Department, Ben-Gurion University the Earth’s atmosphere, is the wave propagation angle to the
of the Negev, Beer-Sheva 84105, Israel (e-mail: [email protected]). horizon. Note that the alpha angle γ is the supplementary
Arkadi Zilberman is with the Katif Research and Development Center, angle to the local zenith angle. In the Kolmogorov, isotropic
Ministry of Science, Sedot Negev Regional Council, Netivot 842500, Israel
(e-mail: [email protected]). case, α does not depend on the angle and is equal to 2/3.
Digital Object Identifier 10.1109/TAP.2020.2986720 However, in the general case, α could be a function of γ .
0018-926X © 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
See https://www.ieee.org/publications/rights/index.html for more information.

Authorized licensed use limited to: Ben-Gurion University of the Negev. Downloaded on August 05,2020 at 05:12:17 UTC from IEEE Xplore. Restrictions apply.
6288 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 68, NO. 8, AUGUST 2020

depending on the angle of propagation for spherical and plane

waves and compare the obtained results with existing models.
Besides, we investigate the behavior of a Gaussian beam
propagating vertically in the atmosphere.
II. M ODEL
The structural function (2) is 1-D homogeneous (homogene-
ity is equal to the spectral index α(γ ) and depends on the angle
γ ) but is not an isotropic function. Moreover, in the plane
perpendicular to the vertical, we consider the spectrum to be
isotropic. Within the framework of these assumptions, a 3-D
Fig. 1. Solid angle γ , defined by the propagation direction and the revolution
structural function will be determined.
of the plane angle b around the z axis. The area of the solid angle within the Following [5], the structure function spectrum has the form
cross section is denoted by dot-dashed line. in the three subranges:
 L0  L0
1
n (kz , kρ ) ≈ Dn (ρ, z) · sin(kρ ρ)
2π 3 kρ 0 −L 0
We believe that this assumption is the most physical one, and ·ex p(−i kzz)dzdρ (3a)
moreover, such an approach is used for the description of the
turbulent anisotropy in plasma as in [21]. Hence, for γ = 0, in the first and second subranges and
◦  L0  L0
the direction is horizontal, and for γ = 90 , the direction 1
is vertical, and for angles below the horizontal direction, Dn (ρ, z)· J0 (kρ ρ)·ex p(−i kzz)ρdzdρ (3b)
(2π)2 0 −L 0
we assume that α(−γ ) = α(γ ).
in the third subrange.
Fig. 1 shows that the angle, γ , determines the direction of
After all the integrations, we obtain for n (kz , kρ ) in
the electromagnetic wave propagation. Besides, a cone with
different subranges
angle, b, of the wave propagation is shown, which is much
First subrange
smaller than γ , that is, |γ ± (b/2)| ≈ γ .
1
This somewhat artificial fact is justified by the need to keep ≤ α(γ ) < 1
3
the axial symmetry in going from the structure function to its    1−α(γ ) π   2 2
C̃n2 a0 α α(γ ) sin [ 2 ]
k a
spectrum. exp − z4 0
As mentioned in the title of this article, the model is → × (4a)
2π 2.5 kρ(α+2)
experiment based, indicating that the model is built in accor-
1 < α(γ ) < 2
dance with experimental measurements of turbulent spec-    α(γ )−2]π 
tra in atmospheric conditions at various heights in different C̃ 2 a α(α − 1) α(γ ) − 1 cos [
n 0 2
directions. As noted above, it follows from the experimental →−
 2 2 2π 2.5
data that with increasing height, the measured spectral index k a
changes. exp − z4 0
× . (4b)
In order to be able to calculate the parameters regarding the kρ(α+2)
communication channels in the turbulent atmosphere such as
Second subrange
the log-amplitude variance σx2 , we use the structure function
1
spectrum in Fourier space. In this case, two crucial physical ≤ α(γ ) < 1
assumptions were used. We suggest that the correlations along 3    

the propagation direction, and the directions close to it, affect C̃n2 2(α+1) a0 α α(γ ) sin 1 − α(γ ) π2
→ √ α
the wave behavior significantly more than the correlations 2π 2.5 kρ 2
along any other direction. This assumption allows us to define   2 
1 kρ − |kz |
a solid angle which includes all points whose line to the origin ×
(α+1) exp − a0
is at an angle with the horizon, which equals approximately the |kz | + kρ 4

 2 
angle γ of the propagation direction. The second assumption sign |kz | − kρ kρ + |kz |
is that we suggest that dα(γ )/dγ is small enough such that, −  
exp − a0 (5a)
|kz | − kρ  (α+1) 4
in all integrals, the value of α(γ ) inside the solid angle
can be taken as constant. In order to be able to use these 1 < α(γ )< 2
   

assumptions, we divide the range of propagation directions C̃n2 2(α+1) a0 α(α − 1) α(γ ) − 1 cos α(γ ) − 2 π2
→ √ α
into three subranges 1) π/6 ≥ γ ≥ 0; 2) π/3 ≥ γ > π/6; and 2π 2.5 kρ 2
3) π/2 ≥ γ > π/3. 
 2 
sign |kz | − kρ kρ + |kz |
In this article, we will compare the calculations of some ×  
exp − a0
parameters of electromagnetic wave propagation through the |kz | − kρ  (α+1) 4
atmosphere obtained for various approximations of anisotropic  2 
1 kρ − |kz |
turbulence. Within the framework of the proposed model, −
(α+1) exp − a0 . (5b)
we consider the behavior of the log amplitude variance |kz | + kρ 4

Authorized licensed use limited to: Ben-Gurion University of the Negev. Downloaded on August 05,2020 at 05:12:17 UTC from IEEE Xplore. Restrictions apply.
DAKAR et al.: EFFECT OF THE ZENITH ANGLE ON OPTICAL WAVE PROPAGATION 6289

Third subrange spectrum shows that the said function on the perpendicular
plane is used as part of a product, which is being integrated
1
≤ α(γ ) < 1 over the said plane. When expressed by the new variables,
3
   1−α(γ ) π   k2 a 2  the infinitesimal area on the plane for an anisotropic spectrum
C̃n2 a02 α α(γ ) cos [ 2 ] exp − ρ4 0 is dk x̃ ·dk ỹ , so the power law for the effective integrand remains
→− · (6a) −(α+2). However, for an isotropic spectrum, the infinitesimal
4π 2 kz(α+1)
area on the plane is 2πk ·dk, so the power law for the effective
1 < α(γ ) < 2
   2−α(γ ) π  integrand changes from −(α + 3) to −(α + 2). Thus, the two
C̃n2 a02 α(α − 1) α(γ ) − 1 sin [ 2 ] power laws are effectively the same. In the second range, one
→− can see that the spectrum is a product of power functions with
 k2 a 2  4π 2
different bases. However, in this region, kz and kρ are of the
exp − ρ4 0 same order so that, as a first approximation, we can assume
· . (6b)
kz(α+1) that the spectrum is determined by one base. The result is
that the power law of this range is also −(α + 2). In the
In the first range, the expressions are true only for z which third range, where ρ z [5], the dependence on kρ is in
is small in comparison to ρ, and the integrand in (3) must short wavelengths and has the form exp(−(kρ2 a02 /4)). Only for
be multiplied by exp(−(z2 /a02 )), where a0 is a parameter ◦
the vertical direction (i.e., γ = 90 ), does the third range’s
that effectively determines the angle b which is shown in spectrum become cylindrically symmetric, so kρ becomes the
Fig. 1. This is because similar to this calculation, in the single variable of the integration area, while kz is set to the
second range calculation, the integrand is multiplied by constant value of (1/L0 ). The function of kρ in the effective
exp(−((ρ − z)2 /a02 )), and in the third range calculation, it is integrand is kρ · exp(−(kρ2 a02 /4)), so due to the fast decline of
multiplied by exp(−(ρ 2 /a02 )) (see also [5]). The gradient of the exponential function, no power law can be attributed to
(z, ρ − z, r) is perpendicular to the gradient of (ρ, ρ + z, z), the integrand.
respectively. Within each range, (z, ρ − z, r) is very small in This fact is balanced, however, by the multiplication of the
comparison to (ρ, ρ + z, z), respectively. Thus, multiplying the function by the large constant (1/L0 )−(α+1) .
integrand with the abovementioned exponent, one is effectively For the dependence on altitude, we will use the α parameter
limiting the integration to the part of space defined by the angle as a function of the altitude, in the form (see [22] and
b and its revolution around the vertical axis. Within these parts references therein)
of space, each of the “smallness” properties of each range is,  b1
respectively, correct. The larger a0 is, the slower is the descent α2 Hh1
α1 1
of each exponent along the gradient for each variable. The α(h) =  b1 +  b1 ·  b2
result is that the angle b grows as a0 grows. Here, we should 1 + Hh1 1 + Hh1 1 + Hh2
note that, in expressions (7), we corrected the types that were  b2
made in [5] in formula (10). α3 Hh2
In the first subrange, the spectrum is a product of an +  b2 − 3 (7)
exponential function (of −(kz2 a02 /4)) and a power function 1 + Hh2
(of kρ ) with an exponent equal to −(α+2) [where α is defined
by (2)]. The power law of the first range’s spectrum can thus where parameters αi are typical values of the isotropic spec-
be taken as −(α + 2). trum exponents in the respective layer of the atmosphere.
A note is needed regarding the actual usage of a spectrum H1 (H2 ) is the altitude of transition between the layers, and
in communication channel calculations. In many of these b1 and b2 are numerical coefficients.
calculations (e.g., for the log-amplitude variance), the spec- Further, we generalize the three-layer model [22] to condi-
trum is being used at every point along the propagation path, tions when the spectral index may depend on the angle. In the
through its function on a plane in Fourier space, which is anisotropic case, α(h) = α(γ , h)
perpendicular to the propagation direction. In the case of an 
isotropic spectrum, the form of this function is identical to α(γ , h) = αρ2 (h) · cos2 (γ ) + αz2 (h)· sin2 (γ ). (8)
the form of the 3-D spectrum. On the contrary, in the case In this article, we assume that αz (h) equals to α(h) from (7),
of an anisotropic spectrum, the form of the function on the and we have to determine αρ (h). We assume that αρ is
perpendicular plane is determined by the change of variables. determined by the same function from (7) but with other αi .
To do this, we introduce a new coordinate system: k x̃ = k x , In the first-order approximation, for most of the altitudes up to
kỹ sinγ = k y , and kz = k ỹ · cos γ . Then, kρ = the second transition, the turbulence is almost isotropic, so α3
k 2ỹ · sin2 γ + k x̃2 and kz = k ỹ · cos γ , and to avoid a situation is set equal to α1 for αρ (h), with α1 , α2 having the same values
in which kz = k ỹ · cos γ = 0 for γ = (π/2) and thus is a they have for αz (h).
singularity, the actual numerical calculation uses the equation In order to clarify this model’s function of α(γ , h), Fig. 2
kz = k ỹ · cos γ + (1 − cos γ /L 0 ). shows graphs of three α’s as functions of altitude. The
◦
The transition to the new variables does not change the parameters αρ (γ = 0) and αz (γ = 90 ) used in this
power law because the power laws for both kρ (k x̃ , k ỹ ) and article were chosen from [22]. Each graph shows α(γ , h)
kz (k x̃ , k ỹ ) are equal to 1. A further look into the usage of the for a specific direction, and each direction is defined by the

Authorized licensed use limited to: Ben-Gurion University of the Negev. Downloaded on August 05,2020 at 05:12:17 UTC from IEEE Xplore. Restrictions apply.
6290 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 68, NO. 8, AUGUST 2020

III. L OG A MPLITUDE VARIANCE FOR

P LANE AND S PHERICAL WAVES
In [5], we briefly estimated the behavior of the log-
amplitude variance in the dependences between different
heights. Further, we will discuss in more detail the result of
calculations under the new model for different conditions.
The log-amplitude variance equals (see [8])
1  
σx2 (
r , L) = · Re E2 (
r , r
) + E3 (

r , r
) (11a)
2  L  ∞  
E2 (
r1 , r
2 ) = 2πK2 dη d2κ  ˜ n K,
l
Fig. 2. Three α’s of a specific direction versus the altitude h[m]. Legend: 0 −∞
α(γ =00 ) α(γ =450 ) α(γ =900 )  
–-–-– ——— - - - - - - -.
iκ 2

· exp iK · δ
∗
r1 − δ r
2 − ∗
δ − δ (L − l)
2K
 L  ∞  
angle above the horizon γ . It should be noted that, as in E3 (
r1 , r
2 ) = −2πK2 dη d2 κ 
l
˜ n K,
−∞
the works where the anisotropy was described using the  0 
parameter η [for example, (1)], we come to the anisotropic
iκ 2 δ
· exp iδ K · (
r1 − r
2 ) − (L − l) (11b)
structures elongated in the horizontal direction with the chosen K
parameters αρ and αz (see Fig. 3 in [5]).
where 
η) is the spectrum of the covariance of the
˜ n (K,
As one can see in Fig. 2, in the framework of the model
of (7), up to the altitude of 5000 [m] (the lower and middle refractive index fluctuations, and l is a variable taken as
troposphere), there is no noticeable difference between αρ and approximately equal to the propagation variable, so the inte-
αz , which turns α(γ , h) from (9) into a function with almost gration over it goes from the transmitter where l = 0 to the
receiver, where l = L; K
is a 2-D vector in Fourier space
no dependence on γ . Then, the model is effectively isotropic
for altitudes below 5000 [m]. However, it should be noted inside the plane perpendicular to the propagation direction
that generally atmospheric turbulence at these altitudes can be which includes the origin (despite the fact that the integration
over K
does not include the origin of course because each
anisotropic too. If we can determine α1 and α2 for the function

|K| < (1/L0 ) is excluded from the integration). Furthermore,

of αρ , distinct from the α’s for the function of αz , then there
is no problem to generalize the model for all the heights.
def
κ = |K|;
δ is the path amplitude ratio characterizing the
Another turbulent parameter is the outer scale, L0 . It should propagating wave, where for a plane wave δ = 1 and for
be noted that all the spectra of the structural function con- a spherical wave δ = (η/L).
sidered in this article (and not in it only) diverge on the The dependence of 
η) on η comes from the fact
˜ nn (K,
large scales (or k 1). For this, in the theory of turbulence, that the values of α and C̃n change with altitude, and α also
2

the outer scale is introduced, above which the structural changes with the angle between the propagation direction and
function rapidly decreases. Therefore, in all integrals, inte- the horizon at each point along the propagation path.
gration is in scales from zero to L0 , which can be determined Combining (12) together, one gets the result for a real δ
from [11], [23]  L  ∞
2
r , L) = πK ·
σx (
2
dl d2 κ 
l)
˜ n (K,
4
L0 (h) =  h−8500 2 [m]. (9) 0 −∞
 
1+ L−l
2500 · 1 − cos κ 2 · ·δ (12)
K
For calculations involving C̃n2 from (2), we will use the
H − V5/7 turbulence model [8], [22] in which the “normal- which makes clear the importance of δ in the resulting
ized” C2n is σx2 (
r , L).
The log-amplitude variances were calculated for various
8
π
 (α− 23 ) initial angles γ0 of the propagation direction γ , either above
· sin K for an uplink channel or below for a downlink channel of the
C̃n2 = C2n 3
3
(10)
(α + 2)· sin α · π2 L horizon as defined at the location of the transmitter (the slant-
path height trough atmosphere decreases from left to right in
where K is the wavenumber of the electromagnetic wave,√ and the figures). The calculations were performed for wavelength
L is the length of the propagation path. We see that (K/L) λ = 1.55 [μm] and spherical geometry. In the spherical geom-
functions as the normalizing scale for z and ρ in the structure etry approximation, one calculates the propagation direction
function from (2). angle to the horizon for each point along the propagation path
On the other hand, under real conditions in the non- while modeling Earth and its atmosphere as a sphere. In this
Kolmogorov case, the parameter α(γ , h) in expressions— case, angle γ changes with an altitude. The horizon for each
(7) and 10 may be ≈1 or ≈2. In these cases, we use the point (altitude) is thus taken as the perpendicular line to the
values 0.99 or 1.99 for calculations. radius, from the center of the Earth to the point on the path.

Authorized licensed use limited to: Ben-Gurion University of the Negev. Downloaded on August 05,2020 at 05:12:17 UTC from IEEE Xplore. Restrictions apply.
DAKAR et al.: EFFECT OF THE ZENITH ANGLE ON OPTICAL WAVE PROPAGATION 6291

Further, we describe a common spectrum applied in the

calculations for the entire range of angles. The determination
of the values for each uplink/downlink channel was done by
◦ ◦
first finding the value for a0 for the subrange 80 ≤ γ0 < 90
so that the result for σx2 will be similar to the one in [24]. This
value and all the other values were set as a factor multiplying
the parameter L0 . The next step was to determine the value
◦ ◦
for the subrange 60 ≤ γ0 < 80 so that the gap in the
value of σx2 in the transition between the two subranges will
be minimal. Another important consideration was to avoid a
situation, in which both the graphs near the transition point
are descending (ascending) while transition itself is ascending
(descending). The determination for ranges with smaller γ0
was performed using the same two considerations, where the
smallest γ0 was taken to be 10◦ .
The single spectrum is a kind of a unification of the spectra
from the three ranges. For this unification, we multiply each
term from each range’s spectrum with an exponential function
of the angle. The said exponential function decreases the
further angle from the (term’s) respective range of angles. For
the second range’s spectrum, which includes two terms, one
takes only the term, which decreases more slowly further away
from the Fourier space origin. All the term are taken without
their coefficients, so there is, now, a difference in the unified
spectrum whether 1/3 ≤ α(γ ) < 1 or 1 < α(γ ) < 2. Then,
the general form of the unified spectrum is thus

n k z , k ρ , γ
⎡  2 ⎤
(1)

k a Log-amplitude variance versus γ0 for the spherical waves.
⎢ z 0 ⎥ −(α+2) Fig. 3.
= Aex p −v 1 γ 2 ex p⎣− ⎦k ρ Solid curve—present model; Dots—Kolmogorov model; Dash—q-model.
4 (a) Uplink channel; receiver at 15 [km] altitude. (b) Downlink channel;
15 [km] transmitter height.
⎧ !
(2) "2 ⎫
    ⎪
⎨ − ⎪
⎬
π 2 k ρ k z · a0
+ Bex p −v 2 γ − ex p −
4 ⎪
⎩ 16 ⎪
⎭ The parameters v 1 = v 2 = v 3 equal to 1 for all the
calculations. For the uplink channels of the plane wave case,

−(α+1)   
kρ + kz π 2 B = 0.005 for H0 = 15 [km], respectively, and A = C = 0.1.
· + Cex p −v 3 · γ − The parameters (a0(1), a0(2) , a0(3) ) were chosen for the altitude
kρ 2
⎡  2 ⎤ H0 = 15 [km]: the plane wave for the downlink channel—
(3)
kρ ·a0 (6, 3, 4) and uplink channel—(2, 20, 2.6); for the spherical
⎢ ⎥ −(α+1)
· ex p⎣− ⎦k z (13) wave, (6, 11, 9)—downlink channel and (8, 20, 5.5)—uplink
4
channel.
Calculations of the log-amplitude variance versus the initial
where A, B, and C are the coefficients in the single spectrum, angle γ0 involving α(γ , h), which is determined by (7) and is
which fix the dimensions of all terms to be identical and also shown in Fig. 2, are shown in Figs. 3 and 4. On the other
give every term its fitting weight in the spectrum. v 1 , v 2 , and v 3 hand, calculations based on plane geometry, in which the
are the parameters inside the exponential functions of the angle between the propagation direction and the horizon at
angle, which determine the rate, in which the functions any specific point along the propagation path is taken to be
decrease further from their ranges of angles. In this spectrum, equal to the initial angle, were also conducted. The difference
they are numbered to clarify that each one may take a different between the results of the two calculations was found to be
value from the other. less or equal to 10−5 due to the relatively low altitudes of all
In all the calculations of the log-amplitude variance for the the receivers/transmitters, so only the results for the spherical
downlink channels, the values of the parameters are A = B = geometry are presented here.
C = 0.1, except for the case of the spherical wave, in which Fig. 3 also shows the results of calculations based on
B = 0.01 for H0 = 15 [km] transmitter height. the Kolmogorov and q-models [see-(1)]. In the Kolmogorov
The same is true for the uplink channels of the spherical model, we applied 
η) = 0.033 · C2n (h) · |K|
˜ n (K,
−(11/3) with
wave case, except for B = 0.05 for H0 = 15 [km] receiver L0 (h) from (9) and C2n (h) from [23]. At the same time, the
height. q-model (1) was used with η = 30 (one of the values proposed

Authorized licensed use limited to: Ben-Gurion University of the Negev. Downloaded on August 05,2020 at 05:12:17 UTC from IEEE Xplore. Restrictions apply.
6292 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 68, NO. 8, AUGUST 2020

in work [21]) and μ = −11/3. Besides, for some initial γ0 In our calculations, we use the height distribution of the
angles, the calculated value of the log-amplitude reaches the alpha parameter in accordance with (7) and its dependence on
saturation range. Hence, to account for saturation, we use the gamma angle (9).
(sat ) (cal) (cal)
σx2 = 0.25 · ln[2 − ex p(−4 · σx2 )], where σx2 is the The 3-D spectrum of the atmospheric turbulence in [24] was
initial result which we obtain in the numerical calculations, adopted from the anisotropic von-Karman model
(sat )
and σx2 is the result presented in graphs after the saturation for 3 < α̃< 5
has been taken into account. It should be noted that for values 2

α̃
˜ n (κ, α̃) = A(α̃)C̃n2 ζ 2 κ 2 + k02 − 2 exp − κ
(cal) (sat )
far from the saturation range, σx2 ≈ σx2 .  (14)
2
km
As shown in Figs. 3 and 4, the values for the spherical
waves are smaller in comparison to the parallel plane waves. 
It is due to the difference in the path amplitude ratio δ in (12) where κ ≡ ζ 2 kxy 2 + k2, k
z xy ≡ kρ , k0 ≡ (2π/L0 ), km ≡
between the two cases. In all the cases, the graph for the (C(α̃)/l0 ), ζ is the effective anisotropic factor, and
isotropic Kolmogorov calculation is decreasing monotonically.
(α̃ − 1)  π
The q-model, which is decreasing monotonically for all the A(α̃) ≡ cos α̃ (15a)
plane wave cases, has graphs with minimum points in all the 4π 2 2
   α̃−5
1
spherical wave cases. The reason for that q-model behavior 3 α̃ 3 − α̃
C(α̃) = πA(α̃) − . (15b)
connects with the difference in the path amplitude ratio (δ). 2 2 3
The amplitude in the spherical wave case reduces the influence
that the length of propagation has on the value of the log- Since all the calculations in [24] are performed for vertical
amplitude variance. At the same time, the anisotropic model propagation, the Markov approximation [10] implies that the
has a minimum point for all the cases. only relevant spectrum is the 2-D spectrum in the kz = 0
The presence of the minimum is associated with a change plane. Because of that, the spectrum which is being used in
in the value of the spectral index depending on the height and the calculations depends only on kρ . For the same reasons,
the γ , as shown in Fig. 2. the relevant spectrum for vertical propagation path is the
In the plane wave case (for the downlink/uplink chan- spectrum obtained for kz = (1/L0 ) in (7) (because assigning
nels), the value of γ0min , which corresponds to the minimum kz = 0 leads to a singularity). As a result, the relevant
of σx2 , decreases when the altitude of the transmitter/receiver, spectrum also depends only on kρ .
H0 , increases. Using [24] for calculations of the cross-spectral density
For the spherical wave, γ0min does not always decrease as H0 matrix, we write the convolution integral (i, j = x, y)
increases. Parameter δ reduces the influence of the length of 
K2 0 0
0 0

propagation on the value of the log-amplitude variance. At the Wi,j (

r1 , r
2 ; ω) = 0
Wi,j r
1 , r
2 ; ω Ki,j r
1 , r
2 , r
1 , r
2 ; ω
same time, the effect of the δ parameter is manifested in the 4(π z) 2

fact that, for the spherical wave, the value of the log amplitude × d2r
10 d2r
20 (16)
variance for each downlink channel is slightly higher than the
value of the corresponding uplink channel. where the superscript “0” denotes points or functions in the
This is in contrast to the plane wave, where the value for source plane z = 0, and Ki,j (
r10 , r
20 , r
1 , r
2 ; ω) is the Green
each uplink channel is substantially higher than the one for function for the cross-spectral density matrix. The first general
the respective downlink channel. In the case of the spherical result is for this Green function, and it gives
wave, the factor [1 − cos(κ 2 · (L − η/K) · δ)] from (12) is

Ki,j r
10 , r
20 , r
1 , r
2 ; ω
close to zero for both the ends of the propagation path, thus 
2
2
reducing the effect the values of C̃n2 near the starting point of r
1 − r
10 − r
2 − r
20
= exp − i K
the propagation have on the resulting σx2 (
r , L). 2z
Therefore, there are three main factors that determine π 2 K2 z !

2 "
σx2 (
r , L) in the model: the propagation length, path amplitude − r1 −
r2 )2 +(
r1 −
r2 ) r
10 −
r20 + r
10 − r
20
(

ratio δ, and anisotropy. For example, when we use (7) and  ∞3 

(9) for α(γ , h) determination, the propagation length and the · κ 3 n (κ, α̃)dκ (17)
0
anisotropy affect the results in opposing manners, which leads
to the minimum in the dependence of σx2 on γ0 . C̃n (h) →C̃n2
2
 H
1
≡ C̃n2 (h)dh. (18)
IV. G AUSSIAN B EAM H − h0 h0

Further, we consider the vertical propagation of a Gaussian This is the average of C̃n2 (h) set by (10), where h0 and H are
beam through the anisotropic turbulent atmosphere. This prob- the altitudes above ground where the transmitter and receiver
lem was studied in [24] and [26]. Authors in that work are located, respectively. In addition to this, we have the
presented the results of calculations for both fully coherent average value of L0 (h)
and partially coherent beams. In this article, we describe the  H
fully coherent beams only. Besides, in this model, we do not 1
L0 (h) → L0 ≡ L0 (h)dh. (19)
use the anisotropy parameter as was done in [24]. H − h0 h0

Authorized licensed use limited to: Ben-Gurion University of the Negev. Downloaded on August 05,2020 at 05:12:17 UTC from IEEE Xplore. Restrictions apply.
DAKAR et al.: EFFECT OF THE ZENITH ANGLE ON OPTICAL WAVE PROPAGATION 6293

Fig. 5. Normalized on-axis spectral density SN versus the propagation

distance z. (a) λ = 1.55 [μm]. (b) λ = 632.8 [nm]. Here: ζ = 10; xxx
ζ = 2; ζ = 1.5; —ζ = 1; ∗∗∗ Present work.

and
Fig. 4. Log-amplitude variance versus γ0 for the plane waves (the same nota-
 
z2 1 1 2π 2 z3 I
tions as in Fig. 3). Solid curve represents present model; Dots represent Kol- 2ij (z)≡ 1+ 2 2 + +
mogorov model; and Dash represents q-model. (a) Uplink channel; Receiver K σ 4σ 2 δij2 3σ 2
at 15 [km] altitude. (b) Downlink channel; 15 [km] transmitter height.  
z2 1 1 σ 2 2ij z
0 0
Rij (z) ≡ 2 2 + +
K σ 4σ 2 δij2 σ 2 2ij (z) + π z 3I−σ
2 3 2
We assume that Wi,j 0
r
1 , r
2 ; ω has a Gaussian–Schell form  ∞

2
2 
0 0
√ √ r
0 + r
20 I≡ κ
˜ n (κ, α̃)dκ. (22)
Wi,j r
1 , r
2 ; ω = Bij Ii Ij · ex p − 1
0
(20) 0
4σ 2 In the framework of this model, we rewrite the integral (22) as
 ∞  ∞
where Bij are complex coefficients characterizing the electric

field, Ii and Ij are the intensities along the x- and y-axis, I≡ κ n (κ, α̃)dκ → I ≡
˜ ˜ kρ ˜ n kρ , α dkρ . (23)
0 0
respectively, σ is the R.M.S. width of the beam, and δij are
the correlation distances between electric field components i This integral (see [25]) has a form
 2
and j so for a fully coherent beam 2−α̃ A(α̃) 2−α̃ k
I=ζ C̃n k̃m ·β · exp 20
2
1 2(α̃ − 2) km
δij → ∞ <=> → 0.  
δij α̃ k 2
× 2 − , 20 −2k̃04−α̃ (24)
From (16), (17), and (20), we obtain 2 km
√ √
Bij Ii Ij where β ≡ 2(k̃02 − k̃m 2
) + α̃ k̃m
2
, k̃02 ≡ (k02 /ζ 2 ), k̃m
2
≡ (km2
/ζ 2 ),
r1 , r
2 ; ω) =
Wi,j (
exp[A + B] (21)
2ij (z) and (,) is the incomplete upper Gamma function.
However, for I˜, we obtain
where ⎧  (α+1)

(
r1 + r
2 )2 iK r
12 − r
22 ⎪
⎪ α 2 1−α L
A=− 2 2 + ⎪
⎪ ·C̃ · (α) · cos π · 02
8σ ij (z) 2Rij (z) ⎪
⎪ π 2 n
2 ā0
⎪
⎨
     for 0 < α ≤ 1
1 1 1 1 2 2 2 I˜ =  (α+1) (25)
B = + 2 + π K zI 1 + 2 ⎪
⎪ α(α − 1) 2−α L
22ij (z) 4σ 2 δij 3 ij (z) ⎪
⎪ · C̃ 2 · (α − 1) · sin π · 02
⎪ π2
⎪ n
2 ā0
 ⎪
⎩
π 4 K 2 z4 I 2 for 1 < α < 2
− (
r1 − r
2 )2
18σ 2 2ij (z) where ā0 ≡ a0 (L0 ).

Authorized licensed use limited to: Ben-Gurion University of the Negev. Downloaded on August 05,2020 at 05:12:17 UTC from IEEE Xplore. Restrictions apply.
6294 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 68, NO. 8, AUGUST 2020

parameter in the determination of wave vector, is introduced.

However, α remains constant and independent of the angle of
propagation.
In this article, we study the new approach for describing the
anisotropy of turbulence relative to the vertical and horizontal
directions in the atmosphere, through the anisotropy of the
structure function exponent, proposed in [5]. It should be noted
that the approach to the description of the anisotropy of the
atmospheric turbulence is based on atmospheric experimental
data (see [5] and references therein).
Thus, for the calculations of the log amplitude variance for
plane and spherical waves, there are three key factors: the
propagation length, path amplitude ratio, δ, and dependence
of the index of the spectral function α on the angle to the
horizon.
In parallel, we have applied our model for the calculations
of Gaussian beams. The obtained results show that the model
is well founded for the purpose of such calculations, and they
show that the model can be considered as having an effective
anisotropic factor of ζ = 2 in comparison to the previous
model. Since this model is based on physical parameters,
they can be determined from experiments on the optical wave
propagation of laser beams in the atmosphere. Thus, we
Fig. 6. rms beamwidth σ (z) versus the propagation distance z.
a) λ = 1.55 [μm]. b) λ = 632.8 [nm]. Notations as in Fig. 4. believe that future experiments will measure the properties of
communication channels in the atmosphere whose propagation
path is neither vertical nor horizontal. Such measurements will
The spectral density S(
r ; ω) connects with the convolution make it possible to determine the validity of the model, and
integral as to what extent, if any, it should be modified. As for further
& ' theoretical and/or numerical research of the model, the authors
r ; ω) ≡ Trace Wi,j (

S(
r , r
; ω) . (26) plan to use the model in order to calculate the modulation
The normalized on-axis spectral density SN of a fully coherent transfer function associated with image transmission through
beam versus the propagation distance z is shown in Fig. 5. turbulent media, where the turbulent media are described by
Each curve in that figure stands for a different value of ζ for the model.
the spectra (14) from [24], where the value for α̃ is taken Typically, there are different combinations of parameters
as α̃ = 3.5 <=> α = 0.5. All the calculations were made that give the same variance. The systematic variations with
for two wavelengths λ = 1.55 [μm] and λ = 632.8 [nm]. propagation direction have always been difficult to isolate
At the same time in our calculations, we used a fully coherent meaningfully. The spatial variation of the field is the best
beam propagating vertically from h0 = 0 to H = 30 [km], indicator of anisotropy but even that is difficult to interpret
σ = 0.05 [m], l0 = 10−3 [m], and L̄ 0 = 1 [m] and calculated because of the cumulative effects of propagation. Moreover,
C̃n2 = 5.1 · 10−17 [m]3−α̃ . We use L0 as a constant, and ā0 ≡ the limiting values of αρ (h, γ ) and αz (h, γ ) may depend on
a0 (L0 = 1 [m]) = (L0 /3) = (1/3) [m]. the time of day, latitude and longitude, time of year, solar
One can see from the figure that the value of SN obtained activity, etc. Therefore, their determination is a separate and
in this model corresponds more or less to the values obtained not simple task.
in [25] for ζ ≈ 2. Thus, the study using various models of the effect of
the anisotropy of the atmospheric turbulence field on the
 (21) and (22), for S(
r = 0; ω), the expres-
As follows from
sion σ (z) ≡ σ 2 2ij (z) functions as the rms beamwidth, propagation at an angle to the horizon can be important in
modeling and interpreting the results of various measurements.
which changes over the propagation distance z. The function
σ (z) is shown in Fig. 6 (the same notations as in Fig. 5).
R EFERENCES
V. C ONCLUSION [1] R. L. Phillips and L. C. Andrews, Laser Beam Propagation Through
Random Media. Bellingham, WA, USA: SPIE, 2005.
It is known that the spectral properties of the structural func- [2] M. S. Belen’kii, S. J. Karis, C. L. Osmon, J. M. Brown, II, and
tion of atmospheric turbulence can differ in the vertical and R. Q. Fugate, “Experimental evidence of the effects of non-Kolmogorov
turbulence and anisotropy of turbulence,” in Proc. 18th Congr. Int.
horizontal directions. Usually, when modeling the propagation Commission Opt., Jul. 1999, pp. 50–51.
of an electromagnetic signal in a turbulent atmosphere, either [3] M. S. Belen’kii, J. D. Barchers, S. J. Karis, C. L. Osmon,
the Kolmogorov isotropic model with the spectral index of J. M. Brown, II, and R. Q. Fugate, “Preliminary experimental evidence
of anisotropy of turbulence and the effect of non-Kolmogorov turbu-
the structural function α = −5/3 or the anisotropic model lence on wavefront tilt statistics,” Proc. SPIE, vol. 3762, pp. 396–406,
in which the effective anisotropic parameter, η, as a weigh Sep. 1999.

Authorized licensed use limited to: Ben-Gurion University of the Negev. Downloaded on August 05,2020 at 05:12:17 UTC from IEEE Xplore. Restrictions apply.
DAKAR et al.: EFFECT OF THE ZENITH ANGLE ON OPTICAL WAVE PROPAGATION 6295

[4] M. S. Belen’kii, E. Cuellar, K. A. Hughes, and V. A. Rye, “Experimental Elad Dakar is currently pursuing the Ph.D. degree
study of spatial structure of turbulence at Maui Space Surveillance Site with the Faculty of Electrooptic and Photonic
(MSSS),” Proc. SPIE, vol. 6304, Sep. 2006, Art. no. 63040U. Engineering, Faculty of Electric and Computer
[5] E. Dakar, E. Golbraikh, N. Kopeika, and A. Zilberman, “Electromagnetic Engineering, Ben-Gurion University of the Negev,
wave propagation in the turbulent atmosphere with an anisotropic Beer-Sheva, Israel.
exponent of the spectrum,” IEEE Trans. Antennas Propag., vol. 65, He is currently the Data Scientist with the
no. 10, pp. 5654–5657, Oct. 2017. Ben-Gurion University of the Negev.
[6] R. E. Hufnagel and N. R. Stanley, “Modulation transfer function
associated with image transmission through turbulent media,” J. Opt.
Soc. Amer. A, Opt. Image Sci., vol. 54, no. 1, pp. 52–61, 1964.
[7] I. Toselli, “Introducing the concept of anisotropy at different scales for
modeling optical turbulence,” J. Opt. Soc. Amer. A, Opt. Image Sci.,
Ephim Golbraikh received the Ph.D. degree
vol. 31, no. 8, pp. 1868–1875, 2014.
in physical and mathematical sciences from the
[8] A. Zilberman, E. Golbraikh, N. S. Kopeika, A. Virtser, I. Kupershmidt,
Applied Geophysics Institute, USSR, Moscow,
and Y. Shtemler, “Lidar study of aerosol turbulence characteristics in
Russia, in 1988.
the troposphere: Kolmogorov and non-Kolmogorov turbulence,” Atmos.
He is currently Professor with the Physics
Res., vol. 88, no. 1, pp. 66–77, Apr. 2008.
Department, Ben-Gurion University of the Negev,
[9] F. Dalaudier and A. S. Gurvich, “A scalar three-dimensional spec-
Beer-Sheva, Israel. He has coauthored with Alexan-
tral model with variable anisotropy,” J. Geophys. Res., vol. 102,
der Eidelman, Herman Branover, and Semen Moi-
pp. 19449–19460, Aug. 1997.
seev the book Turbulence and Structures published
[10] V. S. R. Gudimetla, R. B. Holmes, and J. F. Riker, “Analytical expres-
by Academic Press, NY, Boston, London, in 1999.
sions for the log-amplitude correlation function for plane wave propaga-
His current research interests include the cover the-
tion in anisotropic non-Kolmogorov refractive turbulence,” J. Opt. Soc.
ory of turbulence in hydrodynamic (HD) and magnetohydrodymanic (MHD)
Amer. A, Opt. Image Sci., vol. 29, no. 12, p. 2622, Dec. 2012.
media, propagation of the electromagnetic waves and sound through the
[11] O. G. Chkhetiani, A. Eidelman, and E. Golbraikh, “Large- and small-
turbulent atmosphere and media.
scale turbulent spectra in MHD and atmospheric flows,” Nonlinear
Processes Geophys., vol. 13, no. 6, pp. 613–620, 2006.
[12] F. Dalaudier, M. Crochet, and C. Sidi, “Direct comparison between Natan S. Kopeika was born in Baltimore
in situ and radar measurements of temperature fluctuation spectra: in 1944 and raised in Philadelphia. He received
A puzzling result,” Radio Sci., vol. 24, no. 3, pp. 311–324, May 1989. the B.Sc., M.Sc., and Ph.D. degrees in electrical
[13] F. Dalaudier, C. Sidi, M. Crochet, and J. Vernin, “Direct evidence of engineering from the University of Pennsylvania,
‘sheets’ in the atmospheric temperature field,” J. Atmos. Sci., vol. 51, Philadelphia, PA, USA, in 1966, 1968, and 1972,
no. 2, pp. 237–248, Jan. 1994. respectively.
[14] E. Golbraikh, H. Branover, N. S. Kopeika, and A. Zilberman, “Non- He and his family moved to Israel, and he began
Kolmogorov atmospheric turbulence and optical signal propagation,” his career with the Ben-Gurion University of the
Nonlinear Processes Geophys., vol. 13, no. 3, pp. 297–301, 2006. Negev as a Lecturer in 1973. He was the Chair
[15] M. Lilley, S. Lovejoy, K. Strawbridge, and D. Schertzer, “23/9 dimen- with the Department of Electrical and Computer
sional anisotropic scaling of passive admixtures using lidar data of Engineering, Ben Gurion University of the Negrv,
aerosols,” Phys. Rev. E, Stat. Phys. Plasmas Fluids Relat. Interdiscip. from 1989 to 1993, and was named Reuven and Francis Feinberg Professor
Top., vol. 70, no. 3, pp. 036307-1–036307-7, Sep. 2004. of electrooptics in 1994. He was the first Head of the Department of
[16] M. Lilley, S. Lovejoy, K. B. Strawbridge, D. Schertzer, and Electrooptical and Photonics Engineering, Ben Gurion University of the
A. Radkevich, “Scaling turbulent atmospheric stratification. II: Spatial Negrv, from 1999 to 2005, which grants graduate degrees in electrooptical
stratification and intermittency from lidar data,” Quart. J. Roy. Meteorol. engineering. He has published about 200 articles in international reviewed
Soc., vol. 134, no. 631, pp. 301–315, Jan. 2008. journals and over 160 articles at various conferences. He has authored the
[17] H. Luce, F. Dalaudier, M. Crochet, and C. Sidi, “Direct comparison textbook A System Engineering Approach to Imaging published by SPIE Press
between in situ and VHF oblique radar measurements of refractive (first printing 1998, second printing 2000). He has coauthored with Natan
index spectra: A new successful attempt,” Radio Sci., vol. 31, no. 6, Blaunstein, Shlomi Arnon, and Arkadi Zilberman the book Applied Aspects of
pp. 1487–1500, Nov. 1996. Optical Communication and Lidar, published recently by CRC Press. Recent
[18] G. D. Nastrom, K. S. Gage, and W. H. Jasperson, “Kinetic energy spec- research involves the development of a novel inexpensive focal plane array
trum of large- and mesoscale atmospheric processes,” Nature, vol. 310, camera for terahertz imaging. Other areas of research include interactions of
pp. 36–38, Jul. 1984. electromagnetic waves with plasmas, the optogalvanic effect, environmental
[19] G. D. Nastrom and K. S. Gage, “A climatology of atmospheric effects on optoelectronic devices, imaging system theory, propagation of light
wavenumber spectra of wind and temperature observed by commercial through the atmosphere, imaging and wireless communication through the
aircraft,” J. Atmos. Sci., vol. 42, pp. 950–960, May 1985. atmosphere, image processing and restoration from atmospheric blur, imaging
[20] A. T. Waterman, T. Z. Hu, P. Czechowsky, and J. Röttger, “Measurement in the presence of motion and vibration, lidar, target acquisition, and image
of anisotropic permittivity structure of upper troposphere with clear-air quality in general. He retired officially in 2013 and is still active.
radar,” Radio Sci., vol. 20, no. 6, pp. 1580–1592, Nov. 1985. Dr. Kopeika is a fellow of SPIE—The International Society for Optical
[21] C. H. K. Chen, R. T. Wicks, T. S. Horbury, and A. A. Schekochihin, Engineering. He and Shlomi Arnon was received the JJ Thomson Award by
“Interpreting power anisotropy measurements in plasma turbulence,” the IEE in 1999 and the Glant Prize for Excellence in teaching in 2001. He is
Astrophys. J., vol. 711, no. 2, pp. L79–L83, Mar. 2010. currently a Topical Editor for Marcel Dekker for the topic atmospheric optics
[22] A. Zilberman, E. Golbraikh, and N. S. Kopeika, “Propagation of elec- in their Encyclopedia of Optical Engineering.
tromagnetic waves in Kolmogorov and non-Kolmogorov atmospheric
turbulence: Three-layer altitude model,” Appl. Opt., vol. 47, no. 34, Arkadi Zilberman received the M.Sc. degree
p. 6385, Dec. 2008. in physics from Tomsk State University, Tomsk,
[23] V. P. Lukin, E. V. Nosov, and B. S. Fortes, “The efficient outer scale of Russia, in 1994, and the M.Sc. and Ph.D. degrees
atmospheric turbulence,” in Proc. ESO Conf. Workshop, vol. 56, 1999, in electrical engineering and electro-optics from the
pp. 521–619. Ben-Gurion University of the Negev, Beer-Sheva,
[24] M. Yao, I. Toselli, and O. Korotkova, “Propagation of electromagnetic Israel, in 2003 and 2007, respectively.
stochastic beams in anisotropic turbulence,” Opt. Express, vol. 22, He is currently a Researcher with the School of
no. 26, pp. 31608–31619, Dec. 2014. Electrical and Computer Engineering, Ben-Gurion
[25] I. Toselli and O. Korotkova, “Spread and wander of a laser beam University of the Negev, and the Katif Research
propagating through anisotropic turbulence,” Proc. SPIE, vol. 9614, and Development Center, Ministry of Science, Sedot
Sep. 2015, Art. no. 96140B. Negev Regional Council, Netivot, Israel. His current
[26] I. Toselli and O. Korotkova, “Optical turbulence with anisotropy research interests include the applications of electro-optical and imaging
at different scales and its effect on laser beam propagation systems in remote sensing, especially in agriculture (precision and smart
along vertical paths,” Proc. SPIE, vol. 9465, May 2015, agriculture), bio-meteorology, optical wave propagation, target acquisition,
Art. no. 94650P. and optical spectroscopy.

Authorized licensed use limited to: Ben-Gurion University of the Negev. Downloaded on August 05,2020 at 05:12:17 UTC from IEEE Xplore. Restrictions apply.