Doppler Broadening of

Acoustic Waves Scattered by

Turbulence Flowing With a

Horizontal Wind

ARL-TR-2131 July 2000

George H. Goedecke, Roy C. Wood, and Harry J. Auvermann

Approved for public release; distribution unlimited.
The findings in this report are not to be construed as an

official Department of the Army position unless so

designated by other authorized documents.

Citation of manufacturer’s or trade names does not

constitute an official endorsement or approval of the use

thereof.

Destroy this report when it is no longer needed. Do not

return it to the originator.
ARL-TR-2131 July 2000

Army Research Laboratory

Adelphi, MD 20783-1197

Doppler Broadening of

Acoustic Waves Scattered by

Turbulence Flowing With a

Horizontal Wind

George H. Goedecke

New Mexico State University

Roy C. Wood

Armstrong Atlantic State University

Harry J. Auvermann

Computational and Information Sciences Directorate, ARL

Approved for public release; distribution unlimited.
Abstract

This report investigates the scattering of a monochromatic acoustic wave

by an assumed unchanging atmospheric turbulence distribution that moves

with a uniform constant horizontal wind. The source and detector are at

rest at different positions in a ground-fixed frame. Theoretical expressions

are derived for the scattered acoustic pressure at the detector. Because distances and scattering angles change with time as any turbulent fluctuation (turbule) flows through the scattering volume, the detector signal has

a time-dependent Doppler frequency shift. A simple formula for the consequent frequency broadening is derived. A computer code is developed

that calculates the scattered signal and its fast Fourier transform scattered

by a single turbule, or by a collection of turbules of many different scale

lengths. The code uses a new approximate but quite accurate time-shift algorithm that reduces the calculation time substantially. Several numerical

results from this code are presented for a model scenario that is similar to a

recent experimental scenario, with good agreement among spectral widths

and shapes.

ii
Contents

1 Introduction 1

2 Technical Discussion 3

2.1 Scattering by One Turbule .................... 3

2.1.1 Geometry ......................... 3

2.1.2 Scattered Acoustic Pressure at Detector ........ 4

2.1.3 Turbule Model ...................... 6

2.1.4 Spectral Broadening ................... 7

2.2 Multiturbule Model ........................ 8

2.2.1 General Population .................... 8

2.2.2 Time-Shift Algorithm ................... 8

2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

3 Conclusions and Recommendations 24

References 27

Appendixes 29

A Theoretical Derivations 29

A-1 Geometry and Model . . . . . . . . . . . . . . . . . . . . . . . 29

A-2 Acoustic Equations . . . . . . . . . . . . . . . . . . . . . . . . 30

A-3 Source Wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

A-4 Scattered Waves . . . . . . . . . . . . . . . . . . . . . . . . . . 34

A-5 Temperature Scattering . . . . . . . . . . . . . . . . . . . . . . 35

A-6 Velocity Scattering . . . . . . . . . . . . . . . . . . . . . . . . 38

B Computer Algorithm 41

B-1 Detector and Source Locations and Wind . . . . . . . . . . . 41

B-2 Spectral Analysis . . . . . . . . . . . . . . . . . . . . . . . . . 42

B-3 Turbule Model Parameters . . . . . . . . . . . . . . . . . . . 44

iii
C Quick Doppler Formula 47

D Computer Codes 49

Distribution 63

Report Documentation Page 67

Figures

1 Geometry in ground-fixed reference frame ............ 3

2 Pressure signal amplitude scattered by temeprature and velocity fluctuations for one 1.5-m turbule . . . . . . . . . . . . . . . . 12

3 Frequency spectra due to temperature and velocity fluctuations

for one 1.5-m turbule . . . . . . . . . . . . . . . . . . . . . . . . . 13

4 Pressure signal amplitude scattered by temperature and velocity fluctuations for one 4.7-cm turbule . . . . . . . . . . . . . . . 14

5 Signal frequency spectra scattered by temperature and velocity

fluctuations for one 4.75-cm turbule . . . . . . . . . . . . . . . . 14

6 Pressure signal amplitude scattered by velocity fluctuations in

one tube of 0.75-m turbules, equally spaced, with spins in same

direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

7 Pressure signal amplitude scattered by velocity fluctuations for

one tube of 0.75-m turbules, randomly spaced, with spins in

same direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

8 Signal frequency spectrum scattered by velocity fluctuations for

one tube of 0.75-m turbules, equally spaced, with spins in same

direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

9 Signal frequency spectrum scattered by velocity fluctuations for

one tube of 0.75-m turbules, randomly spaced, with spins in

same direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

10 Signal frequency spectrum scattered by temperature fluctuations

for one tube of 0.19-m turbules, equally spaced, using direct

method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

11 Signal frequency spectrum scattered by temperature fluctuations

for one tube of 0.19-m turbules, randomly spaced, using direct

method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

12 Signal frequency spectrum scattered by temperature fluctuations

for one tube of 0.19-m turbules, randomly spaced, using timeshift method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

13 Signal frequency spectrum scattered by velocity fluctuations for

one tube of 1.5-m turbules, randomly spaced, with spins in same

direction and random spin directions . . . . . . . . . . . . . . . . 19

iv
14 Signal frequency spectrum scattered by temperature fluctuations

for one tube of 1.5-m turbules, randomly spaced, all with ΓT > 0

and with ±ΓT equally likely . . . . . . . . . . . . . . . . . . . . . 20

15 Total frequency spectrum due to temperature and velocity fluctuations for two full tubes with sizes 1.5 m down to 4.7 cm . . . 21

16 Total frequency spectrum due to temperature and velocity fluctuations for two full tubes . . . . . . . . . . . . . . . . . . . . . . 22

17 Total frequency spectrum due to temperature and velocity fluctuations for all four full tubes . . . . . . . . . . . . . . . . . . . . 22

18 Comparison of experimental data with model . . . . . . . . . . 23

B-1 “Primed” coordinate system . . . . . . . . . . . . . . . . . . . . . 41

B-2 Estimating minimum measurement time duration . . . . . . . . 43

B-3 Estimating maximum measurement time duration . . . . . . . . 44

D-1 Flow chart for main program . . . . . . . . . . . . . . . . . . . . 58

D-2 Flow chart for subroutine CHOICE . . . . . . . . . . . . . . . . . 59

D-3 Flow chart for subroutine MULTI . . . . . . . . . . . . . . . . . . 60

D-4 Flow chart for subroutines ROTATE, DIRECT, and TAU . . . . . 61

D-5 Flow chart for subroutines ONETRB, OUTPUT, and DOT . . . . 62

Tables

D-1 Intermediate save operations . . . . . . . . . . . . . . . . . . . . 62

v
1. Introduction

This report investigates acoustic scattering by atmospheric turbulence that

is totally frozen (stationary) in a frame comoving with a uniform constant

horizontal wind. The ambient temperature is also assumed constant and

uniform. For convenience, the turbulence is modeled as an appropriate

collection of self-similar eddies or ”turbules,” as described in past literature [1]. A turbule is a localized region of temperature and/or flow velocity

fluctuations with respect to ambient values.

An acoustic detector and a monochromatic source are at rest in a groundfixed frame at different positions. The detector and source are assumed to

be in the far-field region of a turbule at any instant but in the near-field

region of the complete stationary volume that the turbules traverse. This

is almost always the case in experimental situations. The theory developed

during this work first obtains a Born approximation expression for the time

dependence of the acoustic pressure at the detector caused by scattering

of the source pressure wave by one turbule. A computer algorithm developed during this work calculates this expression. This algorithm also calculates the fast Fourier transform (FFT) of this pressure wave and reveals that

the wave’s frequency spectrum is both shifted and broadened appreciably

around the source frequency. The broadening occurs because the scattering

angles and distances for any turbule are slowly varying functions of time

for realistic wind speeds. A simple formula for the received frequency as

a function of time is also derived by two different methods; it predicts the

same broadening as the computer algorithm.

The conventional theory [2] does not allow for this kind of time dependence, because it assumes that the detector and source are in the far-field

region of the whole scattering volume. Thus it predicts spectral broadening

only if there are time-dependent fluctuations in the comoving frame. Establishing the amount of spectral broadening and shifting for turbulence completely frozen in the comoving frame but for realistic near-field detector

and source locations thus became one of the most important objectives and

results of this work. The experimentally observed spectral broadening used

for comparison has a fairly narrow large peak around the source frequency

that is very jagged but overall appears Gaussian, with wide noisy-looking

wings at smaller amplitude. The theory developed in this report predicts

that the large narrow peak is due to large turbules, but the wide wings are

due to small ones, and also predicts that the jaggedness is due to randomness in turbule locations in the comoving frame. The computer algorithm

based on the theory yields spectra that exhibit the same general features as

some experimental data, including the spectral widths and jaggedness, for

realistic populations of turbules of many length scales.

1
Also developed in this work is an efficient approximate method for adding

contributions of many small turbules of each size that are approximately

in-line with the wind. The method is demonstrated to be accurate, and

it reduces computer running time significantly, enough so that threedimensional modeling is indeed practical.

2
2. Technical Discussion

2.1 Scattering by One Turbule

2.1.1 Geometry

The scattering geometry is defined by figure 1. The letters S, D, and T stand

for source, detector, and turbule, respectively. Coordinate z is vertical, z = 0

is the ground plane, and coordinates x and y are horizontal, ground-fixed.

The uniform constant horizontal wind is w, given by

w = w(ex cos φw + ey sin φw) , (1)

where (ex, ey, ez) are the Cartesian unit basis vectors and φw is the azimuthal angle of the wind as shown. The source S and detector D are at

locations given by position vectors (RS, RD), respectively (fixed with respect to the ground-fixed reference frame) and the turbule center T has

time-varying position vector RT (t).

The location of the “center” of the turbule T is given by

RT (t) = b ̄ + wt , (2)

so b ̄ is the location of the turbule at t = 0 in the ground-fixed frame.

The vectors RT S, RDT are defined by

RT S(t) ≡ RT (t) − RS, RDT (t) = RD − RT (t) . (3)

In the comoving (barred) frame, defined by

x ̄ = x − wxt, y ̄ = y − wyt, z ̄ = z , (4)

so that the origins coincide at t = 0, we have

R ̄ T = b ̄ , R ̄ S(t) = RS − wt , R ̄ D(t) = RD − wt (5)

so that

R ̄ T S(t) = R ̄ T − R ̄ S(t) = b ̄ − (RS − wt) = RT S(t) and (6)

R ̄ DT (t) = R ̄ D(t) − R ̄ T = RD − wt − b ̄ = RDT (t) . (7)

Figure 1. Geometry in

ground-fixed reference

frame. z

x

D

RDT(t)

w

RD

RT(t)

T

RTS(t)

RS

S

y

x

D

RDT(t)

φw

w

T

RTS(t)

RS S

RT(t)

3
2.1.2 Scattered Acoustic Pressure at Detector

We seek expressions for the acoustic pressure waves ( ̄pT

as(t), p ̄vas(t)) scattered by the turbulent (temperature, velocity) fields of the turbule, as functions of time at the detector. These expressions are derived in detail in appendix A. Here, we quote the results valid to all orders of β ≡ w/c, where c

is the ambient adiabatic sound speed, taken to be constant and uniform in

the scattering volume in this work. The source is assumed to be monochromatic and isotropic in its rest (ground-fixed) frame. Reflections from the

ground are ignored. From equation (A-71), we have for p ̄T

as(t) the form

p ̄T

as(t) = AT (t) exp(−iωst) + c.c. , (8)

where ωs is the source angular frequency and AT (t) is a (slowly varying)

complex-valued amplitude, given by equation (A-72):

AT (t) = −Pe

k2

sa

4πT0 γ2

 (R ̄0 )2

(D ̄s)3RDT  

R ̄ 0 · R DT  

T(K) exp[iks(R ̄0 + RDT )] , (9)

where a hat over a vector indicates that it is a unit vector. Here, ks = ωs/c,

Pe is the pressure of the source wave a distance a from the source in a nonmoving atmosphere, T0 is the constant uniform ambient temperature, and



T(K) is the spatial FFT of T ̄( ̄r − b ̄), the stationary temperature variation of

the turbule in the comoving frame, given by



T(K) ≡



d3ξT ̄(ξ) exp(−iK · ξ) . (10)

Also, the various vectors and distances are defined by equations (A-65) to

(A-70), as follows:

Ki = ks

	

RˆDT

i − ∂R ̄0

∂ ̄bi




, (11)

where

∂R ̄0

∂ ̄bi

= (D ̄s)

−1(δij + RˆDT

i β)(R ̄0 ) , (12)

with

R ̄0 = R ̄s

 (t) − βR ̄0 (13)

R ̄s

 = RT S (t) − βRDT (t) (14)

D ̄s = 

(1 − β2)(R ̄s)

2 − (β · R ̄ s

)

2 1

2 (15)

R ̄0 = (R ̄0R ̄0)

1

2 = γ2 


D ̄s − β · R ̄ s (16)

βi ≡ wi/c , γ ≡ (1 − β2)

− 1

2 (17)

and with (RT S(t), RDT (t)) given by equations (6) and (7). Note that all the

quantities in equations (10) to (16) and thus AT (t) are functions of time that

vary slowly compared to exp(−iωst).

4
These results are quite complicated, primarily because of retardation effects. Only when β → 0 are they simple. Then we get

R0 ≈ Rs ≈ RT S

 , D ̄s ≈ R ̄s ≈ RT S, (18)

Ki ≈ ks(RˆDT

i − RˆT Si ) , and (19)

AT = −Pe(k2

sa

/4πT0 )

cos θ

RDT RT S  

T(K) exp[iks(RT S + RDT )] , (20)

where θ is the scattering angle, given here by

cos θ = Rˆ DT · Rˆ T S . (21)

These quantities are still weakly time dependent for β <<< 1. If we define

kˆs = Rˆ T S , Apw = −Pea

(RT S)−1 exp(ikRT S) ,

ˆr ≡ Rˆ DT , r = RDT , (22)

then the amplitude AT of equation (20) has the form

AT = Apw exp(iksr)

r

 

k2

s



T(K)

4πT0



cos θ , (23)

with

cos θ = kˆs · ˆr , K = ks(ˆr − kˆs) . (24)

This is exactly the form of the Monin [3] result for temperature scattering from stationary turbulence, where Apw is the amplitude of an incident

plane wave traveling in direction kˆs.

From equation (A-81), the contribution pv

as(t) scattered by the solenoidal

turbulent velocity field of the turbule is

pv

as(t) = Av(t) exp(−iωst) + c.c. , (25)

where the complex amplitude Av(t) is defined by

Av(t) ≡



−k2

sPea

2πcγ2

 

R ̄3

0 (R ̄ 0 · Rˆ DT )(Rˆ DT ·



v (K))

(D ̄s)4 RDT 

exp[iks(R ̄0 + RDT )] , (26)

with 

v (K) the spatial FFT of the turbule velocity field v ̄( ̄r − b ̄), stationary

in the comoving (barred) frame, given by



v (K) ≡



d3ξ exp(−iK · ξ)v ̄(ξ) . (27)

For β → 0, Av(t) also reduces to the Monin result for scattering by turbulent

velocity fluctuations.

Equations (9) and (26), with the accompanying equations (10) to (17) and

(27), are used in the computer algorithm described in appendix B. To obtain 

T(K) and 

v (K) specifically, we must use a specific turbule model

described in the following section.

5
2.1.3 Turbule Model

For convenience, we use Gaussian envelopes and spherical symmetry, as in

a past paper [1]. A temperature turbule of scale length or “size” a, fixed in

the comoving (barred) frame, is written

T ̄( ̄r − b ̄)=∆Ta exp 

−( ̄r − b ̄)2

a2



, (28)

and a velocity turbule of size a is written

v ̄( ̄r − b ̄)=[Ωa × ( ̄r − b ̄)] exp 

−( ̄r − b ̄)

2

a2



. (29)

Note that Ωa is an angular velocity, and v ̄ is solenoidal.

The amplitudes (∆Ta, Ωa = |Ωa|) are scaled with a, as in a previous paper

[1]. Specifically, we designate a1 as the size of the largest turbule to be used,

the outer scale length for this problem, and choose values for

∆T1

T0

≡ ΓT ,

Ω1 a1

c

≡ Γv . (30)

Typical values are ΓT ≈ Γv ≈ 0.01. Here, Ω1 ≡ Ωa1 , ∆T1 ≡ ∆Ta1 . (We must

also choose the direction Ωˆ a.) Then, from the previous paper [1] and also

in accordance with conventional models [2], we choose for any size a

∆Ta = (a/a1)

1/3∆T1 , Ωaa = (a/a1 )1/3Ω1 a1 , (31)

or

Ωa = (a1 /a)

2/3Ω1 . (32)

From the definitions of equations (10) and (27) and from equations (28) and

(29), we get the same expressions as in appendix A,



T(K) = π3/2 a3 ∆Ta exp 

−K2a2

4



, (33)



v (K) = −i

2

π3/2a5Ωa × K exp 

−K2a2

4



. (34)

Note that there is no velocity scattering if Ωa 
 K .

Consider the lowest-order approximation for the scattering vector K, i.e.,

just like equation (19):

K(t) ≈ ks(Rˆ DT (t) − Rˆ T S(t)) . (35)

This is a good approximation if β  1. Then

K2 = 2k2

s (1 − cos θ(t)) = 4k2s sin2

θ(t)

2



, (36)

with the scattering angle θ(t) defined by equation (21). Therefore, if ksa  1,

then 

T(K) and 

v (K) are extremely small unless θ is near 0◦, so the scattering is negligible at large angles. But for ksa  1, the Gaussian envelope in

equations (33) and (34) is essentially unity for all scattering angles, so backward scattering is comparable to near forward scattering in magnitude.

6
2.1.4 Spectral Broadening

We note that both AT (t) and Av(t) contain the factor

exp{−iωs[t − (R ̄0(t) + RDT (t))/c]} ≡ exp(−iΦ(t)) , (37)

where the equality defines the phase Φ. Consider any time t = t0 . Then

write for t near t0

Φ(t) ≈ Φ(t0 )+(Φ( ̇ t0 ))(t − t0 ) + neglect , (38)

where we can neglect terms of higher order in (t − t0 ) if R ̄0(t) and RDT (t)

are “slowly varying,” which is the case for β  1. Thus Φ( ̇ t0 ) is the effective

angular frequency ω(t0 ) at time t0 . Since t0 is any time, we have

ω(t) = ωs[1 − 1

c

(  ̄ ̇

R0 (t) + R ̇ DT (t0 ))] . (39)

When equations (12), (3), and (13) to (17) are used, obtaining an expression

for ω(t) valid to first order in β is very simple. Since c−1 times the time

derivative of any position vector or distance is proportional to β, then for

equations (13) to (17), we may put, to order β,

 ̄ ̇

R0 (t) ≈ R ̇ T S(t) = dt(RT S · RT S)

1/2 = Rˆ T S · R ̇ T S = w · Rˆ T S (40)

R ̇ DT (t) = dt(RDT · RDT )

1/2 = Rˆ DT · R ̇ DT = −w · Rˆ DT . (41)

Therefore, to order β,

ω(t) ≈ ωs[1 + β · (Rˆ DT (t) − Rˆ T S(t))] . (42)

Note that this can be written in terms of the approximate scattering vector

K(t) of equation (35) as

ω(t) = ωs + w · K(t) . (43)

This relation is actually correct for the exact K(t) given by equation (11).

From equations (32) or (33), (36), and (33) and (34), and from figures 1(a)

and 1(b), it is evident that small turbules may produce much wider spectral

broadening than large ones. This is because the Gaussian envelope cuts off

the scattering from large turbules when |w · K| is “large,” but does not

do so for small turbules. Thus from small turbules, we can get appreciable

backscattering, so the relative frequency shift w · Kωs may vary from zero

to ±2β, while for large turbules, |w · K|ωs  2β for any scattering angles

that produce appreciable scattering.

It should be emphasized that this general behavior is not peculiar to a turbule model, or to Gaussian model turbules, but is common for any model

of turbulence that contains both large and small length scales and is frozen

in a comoving frame, provided only that the source and detector are in the

near-field region of the complete scattering volume Vs. This applies to the

scenarios of interest in this work. If (and only if) the source and detector

are both in the far-field region of Vs at all times, then the conventional result ensues that moving frozen turbulence may produce only a frequency

shift but no broadening [2].

In appendix C, a different simple derivation of equation (42) is provided.

7
2.2 Multiturbule Model

This section contains a summary description of the scheme used later in

the report for populating the scattering volume with turbules. Rather than

having a continuous distribution of sizes for the turbule size parameter a, a

discrete set of sizes an will be selected. For each size class index n, a cubical

cell containing one turbule is defined whose linear dimension dn is an appropriate multiple of an. Suppose that the time interval for which data are

desired is T and that the average wind velocity is w. Then an approximate

scattering volume length L is equal to wT. A tube is defined to be a region of space consisting of a sufficient number of cells Mn aligned along the

wind direction so that the total length of the tube is L or greater. The several

tubes for each size class n are stacked laterally and vertically to overspread

the scattering volume. For this discussion, consider that an+1/an = 1/2.

Then dn+1/dn = 1/2 and Mn+1/Mn = 2. The number of tubes for class

n + 1 is four times the number for class n.

2.2.1 General Population

Our interest is in the scattered acoustic pressure produced by a multiturbule model that contains appropriate numbers of turbules of many sizes,

from the largest (outer scale) a1 to the smallest (inner scale) aNs , with Ns

different scales included. As shown by Goedecke et al [1], if we define the

packing fraction

φ = naa3 , (44)

where na is the number density of turbules of size a, then φ must be independent of a so that the energy transfer rate per unit mass in a Kolmogorov

cascade can be scale-independent as well. On average, one turbule of size

a is somewhere in each cubical cell of size da, where

da = n−1/3 a = a/φ1/3 . (45)

For a typical φ ≈ 10−3 , da = 10a .

If we had unlimited computer resources, we could populate a defined scattering volume Vs with turbules of different sizes (one “randomly” located

in each cell for each size) and simply add up the (AT (t), Av(t)) from each.

This is actually what we do for the larger turbules with ksa  1. But if there

are perhaps 100 turbules with a = 2 m in Vs, then for a = 0.02 m = 2 cm, according to equation (44), there will be (100)(2/0.02)3 = 108 turbules, a truly

formidable number.

2.2.2 Time-Shift Algorithm

Fortunately, there seems to be a good approximate way to handle the multitude of smaller turbules that reduces the calculation time significantly.

Imagine two identical turbules, 1 and 2, more or less in-line along the wind

direction, chosen to be the x-direction. That is, let the centers of their two

8
cells have coordinates ( ̄b

(1)

cx ,  ̄b

(1)

cy ,  ̄b

(1)

cz ) and ( ̄b

(2)

cx ,  ̄b

(1)

cy ,  ̄b

(1)

cz ); i.e., their cells are

in-line with the wind.

Now, consider equation (9) for AT

(1)(t), from turbule 1, and assume that the

source and detector are always in the radiation zone of a comoving cell of

side d = a/φ1/3 in which the turbule is located. This implies RDT

(1) (t) and

RT S

(1) (t)  (d, λ, d2/λ), whichever is greater, at all (or most) times. For such

cases, vectors such as (RT S(1),...) can be replaced by (RT S(1) c ,...), the vectors to the center of the cell, except in the factor exp[iks(R ̄(1)

0 + RDT(1))]. In

this factor, we may write

exp[iks(R ̄(1)

0 + RDT(1))] ≈ exp[iks(R ̄(1)

0c + RDT(1) c )] exp(−iK(1)

c · η1) , (46)

where η1 is the displacement of the center of turbule 1 from its cell center.

So we get

AT

(1)(t) = AT(1)(t) exp(−iK(1)

c (t) · η1) , (47)

where from equations (9) and (46), we have

AT

(1)(t) ≡ −Pe(k2

s a/4πT00 γ2)[(R ̄(1)

0c )

2/(D ̄s(1)

c )

3RDT(1) c ]

×(R ̄ (1)

0c · Rˆ DT(1) c )T

0 (K(1)

c ) exp[iks(R(1)

0c + RDT(1) c )]. (48)

That is, AT

(1)(t) is the old AT(1)(t) but calculated with vectors to the cell center

times the phase factor exp[iks(R(1)

0c + RDT(1) c )].

Now, for turbule 2, clearly (AT

(2)(t), K(2)

c (t)) are simply time-shifted from

(AT

(1)(t), K(1)

c (t)). That is, we have

AT

(2)(t) = AT(1)(t + τ2) , K(2)

c (t) = K(1)c (t + τ2) , (49)

where τ2 ≡ w−1( ̄b

(2)

cx −  ̄b

(1)

cx ). For example, suppose turbule 2 is one cell

ahead of turbule 1, so that b

(2)

cx − b

(1)

cx = d. Then τ2 = d/w, so AT

(2)(t) is just

like AT

(1)(t), but occurs earlier.

Thus, the time-shift algorithm for a whole tube of M turbules of a given

size, with cell centers equally spaced by d along the wind direction, is

obviously

AT

tube(t) = 

M

=1

AT

(1)(t + τ) exp[−iK(1)

c (t + τ) · η] , (50)

where

τ = w−1( ̄b()

cx −  ̄b(1)cx )=() − 1)(d/w) , (51)

and η is chosen “randomly” for each ); i.e., −d/4 ≤ (ηx, ηy, ηz) ≤ d/4,

with uniform independent probabilities. So, for a given tube, we need to

9
compute the complicated functions (AT

(1)(t), K(1)

c (t)) only once; then, for

each turbule ) in the tube, we choose an η and compute

AT

()(t − τ) = AT(1)(t) exp(−K(1)

c (t) · η) , (52)

then simply time-shift this to get AT

()(t), and add that to the accumulated

AT

tube(t).

By comparing with direct computation of the AT

tube(t) for the same tube

and using the full equations (9) to (16) for each turbule with its center

at the same set of random locations as η from above, we found that for

νs = 500 Hz, the time-shift algorithm is a good approximation for d 

4 m. However, the algorithm is not as good for d ≥ 8 m, with RDT

c and

RT S

c  10 m minimum but actually  100 m during the times of largest signal production. Since a direct computation involves ∼100 operations per

turbule and a time-shift of only ∼5 operations, the computing time advantage of the time-shift algorithm is major.

These same considerations apply for the Av(t) of equation (26), except that

we need a separate Ωˆ for each turbule. For any turbule, we may write in

general

Av(t) = Ωˆ · Bv

(t) , (53)

where from equations (26) and (34), we have

Bv(t) = −



iπ1/2a5ΩPea

4cγ2

 R ̄3

0 (R ̄ 0 · Rˆ DT )

(D ̄s)4RDT 

exp(−K2a2/4) exp[iks(R ̄0 + RDT )](Rˆ DT × K) .

(54)

The time-shift algorithm for a tube of M turbules along the wind is now

Av

tube(t) = m

=1

Ωˆ  · Bv(1)

(c) (t + τ) exp[−iK(1)

c (t + τ) · η] , (55)

where (Bv(1)

(c) , K(1)

c (t)) are calculated with vectors to the center of the cell

containing and comoving with turbule 1. We calculate the same way as for

AT

tube(t), except that we must choose an Ωˆ  as well as an η for each turbule.

2.3 Results

Many of the theoretical results are presented and discussed in the previous

section and in appendixes A to C, so we limit this section primarily to presenting and discussing graphs that illustrate some of the most interesting

predictions of the computer code (Slow Eddy) that was developed on the

basis of the theoretical results. However, we will also show how the theory

easily predicts the major features of these computational results.

10
For all the following results, the detector was placed at (−150, 0, 0) m and

the source at (150, 0, 10) m, and the wind was chosen with components

(3.1, 0, 0) m/s. Thus, the primed and unprimed ground-fixed frames described in section 2.1.1 and discussed in appendix B were the same, and

the wind velocity was blowing from the detector toward the source, in the

(+x) direction. The total time T was 24 s, 0 ≤ t<T, tube length was

L = wT = 756.4 m centered at the origin, the time steps ∆t were 0.0403 s,

and the total number of points used in the calculation was N = 6048. (The

point n = 6049 ↔ t = T is “equivalent” to t = 0 because of the inherent periodicity of the discrete FFT.) The adiabatic sound speed was taken

to be c = 340 m/s, corresponding to an ambient temperature T0 ≈ 15 ◦C

≈ 288 K ≈ 60 ◦F. For multiturbule computations, each turbule angular

velocity unit vector was generally chosen to have uniform random distribution in 4π solid angle, independently of other turbules. This orientation

distribution produces a model of isotropic turbulence, as was discussed in

another paper [1]. The ∆Ta in the scattering volume Vs were generally chosen to have the same sign (positive) because temperature turbules are air

pockets warmer than ambient that have risen from the ground because of

buoyancy forces, such as those over a paved runway that was used in the

relevant experiment. Over a larger volume, the mean ∆T probably should

be zero. In the following results, one calculation was made with each turbule ∆T either > 0 or < 0 with equal probability to contrast the results of

the different conditions. Also, some calculations were made with the angular velocity vectors of the larger turbules all in the +y-direction, for the

same reason.

The ratio of cell size to turbule size was chosen to be d/a = 8, yielding

a packing fraction φ = 1.95 × 10−3. The number of sizes in the cascade

was chosen to be Ns = 6, and the (largest, smallest) turbule sizes (a1 , aNs )

were chosen as (1.5 m, 0.046825 m), respectively. Thus the parameter Nd

of equation (B-18) was equal to 2 (whereby μ = ln 2), such that each successively smaller size was half the previous size. The relative fluctuation

strengths (∆T1/T0 ≡ ΓT , Ω1 a1 /c ≡ Γv) of equation (30) were chosen as

(ΓT = 0.013, Γv = 0.027). From our previous paper [1] or equation (B-17)

of this report, these values yield C2

v = 0.088 m4/3/s2, C2

T = 0.044 K2/m2/3.

These values are in good agreement with those given by Brown and Clifford [4]. Also, this value of C2

v was inferred∗ from wind data that were

measured, along with the time-dependent acoustic pressure received by

many different detectors, during an experiment by Galindo and Havelock

[5]. All the numerical values chosen in this report are supposed to provide

a good approximation to those that existed in that experiment. (As mentioned earlier, a main objective of this work is to determine whether a realistic model of acoustic scattering by turbulence that is frozen in a frame and

comoving with a wind can possibly simulate actual experimental results.)

Using these values of the input parameters and the method and computer

code described in this subsection and in appendix B, we calculated total

∗D. K. Wilson of ARL computed this value.

11
time-dependent complex acoustic pressure amplitudes (AT (t), Av(t)) scattered by individual turbules, tubes of turbules, and all tubes of turbules in

several different scattering volumes and then Fourier transformed to get

the spectra (A ̃T (ω), A ̃v(ω)). For convenience, each (AT , Av) of equations (9)

and (26) was “normalized” by division by the factor Pe with a = 1 m. Thus,

Pe is now the pressure amplitude of the source wave at a distance a = 1 m

from the source in a still atmosphere. In what follows, the (AT , Av) are the

normalized ones, which are dimensionless.

In figure 2, we plot the real and imaginary parts and the magnitudes of

(AT (t), Av(t)), respectively, for a single turbule with a = a1 = 1.5 m,

with  ̄bx = −L/2,  ̄by = 0,  ̄bz = 20 m. Here,  ̄bz = 20 m corresponds to

an above-ground base height zb = 14 m, plus d1 /2 = 6 m, half the side

of a cubical cell having d1 = 8a1 . In figure 3, we plot the magnitudes of

the FFTs A ̃T (ω) and A ̃v(ω). One can easily understand from equation (33)

why the pulse in figure 2 has a fairly narrow width in time: The Gaussian envelope rapidly kills the signal as the scattering angle increases, since

ksa1 ≈ 13.9. The frequency width can be obtained easily from equation (42),

ω(t) ≈ ωs[1+β·(Rˆ DT (t)−Rˆ T S(t))] . That is, calculate or obtain from figure

Figure 2. Pressure signal

amplitude scattered by

(a) temperature and

(b) velocity fluctuations

for one 1.5-m turbule.

0.00010

0.00005

0

–0.00005

0.00010

Im{AT}

Re{AT}

|AT|

Time (s)

0 40 80 120 160 200 240

Relative pressure

(a)

0.00010

0.00005

0

–0.00005

–0.00010

Im{Av}

Re{Av}

|Av|

Time (s)

0 40 80 120 160 200 240

Relative pressure

(b)

12
2 the two times (t1 , t2 ) at which the time-dependent magnitudes | A(t) | are

perhaps 5 percent oftheir maximum values and find Rˆ DT and Rˆ T S at these

times. Then verify from figure 3 that | ω(t1) − ω(t2) | is a good estimate of

the (angular) frequency width of the A ̃(ω). This calculation works because

the complete frequency range νs(1 − 2β) to νs(1 + 2β) is “scanned” by any

turbule, but only a small part of it is manifested in the scattered waves from

large turbules.

In figures 4 and 5, we repeat the presentations of figures 2 and 3, respectively, but for a single small turbule of size a6 = 4.6875 cm, with the same

( ̄bx,  ̄by,  ̄bz). For this size turbule, (ksa6 )2 = 0.188, so the Gaussian envelope

of the A(t) is almost unity for all scattering angles. Therefore, each A(t) is

nonzero over the total time duration, and the frequency width extends to

almost the maximum possible, ±2βνs, obtainable by backscattering from

equation (42). Note also that each A ̃(ω) goes to zero at frequencies ν+ =

4.3 Hz, ν− = −4.2 Hz. Again, this is easily understood from equations (9),

(26), and (42): The time-dependent amplitude A(t) is zero approximately

when Rˆ DT · Rˆ T S = cos θ = 0 (to lowest order in β). From figure 3, obtain the times t3 , t4 at which this is true, and then use equation (42) to get

ω(t3) and ω(t4 ); these agree quite well (within 5–10%) with the frequencies

ν±. (The algebraic results are very sensitive to (t3 , t4 ) values, and the computer algorithm contains higher order terms in β that slightly influence the

results.) Also note the asymmetry of the A(t) and A ̃(ω); this is due to the

source height being different from the detector height, in this geometry.

Figures 6 to 9 contrast the effects of randomness versus regularity in turbule

center locations for many turbules. In figure 6, we present |Av(t)| for a tube

of turbules with a = 0.75 m, with all angular velocities in the y-direction

and each turbule center at the center of its comoving cell of side length

da = 8a = 6 m. This tube forms a moving “diffraction grating.” The tube

centerline is ( ̄by = 0,  ̄bz = 20 m), and there are Ma = L/da = 126 turbules

contributing at any time. (In all calculations done for tubes, each cell in a

Figure 3. Frequency

spectra due to

temperature and

velocity fluctuations for

one 1.5-m turbule.

0.05

0.04

0.03

0.02

0.01

0

|AT|

|Av|

~

~

Frequency (Hz)

–0.5 –0.4 –0.3 –0.2 –0.1 0 0.1 0.2 0.3 0.4 0.5

Relative pressure spectral magnitude

13
Figure 4. Pressure signal

amplitude scattered by

(a) temperature and

(b) velocity fluctuations

for one 4.7-cm turbule.

(b)

3×10–9

2×10–9

1×10–9

0

–1×10–9

–2×10–9

–3×10–9

0 40 80 120 160 200 240

Time (s)

|AT|

Re{AT}

Im{AT}

Relative pressure

(a)

3×10–9

2×10–9

1×10–9

0

–1×10–9

–2×10–9

–3×10–9

0 40 80 120 160 200 240

Time (s)

|Av| Re{Av}

Im{Av}

Relative pressure

Figure 5. Signal

frequency spectra

scattered by

temperature and

velocity fluctuations for

one 4.75-cm turbule.

5×10–7

4×20–7

3×10–7

2×10–7

1×10–7

0

|Av| ~ |AT| ~

Frequency (Hz)

–10 –5 0 5 10

Relative pressure spectral magnitude

14
fixed tube of total length L is kept filled; i.e., as a turbule is leaving at L/2,

another is entering at −L/2. This ensures steady-state conditions.) The time

dependence is extremely regular, as one would expect.

Figure 7 is for the same conditions as figure 6, except that each turbule

center is randomly placed inside its comoving cell, with no correlations

among placements of different turbules. The randomness consists of uniform uncorrelated probabilities of a turbule center’s (x, y, z) locations being

between −d/4 and d/4 from its comoving cell center. It is interesting (but

totally expected) that this randomness produces extremely large relative

fluctuations.

Figure 6. Pressure signal

amplitude scattered by

velocity fluctuations in

one tube of 0.75-m

turbules, equally

spaced, with spins in

same direction.

0 40 80 120 160 200 240

Time (s)

0.00025

0

Relative pressure magnitude

Figure 7. Pressure signal

amplitude scattered by

velocity fluctuations for

one tube of 0.75-m

turbules, randomly

spaced, with spins in

same direction.

0.00025

0

0 40 80 120 160 200 240

Time (s)

Relative pressure magnitude

15
In figures 8 and 9, we plot the magnitudes of the FFTs for the conditions of

figures 6 and 7, respectively, on a semilog plot for good visibility. In figure

8, the principal maxima are striking; they are spaced just as they should

be, i.e., at frequencies νq = qMa/T ≈ (q)(0.516) Hz, where q is an integer. There should be exactly zero contribution for other frequencies, for

the discrete FFT. The principal maxima diminish rapidly as |ν| increases,

probably because of the Gaussian envelope of a single turbule. (The picture is analogous to a multiple-slit interference pattern as modified by the

single-slit diffraction pattern, in optics.) In figure 9, the randomness in turbule placement yields a jagged continuous spectrum. Both figures 8 and

9 display a remarkable straight-line behavior, which corresponds on these

semilog plots to an exponential spectral envelope, exp(−|ν|τa), where τa is

a constant with the dimension of time that most likely depends on the scale

length a. We have not yet been able to predict this behavior theoretically.

Figure 8. Signal

frequency spectrum

scattered by velocity

fluctuations for one

tube of 0.75-m turbules,

equally spaced, with

spins in same direction.

10

0.1

0.001

0.00001–2 –10 12

Frequency (Hz)

Relative pressure spectral magnitude

Figure 9. Signal

frequency spectra

scattered by velocity

fluctuations for one

tube of 0.75-m turbules,

randomly spaced, with

spins in same direction.

1

0.1

0.001

0.0001

0.00001

Frequency (Hz)

–2 –10 1 2

Relative pressure spectral magnitude

16
Figures 10 to 12 are intended to illustrate that the time-shift algorithm gives

results essentially the same as direct calculation. These figures are |A ̃T (ω)|

for a tube of Ma = 504 turbules with a = 1.5/23 = 0.1875 m, with the

centerline as previously mentioned, with all ∆Ta > 0. Figure 10 is for regular placement (for reference), with the use of “either” the direct or the

time-shift calculation; i.e., the results are identical for regular placement

of the turbule centers at their cell centers initially. Figure 11 is for random

placements using the direct calculation. Figure 12 is for the “same” random

placements but using the time-shift calculation, which is eight times faster

than the direct method. Clearly the two methods produce substantially the

same results, although there are small differences, as there should be, because the time-shift method is not exact.

Figure 10. Signal

frequency spectrum

scattered by

temperature

fluctuations for one

tube of 0.19-m turbules,

equally spaced, using

direct method.

10–1

10–2

10–3

10–4

10–5

10–6

Frequency (Hz)

–10 –5 0 5 10

Relative pressure spectral magnitude

Figure 11. Signal

frequency spectrum

scattered by

temperature

fluctuations for one

tube of 0.19-m turbules,

randomly spaced, using

direct method.

10–1

10–2

10–3

10–4

10–5

10–6

Frequency (Hz)

–10 –5 0 5 10

Relative pressure spectral magnitude

17
Figure 12. Signal

frequency spectrum

scattered by

temperature

fluctuations for one

tube of 0.19-m turbules,

randomly spaced, using

time-shift method.

10–1

10–2

10–3

10–4

10–5

10–6

–10 –5 0 5 10

Frequency (Hz)

Relative pressure spectral magnitude

Figures 13 and 14 illustrate the effects of random versus nonrandom angular velocity directions Ωˆ and temperature amplitude signs ∆T /|∆T|. They

are for a tube of Ma = 63 turbules with a = 1.5 m, with the centerline at

by = 6 m, bz = 20 m. Figure 13 plots |A ̃v(ω)|; figure 13(a) is for all the Ωˆ

in the y-direction, while figure 13(b) is for the Ωˆ chosen isotropically at

random. Figure 14 plots |A ̃T (ω)|; figure 14(a) is for all ∆T > 0, while figure 14(b) is for sgn(∆T) = ±1 with equal probability (this would make

T0 equal to the mean temperature in Vs). In each case, the turbule center

locations were randomized with the same set of random numbers, so differences in the comparable graphs must be caused by the different choices

of Ωˆ or sgn(∆T). What is remarkable is that the differences are not large.

This can be understood by noting that, for example, changing the sign of

∆T /|∆T| for a given turbule introduces a phase change of π in its scattered wave, which could have been done equally well by moving its center

about a half wavelength. So, in general, randomization of Ωˆ and ∆T /|∆T|

is nearly equivalent to nonrandomization with a little more “freedom of

movement” of the turbule centers.

18
Figure 13. Signal

frequency spectrum

scattered by velocity

fluctuations for one

tube of 1.5-m turbules,

randomly spaced, with

(a) spins in same

direction and

(b) random spin

directions.

0.8

0.6

0.4

0.2

0

Frequency (Hz)

–0.2 –0.1 0 0.1 0.2

Relative pressure spectral magnitude

(a)

0.8

0.6

0.4

0.2

0

Frequency (Hz)

–0.2 –0.1 0 0.1 0.2

Relative pressure spectral magnitude

(b)

19
Figure 14. Signal

frequency spectrum

scattered by

temperature

fluctuations for one

tube of 1.5-m turbules,

randomly spaced, (a) all

with ΓT > 0 and (b)

with ±ΓT equally likely.

0.16

0.12

0.08

0.04

0

Frequency (Hz)

–0.2 –0.1 0 0.1 0.2

Relative pressure spectral magnitude

(a)

0.16

0.12

0.08

0.04

0

Frequency (Hz)

–0.2 –0.1 0 0.1 0.2

Relative pressure spectral magnitude

(b)

20
Figures 15 to 17 are sums over contributions of many turbules of all sizes

used, i.e., aα = (1.5, 0.75, 0.375, 0.1875, 0.09375, 0.046875) m for α = (1,

..., Ns = 6), respectively. The amplitudes (Ωa, |∆Ta|) were scaled as discussed earlier, so these calculations are an attempt to numerically represent

realistic, isotropic, homogeneous, fully developed turbulence by a turbule

model, following the formulas by Goedecke and Auvermann [1]. This attempt is a compromise with the computing power of a PC with a Pentium II

233 MHz processor. We really should have included both larger and smaller

turbules, plus many more different sizes in between. But what we did already required ∼4.25 hr per large tube (d1 = 8a1 = 12 m), corresponding

to 30 s per turbule, and this was only after we compromised further by

making a reasonable approximation for the smallest size a6 = 0.046875 m.

We randomly put only “one” such turbule instead of eight in each cell of

the next “larger” size, d5 = 8a5 = 0.75 m, and then we multiplied the result by eight. We checked this for a single tube of side length d5 versus the

“correct” method using turbules of scale a6 only, and the results were quite

close. Since so many such tubes are in one big tube ((d1/d5)2 = 256), the

total error induced by this economy should be quite small.

Figure 15 plots |A ̃T (ω) + A ̃v(ω)| for two big tubes with centerlines by =

± 6 m, bz = 20 m. Figure 16 does the same for two big tubes above those

of figure 15: i.e., by = ± 6 m, bz = 32 m. Figure 17 is the total frequency

spectrum from these four tubes. Note that the contributions from the higher

big tubes are significantly smaller, and that the general shapes of the plots

are essentially unchanged, so these shapes are established already by just

one or two tubes. The semilog plots clearly show the contributions of the

small turbules all the way to ±2βvs = ± 10 Hz.

0.6

0.5

0.4

0.3

0.2

0.1

0

Frequency (Hz)

–10 –5 0 5 10

1

0.1

0.01

0.001

0.0001

0.00001

Frequency (Hz)

–2 –10 1 2

(a) (b)

|A|

~

|A| ~

Figure 15. (a) Total frequency spectrum due to temperature and velocity fluctuations for two full tubes with sizes 1.5 m down to 4.7 cm. (b) Log plot of (a) from

±10 Hz.

21
0.6

0.5

0.4

0.3

0.2

0.1

0

Frequency (Hz) Frequency (Hz)

–2 –10 1 2

1

0.1

0.01

0.001

0.0001

0.00001

–10 –5 0 5 10

(a) (b)

|A|

~

|A|

~

Figure 16. (a) Total frequency spectrum due to temperature and velocity fluctuations for two full tubes located above the two tubes in figure 14, with sizes 1.5 m

down to 4.7 cm. (b) Log plot of (a) from ±10 Hz.

0.6

0.5

0.4

0.3

0.2

0.1

0 –10 –5 0 5 10

1

0.1

0.01

0.001

0.0001

0.00001

Frequency (Hz) Frequency (Hz)

–2 –10 1 2

(a) (b)

|A|

~

|A|

~

Figure 17. (a) Total frequency spectrum due to temperature and velocity fluctuations for all four full tubes from figure 14 and 15. (b) Log plot of (a) from ±10 Hz.

22
Finally, figure 18 contrasts some experimental data with our model results.

Figure 18(a) is from a presentation by Auvermann and Goedecke [6]; the

jagged curve is the data. The central spike is about 101.56 Hz; this corresponds to ν = 0 in the previous spectral plots and in figure 18(b). The

ratio of the central spike height to the “mean” height at ν = 102 − 101.56

= 0.44 Hz is about 6/0.23 = 26. These data were for a source detector and

wind velocity scenario very similar (but not identical) to what we modeled.

Figure 18(b) is an expansion of our figure 17(a) along the abscissa to match

the total frequency width in figure 18(a). The model predictions are similar to the experimental results; e.g, the ratio of the central spike height (at

ν = 0) to the “mean” height at ν = 0.44 Hz is ≈0.52/0.02 = 26, in figure

18(b). But the jaggedness in the model results is not as great as in the experimental results. In later investigations, it has been found that the jaggedness, defined as the mean frequency spacing of adjacent peaks and valleys

in the spectrum, is approximately equal to 3/T, where T is the duration of

the time signal. The model jaggedness agrees with the experimental jaggedness, if T is the same for both.

Figure 18. Comparison

of experimental data

with model: (a) Frequency spectrum from

experiment and (b) frequency spectrum from

figure 17(a).

Frequency (Hz)

Relative sound pressure (FFT)

Measured signal

Turbule, rad = 1.8040 m

Turbule, rad = 1.4937 m

Turbule, rad = 1.2368 m

Data set GA3006, 500 Hz

× 104

7

6

5

4

3

2

1

0

101 101.1 101.2 101.3 101.4 101.5 101.6 101.7 101.8 101.9 102

(a)

0.6

0.5

0.4

0.3

0.2

0.1

0

Frequency (Hz)

–0.56 –0.36 –0.16 0.04 0.24 0.44

(b)

|A|

~

23
3. Conclusions and Recommendations

The theory developed in this work allows understanding of the spectral

broadening and shifting in an acoustic wave scattered from moving turbulence on the basis of two simple general physical principles. The first

principle is that, even if the turbulence is frozen in a frame comoving with

a wind, spectral broadening occurs because the scattering angles and the

distances from source to turbulent eddies and from the eddies to the detector are time-dependent. The analog in everyday experience is what you

hear when an approaching car’s horn is blowing at frequency νs and then

recedes while you are standing to one side of the car’s path: You actually

hear the whole Doppler-shifted frequency range, from νs(1+β) to νs(1−β),

as a function of time. (For scattered waves, the factors are 1±2β, where β =

(speed of object/speed of sound), to first order in β.)

The second principle is that large objects (compared to the wavelength)

scatter strongly but mainly in the near-forward direction, while small objects scatter weakly but nearly isotropically. Thus from large eddies, we get

only a small spectral width, while from small eddies, we get ±2βνs but less

amplitude.

Of course, in actuality the turbulent flow in a frame comoving with the

mean wind is also time-dependent. But we have shown via these general

principles and a computer model that this time dependence is unnecessary to explain the general features of some observed acoustic data. This

is markedly in contrast to the conventional treatment of spectral broadening, which attributes all of spectral broadening to time dependence in the

comoving frame. In that treatment, the source and detector are assumed to

be in the far field of the whole turbulent scattering volume, so that, effectively, scattering angles and distances are not time-dependent. This is an

unrealistic assumption for acoustics, where both the source and detector

are invariably in the near field of the scattering volume.

One of the approximations apparently made in the theory developed in this

report was the assumption of a “uniform” mean wind and a “uniform” ambient temperature. (The commonly used ambient temperature is taken as

the mean temperature in the scattering volume Vs; we did not so restrict our

ambient temperature.) But these may only be “apparent” approximations.

For example, suppose we define the mean wind as the time and space averaged wind over the time of an experiment and over the whole scattering

volume from the ground up. Then of course it is by definition uniform and

constant, so “all” wind variations are part of the turbulent flow. This point

of view would introduce anisotropy and/or inhomogeneity in the turbulent flow, which is easy to do in a turbule model. Also, suppose we define

the mean temperature as the same time and space average of the nonacoustic temperature. Then all nonacoustic variations of temperature in Vs

24
are part of the turbulence. For example, suppose we used thermometers to

“measure” a “mean temperature,” and found these readings to be a function of height only, in Vs; typically, they would decrease with height. We

could model this in terms of a constant ambient T0 , the average over height

of the thermometer readings, plus turbules whose amplitudes ∆T did not

cancel on average and depended on height in such a way as to reproduce

the thermometric results. This method would introduce inhomogeneity in

the turbulent temperature variations, and again this is easy to do with a

turbule model. These ideas may be worth further investigation. However,

this approach would not allow for the curved limiting ray paths that create

shadow zones in the source radiation patterns. This curvature is small and

was neglected in this work, because we treated scattered radiation only.

Another approximation that was made in the theory and thus in the developed code was that the source and detector are instantaneously in the

radiation zone of any turbule, i.e., that RT S and RDT are large compared to

the larger of (λ, a2/λ), at all times. This is not true at all times for all turbules

in realistic scenarios, in which a few of the turbules may actually intercept

the source and/or detector and in which a substantial fraction of the large

turbules spend some time close to the source and detector. The radiation

zone expressions were used as a compromise, because they yield the same

results as the conventional theory in the limit of β → 0. We expect that corrections to this approximation would not alter the general character of the

total |A ̃(ω)|, but this expectation should be checked in the future. The fully

correct theoretical expressions for the Born approximation scattered waves

in near-field regions of a turbule can be obtained for the general model

used, if we do not make the approximations mentioned just before equations (A-21), (A-37), and (A-52) and just after equations (A-38) and (A-39).

But then, the expressions for the waves are extremely complicated and difficult to work with, so much so that it seems unlikely a three-dimensional

model could be handled.

The computer code Slow Eddy developed in this work is not optimized.

It can model a three-dimensional problem only roughly on a PC and only

because the new time-shift algorithm is so much faster than direct computation. One improvement that might still be made in the code is to allow the

scattering volume to be different for each turbule size: The larger the turbules, the smaller the volume needed. But a recipe for choosing the optimal

scattering volumes is not so easy to attain, because the choices depend on

virtually all the parameters of a scenario. It is probably best to experiment

with this code to help determine how to modify it to use scale-dependent

scattering volumes. Unfortunately, the effective volume needed increases

inversely with turbule size, so this would not help much in reducing running time.

Other facets of an experimental scenario are that source waves reflected

from the ground will be incident on the turbulence and that the source

and/or the detector may be anisotropic. The latter could be modeled by

suitable modifications of the computer code, along with some further

25
theoretical development. The former would require substantially more theoretical development and code alterations. The general character of scattered wave spectra may be altered appreciably by both effects.

Overall, the frequency spectra predicted by this model of time-dependent

turbulence seem to have the same general characteristics as some observed

spectra. How well this model applies to many different sets of experimental

results needs to be tested. If it seems to work fairly well, then quite likely, it

could be used to analyze acoustical remote sensing data with a view toward

inferring detailed information about atmospheric turbulence and/or wind

speeds.

26
References

1. G. H. Goedecke and H. J. Auvermann, “Acoustic scattering by atmospheric turbules,” J. Acoust. Soc. Am. 102 (1997), pp 759–771.

2. V. I. Tatarskii, The Effects of the Turbulent Atmosphere on Wave Propagation,

Keter, Jerusalem (1971).

3. A. S. Monin, “Characteristics of the scattering of sound in a turbulent

atmosphere,” Sov. Phys. Acoust. 7, No. 130 (1962).

4. E. H. Brown and S. F. Clifford, “On the attenuation of sound by turbulence,” J. Acoust. Soc. Am. 60 (1976), pp 788–793.

5. M. Galindo and D. I. Havelock, “Temporal coherence of a sound field in

the turbulent atmosphere near the ground,” Proc. of the 7th International

Symposium on Long-Range Sound Propagation, Ecully, France (July 1966),

pp 149–159.

6. H. J. Auvermann and G. H. Goedecke, “Influence of the Doppler effect

on the bandwidth of acoustical signals scattered from atmospheric turbulence,” Proc. of the 8th International Symposium on Long-Range Sound

Propagation, Pennsylvania State University (9–11 September 1998).

27
Appendix A. Theoretical Derivations

A-1 Geometry and Model

See figure 1 of the main report. In the derivation that follows—

t = clock time

x, y = horizontal coordinates, z = vertical coordinate

r = position vector

w = ambient wind velocity, uniform, constant

T0 = ambient temperature

p0 = ambient pressure

ρ0 = ambient mass density

RS ≡ s = source location

RD ≡ d = detector location

RT (t) = RT+ wt ≡ b ̄ + wt = location of turbule “center”

RT S ≡ RT (t) − RS

RDT ≡ RD − RT (t)

The DS or ground-fixed frame is the “unbarred ” frame (D, S = detector,

source). The T (barred) frame is comoving with the wind (T = turbule).

Their Cartesian coordinates are related by

x ̄i = xi − wit, t

 ̄= t Galilean transformation. (A-1)

The total (temperature, pressure, mass density, and fluid velocity) fields in

the DS system are

T(r, t) = T0 + Tt(r, t) + Ta(r, t) (A-2)

p(r, t) = p0 + pt(r, t) + pa(r, t) (A-3)

ρ(r, t) = ρ0 + ρt(r, t) + ρa(r, t) (A-4)

v(r, t) = w + vt(r, t) + va(r, t) , (A-5)

where the subscripts (t, a) ≡ (turbulent, acoustic). In the T system, these

fields are

T( ̄r,t

 ̄) = T0 + Tt( ̄r, t) + Ta( ̄r, t) (A-6)

p ̄( ̄r,t

 ̄) = p0 +  ̄pt( ̄r, t)+ ̄pa( ̄r, t) (A-7)

v ̄( ̄r,t

 ̄) = v ̄t( ̄r, t) + v ̄a( ̄r, t) (A-8)

ρ ̄( ̄r,t

 ̄) = ρ0 +  ̄ρt( ̄r, t)+ ̄ρa( ̄r, t) , (A-9)

29
where

T( ̄r,t

 ̄) = Tt(r, t),... , v ̄t( ̄r, t) = vt(r, t)|r= ̄r−wt (A-10)

and similarly for acoustic quantities.

We assume that the source emits uniformly in 4π and that the ground is

totally absorbing; i.e., we ignore ground reflections.

A-2 Acoustic Equations

In the barred (comoving) frame, the relevant acoustic equations are the

usual ones:

∂tv ̄ai + ρ−1

0 ∂ ̄

ip ̄a = − 

T T /ρ0 T0



∂ ̄

ip ̄a − v ̄a∂ ̄v ̄T i − v ̄T ∂ ̄v ̄ai , (A-11)

∂tp ̄a + γh p0 ∂ ̄

iv ̄ai = −v ̄T i∂ ̄ip ̄a + Sa . (A-12)

The partial derivative symbols are ∂t = ∂/∂t and ∂i = ∂/∂xi and xi =

x, x2 = y, and x3 = z in the stationary frame. The ∂ ̄

r, etc, are partial derivatives in the comoving frame, where γh is the ratio of specific heat at constant

pressure to the specific heat at constant volume.

In these equations, we have dropped terms in p ̄T , since this is second order

in v ̄T , and used the equation of state to put ρ ̄T /ρ0 + T T /T0 =  ̄pT /ρ0 ≈ 0.

Also, we have included a scalar acoustic source function Sa to be chosen.

Note that the right-hand sides of these equations involve the turbulent temperature and velocity fluctuations T T ( ̄r, t), v ̄T i( ̄r, t); in this work, we regard

them as stationary, i. e., T T ( ̄r), v ̄T i( ̄r) only, with no explicit time dependence

in the barred frame. They may be stochastic functions. Thus, in every detail, the description of the turbulence and acoustics in the barred frame is

subjectively the same as in the standard description1 in any reference frame

with no wind.

The acoustic quantities in equation (A-11) and (A-12), the sum of the waves

created by the source and the waves scattered by the turbulence, are written

as

p ̄a =  ̄pae +  ̄pas , v ̄ai =  ̄vaei +  ̄vasi ,

in the Born approximation where the subscript e stands for “external,” i.e.,

for the source waves that act as externally applied fields on the turbulence,

and the subscript s stands for “scattered,” representing waves scattered by

the turbulence.

A-3 Source Wave

The “external” acoustic pressure and velocity waves in the barred (T)

frame, generated by a scalar source, satisfy equations (A-11) and (A-12)

with v ̄T i, T T



= 0:

∂tp ̄ae + γh p0 ∂ ̄

iv ̄aei = Sa( ̄r, t) (A-13a)

1

V. I. Tatarskii, The Effects of the Turbulent Atmosphere on Wave Propagation, Keter,

Jerusalem (1971).

30
∂tv ̄aei + 1

ρ0 ∂ ̄

ip ̄ae = 0 , (A-13b)

where Sa( ̄r, t) is the source term. From equation (A-13b), we may put

v ̄ae = ∂ ̄

iΦ ̄ → p ̄ae = −ρ0 ∂tΦ ̄ , (A-14)

where Φ ̄ is a “velocity potential.” Then from equation (A-13a),

−ρ0 ∂2

t Φ + ̄ γh p0∇2

Φ = ̄ Sa .

Dividing by γh p0 yields

(∇2 − 1

c2 ∂2

t )Φ = ̄ 1

ρ0c2 Sa , c2 ≡ γh p0 /ρ0 , (A-15)

where c is the adiabatic sound speed.

In the unbarred frame, the monochromatic isotropic point source is located

at r = RS = s and has frequency νs = ωs/2π:

Sa(r, t) = A

δ(r − s) exp(−iωst) + c.c. , (A-16)

where A is an amplitude and δ(r) is the Dirac delta. So, using equation

(A-1),

S ̄a( ̄r, t) = A

δ( ̄r − s + wt) exp(−iωst) + c.c. (A-17)

The Green function solution of equation (A-15) is

Φ = ̄ − A

ρ0 c2

exp(−iωst)

4π



d3r ̄ 1

R

δ( ̄r − s + wt − βR) exp(iksR ̄) + c.c. ,

(A-18)

where β ≡ w/c, and

R ≡  ̄r −  ̄r. (A-19)

We note that, with η ≡  ̄r − s + wt − βR,

∂

∂r ̄i

δ( ̄r − s + wt − βR) = ∂δ

∂η

∂η

∂r ̄i

=

∂δ

∂s

β

 

Ri (A-20)

∂tδ( ̄r − s + wt − βR) = ∂δ

∂η

∂η

∂t = −c

∂δ

∂s

β



.

Also, we neglect ∂ ̄

i



1/ R



= −Ri/R2 compared to (1/R)iks∂iRi = iks



Ri/R ,

since we will have ksR  1 always. So we have from equations (A-14), (A18), and (A-20)

p ̄ae( ̄r, t) = ρ0 c (iks + β · ∇s) Φ + ̄ c.c. , (A-21)

v ̄aei( ̄r, t)=(iks + β · ∇s)Ui + c.c. , (A-22)

31
where

Φ ̄ ≡ −

A

4πρ0 c2



d3r ̄ 1

R ̄ δ( ̄r − s + wt − β R) exp(iksR ̄) exp(−iωst) , (A-23)

Ui ≡ −

A

4πρ0 c2



d3r ̄



Ri

R

δ( ̄r − s + wt − β R) exp(iksR ̄) exp(−iωst) . (A-24)

The integrals in equations (A-23) and (A-24) can be done. Transform the

variables by

ηi =  ̄r

i − si + wit − βiR . (A-25)

Then

∂ηi

∂r ̄ = δi + βi



R , since ∂R

∂r ̄



= −∂R

∂r ̄

= −Rj . (A-26)

So the Jacobian is

J

 η

r



=1+ β · 

R . (A-27)

Then the integrals involve δ(η), so we get

Φ = ̄ − A

4πρ0 c2

exp(−iωst)

R + β · R exp(iksR), Ui = 

RiΦ ̄ . (A-28)

Here, from equation (A-25)

R =  ̄r −  ̄r =  ̄r − [η + s − wt + βR]η=0 ≡ Rs − βR , (A-29)

where

Rs ≡  ̄r − s + wt =  ̄r −  ̄s(t) (A-30)

is the instantaneous vector from the source to the observation point  ̄r. Then

taking R · R from equation (A-29) yields

R2 = Rs2 − 2β · RsR + β2R2, or



1 − β2

R2 + 2β · RsR − (Rs)

2 = 0 .

The solution is

R = γ2

	

−β · Rs +



1 − β2

 Rs2

+ β · Rs2




, (A-31)

where

γ2 = 

1 − β2−1 , as in relativity. (A-32)

32
Then

β · R = β · Rs − β2R . (A-33)

So

R + β · R = R

γ2 + β · Rs =



1 − β2

 Rs2

+ β · Rs2≡ Ds. (A-34)

Thus from equation (A-28),

Φ = ̄ − A

4πρ0 c2

exp(−iωst) exp(iksR)

Ds , Ui = − A

4πρ0 c2

exp(−iωst) exp(iksR)

Ds 

Ri . (A-35)

Then from equation (A-21),

p ̄ae( ̄r, t) = ρ0 c



iks + βj

∂

∂sj



Φ ̄ . (A-36)

But again we neglect (∂/∂s)(1/Ds

) compared to (iks/Ds)(∂/∂s)R. Now

from equations (A-30) and (A-31),

R



∂R

∂si

= −δi − β

∂R

∂si



= R∂R

∂si

= −Ri − β · R ∂R

∂si

,

so

∂R

∂si

= −Ri

R + β · R = − Ri

Ds . (A-37)

Thus from equations (A-35) and (A-36)

p ̄ae( ̄r, t) = − A

4πc

iks

(Ds)



1 − β · R

Ds



exp(−iωst) exp(iksR) + c.c. ,

but 1 − (β · R/Ds

) = R/γ2Ds, so finally

p ̄ae( ̄r, t) = − A

4πc

iksR

γ2(Ds)2 exp(−iωst) exp(iksR) + c.c. (A-38)

Again we have neglected (∂/∂si)(1/Ds

) compared to (iks/Ds)(∂R/∂si),

since ksR and ksDs  1 always.

Then from equation (A-22), we get similarly

v ̄aei( ̄r, t) = − A

4πρ0 c2

 iksRi

γ2(Ds)2



exp(−iωst) exp(iksR) + c.c. (A-39)

We will also need (∂/∂r ̄i) ̄pae( ̄r,t

 ̄). Again for the far-field condition ksR >>>

1, we have to find ∂R/∂r ̄i. But just as in equation (A-37) and above,

∂R

∂r ̄i

= Ri

Ds . (A-40)

33
Thus

∂p ̄ae

∂r ̄i

= A

k2

s

4πcγ2

RRi

(Ds)3 exp(−iωst) exp(iksR) + c.c. (A-41)

As β → 0, Ds → Rs

, R → Rs, γ → 1, and from equation (A-38), we get

p ̄ae( ̄r, t) = − A

4πc

iks

Rs exp(−iωst) exp(iksRs

) . (A-42)

Thus, we could put in general

−iksA

4πc

≡ Pea , (A-43)

where a = length, and Pe is thus the pressure amplitude of the wave from

the source at a distance Rs = a from the source in a still atmosphere. This

relation identifies the significance of the arbitrary amplitude A

.

A-4 Scattered Waves

In Born approximation, the scattered waves satisfy equations (A-11) and

(A-12) with ( ̄pas, v ̄asi) on the left-hand sides, but with ( ̄pae, v ̄aei) on the righthand sides. By combining these equations again to eliminate ∂ ̄

iv ̄asi and applying the usual solenoidal condition ∂ ̄

iv ̄T i = 0 on the turbulent velocity,

we get the Monin wave equation



∂ ̄

i∂ ̄i − c2∂2

t



p ̄as = −T −1

0 ∂ ̄

i



T T ∂ ̄

ip ̄ae

− 2ρ0 ∂ ̄

i∂ ̄



v ̄T v ̄aei , (A-44)

where two terms have canceled via equation (A-13a), because we assume

that S ̄a ≡ 0, where the turbulence is nonzero.

We write

p ̄as =  ̄pT

as +  ̄pvas , (A-45)

where p ̄T

as is due to the turbulent temperature fluctuations T T , the first term

on the right-hand side of equation (A-44), and p ̄v

as is due to the turbulent

velocity fluctuations v ̄T , the second term on the right-hand side. Then the

Green’s function solutions of equation (A-44) are

p ̄T

as( ̄r, t) = ∂

∂r ̄i

1

4πT0

 

d3r ̄ T

( ̄r)

| ̄r −  ̄r|



∂p ̄ae( ̄r

, t)

∂r ̄

i



t=t−| ̄r− ̄r|/c

(A-46)

p ̄v

as( ̄r, t) = ∂2

∂r ̄i∂r ̄

2ρ0

4π



d3r ̄ v ̄T

 ( ̄r

)

| ̄r −  ̄r|

v ̄aei 

 ̄r

, t − | ̄r −  ̄r|/c , (A-47)

where  ̄r is now the observation point. After doing the indicated derivatives,

we will put  ̄r = d − wt, d = RD, the detector location.

34
A-5 Temperature Scattering

From equation (A-41), we have



∂p ̄ae( ̄r

, t)

∂r ̄



t=t−| ̄r− ̄r|/c

= A

k2

s

4πcγ2

RRi

(Ds)3 exp(iksR) exp(−iωst) exp(iksRd

) + c.c. , (A-48)

where

Rd =  ̄r −  ̄r (A-49)

Ri = Rs

i − βiR , R = γ2 


−β · Rs + Ds (A-50)

Ds = Rs



(1 − β2) + β · R

s2

, (A-51)

but now,

Rs

i =  ̄r

i − si + wit − βiRd . (A-52)

In the expression for p ̄T

as of equation (A-46), when we take ∂/∂r ̄i, we again

neglect (∂/∂r ̄i)(1/Rd

) compared to (iks/Rd)(∂/∂r ̄i)



R ̄ + R ̄d



, since we will

have ksR ̄d >>> 1 also.

Similarly, we neglect (∂/∂)[ ̄ri



RRi



/(Ds

)3] also. So we will need

∂

∂r ̄i

exp iks



R ̄ + R ̄d

 = iks

∂

∂r ̄i



R ̄ + R ̄d



exp iks



R ̄ + R ̄d

 . (A-53)

But

∂R ̄d

∂r ̄i

= R

d

i , ∂R ̄

∂r ̄i

= ∂R ̄

∂R ̄s



∂R ̄s



∂r ̄i

= −β

∂R ̄

∂R ̄s





R

d

i . (A-54)

Using equations (A-50), (A-51), and (A-54), we get

∂R ̄

∂R ̄s



= 1

Ds




R ̄ s

 − βR ̄ = R ̄

Ds ,

whereby

∂

∂r ̄i

(R ̄ + R ̄d) = R

d

i



Ds − β · R ̄ s

 + β2R ̄

 

R

d

0R ̄

Ds . (A-55)

Then

∂

∂r ̄i

exp iks



R ̄ + R ̄d

 = iks





R

d

i R/ ̄ Ds



exp iks



R ̄ + R ̄d

 . (A-56)

So from equations (A-46), (A-56), (A-41), and (A-43), we get

35
p ̄T

as( ̄r,t

 ̄) = − ks

4πT0

Pea

γ2



d3r ̄ RRi

(Ds)3



R

d

i

R ̄d T T

 ̄r



exp[iks(R ̄ + R ̄d)] exp(−iωst) + c.c. (A-57)

Now, the turbule is localized at  ̄r = b ̄; that is,

T T

 ̄r



= T( ̄r − b ̄) = T (ξ) . (A-58)

Define

ξ =  ̄r − b ̄ , (A-59)

and change variables in the integrand. Discard the ξ terms in the R,  ̄ R ̄i,

etc, in the numerator and denominator of equation (A-57), because R ̄d 

ξ, R ̄  ξ. But we must keep terms of order ξ in the exponent. That is, we

write

exp iks



R ̄ + R ̄d

 ≈ exp 

iks



R ̄ + R ̄d



ξ=0 exp(−iK · ξ) , (A-60)

and this defines the vector K. Then we put

 ̄r = d − wt = d ̄(t),  ̄r = b ̄ + ξ . (A-61)

Then

R ̄ d = d − b − wt − ξ = RDT − ξ , (A-62)

R ̄d = RDT − R DT · ξ , (A-63)

where

RDT ≡ d − b − wt . (A-64)

Now we need R as a Taylor series in ξ. Define R ̄0 by

R ̄ (ξ) = R ̄0 + ξ · ∇bR .  ̄

Here, R ̄0 is R



ξ=0 = 

R



 ̄r=b ̄, with R defined by equations (A-49) to (A-51).

After some algebra, the previous equations yield

Ki = ks

	

RDT

i − ∂R ̄0

∂ ̄bi




, (A-65)

where here

∂R ̄0

∂ ̄bi

= 1

Ds



δi + RDT

i β



R ̄0 , (A-66)

36
with

R ̄s

 = R ̄T S − βRDT , (A-67)

R ̄0 = γ2 


−β · R ̄ s + Ds (A-68)

Ds =



1 − β2 R ̄s2 − β · R ̄ s 1

2

, (A-69)

R ̄0i = R ̄s

i − βiR ̄0 . (A-70)

Then for p ̄T

as ( ̄r = d − wt, t) (at the detector), we get from equation (A-57)

p ̄T

as(detector) = AT (t) exp(−iωst) + c.c. , (A-71)

where AT (t) is a complex-valued time-dependent amplitude,

AT (t) = − k2

s

4πT0

Pea

γ2

R ̄2

0R ̄ 0 ·R DT

(Ds)3RDT exp[iks(R0 + RDT )]

T(K) , (A-72)

with K given by equation (A-65) and the other vectors by equations (A-66)

to (A-70). Here, 

T(K) is the Fourier transform of T(ξ), given by



T(K) ≡



d3ξ exp(−iK · ξ)T(ξ) . (A-73)

This is accurate only if the source and detector are in the far-field region of

the turbule at any given instant. That is, we must have RT S  λs, RT S  a,

and RT S  a2/λs, where a is the scale length of the turbule, λs = 2π/ks,

and similarly for RDT .

This is also valid to all orders of β; computationally, we should use it, although probably only first order in β is needed, since we expect β 0.05

in most realistic situations, except in strong storms. Note that β = 0.01

yields w = 3.4 m/s ≈ 7.6 mi/hr, so β ≤ 0.05 → w ≤ 38 mi/hr, reasonable. But terms of order β2 in K and in exp[iks(R ̄0 + RDT )] might shift

the phases noticeably, since ks = ωs/c = (2π)(500)/340 = 9.24 m−1. It

would be accurate enough to keep only first-order terms in β in the term

R ̄0 (R ̄ 0 · R DT )/(D ̄s)3RDT , but that would save very little time in numerical

computation, so the higher order terms are kept.

For simplicity, we use a Gaussian turbule,

T(ξ) = (∆T) exp(−ξ2/a2) . (A-74)

The size or scale of this is a. Then the Fourier transform is



T(K) = π3/2a3(∆T) exp(−K2a2/4) . (A-75)

37
Consider the result equation (A-72) in the β → 0 limit. Then

AT (t) = 

Pea exp(iksRT S)

RT S 

 ! "



 −k2

s

4πT0

T(K) cos θ

 ! "







exp(iksr)

r



 ! "

, (A-76)

effective incident fT (θ) spherical outgoing

plane wave amplitude scattering amplitude wave from scatterer

where we have defined RDT ≡ r, R DT = ˆr, R T S = kˆs, and, for β → 0,

R0 = kˆRT S, Ds = RT S, R0 · R DT = cos θ ,

with θ the scattering angle. This A ̄T (t) agrees with equation (2) found in

an earlier paper by Goedecke and Auvermann;2 it is the Monin result for

acoustic scattering by stationary turbulent temperature inhomogeneities.

The ∆T is the temperature fluctuation amplitude. For isotropic homogeneous turbulence, it scales with a as |∆Ta| = |∆T1 |(a/a1 )

1/3

, as discussed

in the main text. If T0 is the mean temperature T, then half the ∆Ta must

be taken negative and the other half positive, but with the same magnitude

for each a.

A-6 Velocity Scattering

From equations (A-39), (A-52), and (A-47), we see that we need

∂

∂r ̄i∂r ̄j

exp[iks(R ̄ + R ̄d)] .

Previously, we found in equation (A-56)

∂

∂r ̄i

exp[iks(R ̄ + R ̄d)] = iks





R

d

i R/ ̄ Ds



exp[iks(R ̄ + R ̄d)] .

We obtain the derivative of this here by once again neglecting ∂R/∂  ̄ r ̄ compared to ksR ∂R/∂r ̄, valid since ksR >>> 1. So we get

∂2

∂r ̄i∂r ̄

exp[iks(R ̄ + R ̄d)] ≈ −k2

s



R

d

i 

R

d



R2

(Ds)2 exp[iks(R ̄ + R ̄d)] . (A-77)

So from equations (A-47) and (A-39),

p ̄v

as( ̄r,t

 ̄) = −k2

s

1

2πc 

d3r ̄







(Pea) R ̄i

γ2(D ̄s)2

 R

 ̄

d

i R

 ̄

d

j



R ̄2



(D ̄s)2 

R ̄d





 (A-78)

× v ̄tj ( ̄r

) exp[iks(R ̄ + R ̄d)] exp(−iωst) + c.c.

2

G. H. Goedecke and H. J. Auvermann, “Acoustic scattering by atmospheric turbules,”

J. Acoust. Soc. Am. 102 (1997), pp 759–771.

38
Now we put

v ̄t( ̄r

) ≡ v ̄



 ̄r − b ̄

≡ v ̄ (ξ) , (A-79)

and we put ξ = 0 in the terms in the square bracket in the integrand. So

there R

 ̄

d

i = R

 ̄

DT

i , Ri = R0i , R = R0 , etc, as in equations (A-66) to (A-70),

and also

exp[iks(R ̄ + R ̄d)] = exp[iks(R ̄0 + RDT )] exp(−iK · ξ) (A-80)

as in equation (A-60), where K is given by equation (A-65). Therefore, for

 ̄r = d ̄ − wt,

p ̄v

as (detector) =



− k2

sPea

2πc



R ̄0

3

γ2(Ds)4

 

R0 · R DT  R DT · v ̄(K)



RDT (A-81)

× exp[iks(R ̄0 + RDT )]

exp(−iωst) + c.c.

≡ Av(t) exp[−iωst] + c.c. ,

where Av(t) is defined by the brace. Here,

v ̄(K) ≡



d3ξ exp(−iK · ξ)v ̄ (ξ) . (A-82)

For simplicity, we use a Gaussian turbule

v ̄ (ξ)=(Ω × ξ) exp(−ξ2/a2) (A-83)

as in Goedecke and Auvermann’s paper. Thus

v ̄(K) = Ω × i∇K



d3ξ exp(−iK · ξ) exp(−ξ2/a2)

= ia3π3/2Ω × ∇K exp(−K2a2/4) or

v ̄(K) = −ia5π3/2Ω

2



Ω × K



exp(−K2a2/4) . (A-84)

For isotropic homogeneous turbulence, the energy transfer rate per unit

mass, ε, from one turbule size to the next satisfies

ε ∝ v2

(v/a) = (independent of a) ≈ Ω3a2 .

Thus Ω(a) scales with a as

Ω(a) = Ω(a1 )

a1

a

2/3

, a1 = amax = largest turbule scale length used. (A-85)

39
So we put

v ̄(K) = −ia5π3/2

2



Ω(a1 )a2/3

1



a13/3



Ω × K



exp(−K2a2/4) . (A-86)

In equation (A-84), Ω is a “randomly” oriented unit vector. If the turbulence

is assumed to be isotropic, then Ω has a uniform probability distribution in

4π; that is, if we write

Ω = e1 sin θΩ cos φΩ + e2 sin θΩ sin φΩ + e3 cos θΩ ,

then

Pθ (θΩ ) dθΩ = sin θΩ dθΩ

2 → Pθ (θΩ ) = sin θΩ2 ,

Pφ (φΩ ) dφΩ = dφΩ

2π → Pφ (φΩ ) = 12π ,

where Pθ (θΩ ), Pφ (φΩ )



are the probability distributions of (θΩ , φΩ ). If the

turbulence is not taken to be isotropic, then a spherically symmetric rotating turbule model may not be adequate, and also Ω is not uniformly

distributed in 4π.

Consider the result equation (A-81) in the limit β = 0. Insert b ̄ = 0 = b.

Then R T S = kˆ; R DT = ˆr = dˆ; R ̄ 0 = RT S; Ds = RT S; RDT = r; so we get

p ̄v

as =



Pea exp(iksRT S)

RT S 

 ! "



−k2

s

2πc

cos θ



ˆr · v ̄(K)



 ! "

 

exp (iksr)

r



 ! "

, (A-87)

effective incident fv (θ, φ) spherical outgoing

plane wave amplitude scattering amplitude wave from scatterer

where here

K = ks



ˆr − kˆ



. (A-88)

Note that K · v ̄(K)=0, since ∇ξ · v ̄ (ξ)=0, whereby ˆr · v ̄(K) = kˆ · v ̄(K)

for β = 0. This scattering amplitude is the same as equation (2) of Goedecke

and Auvermann’s paper. (This is the Monin result for scattering by turbulent velocity fluctuations.)

We note that K · v ̄(K)=0 in general, for any β, for the rotating sphere

model (see eq (A-82)). In fact, this must be true for any turbulent eddy that

satisfies ∇ξ · v ̄ (ξ)=0, because then



d3ξ exp(−iK · ξ)∇ξ · v ̄ (ξ)=0= iK · v ̄(K) , (A-89)

where integration by parts and discard of the surface integral at “spatial

infinity” yield the second equality. This surface integral must vanish since

each turbule is localized. (We assume that turbules are “wholly contained”

in the scattering volume.)

40
Appendix B. Computer Algorithm

B-1 Detector and Source Locations and Wind

The calculations are done in the “primed” coordinate system, which is

ground-fixed but in which the wind is along the x axis,

wx = w > 0, wy = wz = 0,

as in figure B-1. Define

xCM ≡

1

2

(xD + xS),

yCM ≡

1

2

(yD + yS).

The coordinate transformation is then

z = z,

x = (x − xCM) cos α + (y − yCM) sin α,

y = −(x − xCM) sin α + (y − yCM) cos α ,

(B-1)

and components of all position vectors are related in the same way. The

user inputs the Cartesian components of (RD, RS) and the magnitude w

and azimuthal angle α of the horizontal wind velocity w, in the “original”

(unprimed) frame. The algorithm calculates the primed components of all

vectors.

Figure B-1. “Primed”

coordinate system. y

w

x yD

D

α

yCM

xD

xCM

x ys s

S

x ́ Turbule

y ́

O ́

O

41
B-2 Spectral Analysis

Consider the Fourier transform of equation (8) from the main report,

p ̃T

as(ω) =  ∞

−∞

dt exp(iωt)



AT (t) exp(−iωst) + AT∗(t) exp(iωst)



or

p ̃T

as(ω) = A ̃T (ω − ωs) + A ̃T∗(ω + ωs) . (B-2)

Therefore, all we need to calculate is

A ̃T (ω) ≡

 ∞

−∞

dt exp(iωt)AT (t) (B-3)

and then use equation (B-2) if we wish. But the spectral broadening and

shifting is all contained in A ̃T (ω), and similarly for the relation of p ̃v

as(ω) to

A ̃v(ω).

Since we use fast Fourier transform algorithms, we must establish the minimum number of time-series points needed for adequate accuracy. First,

equation (42) from the main report shows that the frequencies in A ̃T (ω)

are given (approximately) by νsβ ·



R DT (t) − R T S(t)

. This varies from

a minimum of −2βνs to a maximum of 2βνs; so from the Nyquist criterion,

we need for A ̃T (ω)

νmax = 

(2∆t)

−1

> 2βνs → ∆t < (4βνs)

−1 . (B-4)

For example, for νs = 500 Hz, β = 0.01, we need ∆t 0.05 s → νmax 

10 Hz. The computer code uses an algorithm that yields these numbers for

νs = 500 Hz, β = 0.01. That is, it chooses

νmax = integer value ≥ 2βνs , ∆t = (2νmax)

−1 . (B-5)

Second, we need appropriate frequency resolution. We can obtain an estimate of TR (the minimum total time duration of a measurement that ensures adequate frequency resolution) by examining the scattered wave amplitude associated with large turbules with ksa  1. Consider figure B-2.

For what values x1 of x is the Gaussian reduced to perhaps e−4 of its maximum value? If ksa  1, x will be close to D/2. A little algebra yields

x1 = D

2

±

1

√2ksa

(D2 + 4h2)3/2

Dh . (B-6)

Now use equation (42) from the main report to obtain ν1 when x = D/2

and x = x1 . We easily get

ν1 ≡ νx=D/2 − νx=x1 =

√2w

πa D2 + 4h22Dh 

. (B-7)

42
Figure B-2. Estimating

minimum measurement

time duration.

z

D

x (D – x)

D x

h

RDT

w

Turbule

RTS

S

Essentially, only −ν1 ν ν1 will be in the Fourier transforms of the

amplitudes A(t).

As w/a gets smaller, so does ν1 . The quantity D2 + 4h2/2Dh has minimum

value 2 for h = D/2. So

(ν1 )min =

√2wmin

πamax

. (B-8)

Adequate frequency resolution requires perhaps ∆ν ≡ T −1

R ≤ 10−1 (ν1 )min.

We suppose that amax = 10 m, wmin = 1 m/s. Then

TR  10/ (ν1 )min = 10πamax √2wmin

≈ 200 s . (B-9)

Independently of this, we would like the total time T to be long enough to

contain almost all the scattered pulse from one turbule that passes “symmetrically” over the source/detector. For small turbules, the scattered pulse

amplitudes are approximately proportional to RDT RT S−1

. For source and

detector positions both near the ground and a turbule passing directly overhead, as in figure B-3, the scattered amplitude is essentially a maximum

when the turbule is at position 2, proportional to R−2

min, where

Rmin =



(D/2)2 + h2

 1

2

. (B-10)

At positions 1 and 3, the amplitude is proportional to R−2, where

R2 ≡ 

(R1 + D/2)2 + h2 1

2 

(R1 + D)

2 + h2 1

2 . (B-11)

For most scenarios, h  D for turbules of interest, so we could neglect h in

equations (B-10) and (B-11). Then we could require that we limit the length

of the scattering volume along the wind direction (which is along +x in this

simple analysis) by requiring R2  16R2

min; this would include all signals

within 6 percent of the maximum from a turbule that starts at x = −R1 and

travels the total distance 2R1 + D in time T. For h ≈ 0, this yields

(R1 + D)  2D → R1  D . (B-12)

43
Figure B-3. Estimating

maximum measurement time duration.

3

D x

w

h h

S

z

R1 D R1

h

1

D/2

2

Then, the total length L of a tube along the wind would satisfy

L = D + 2R1  3D = 3|RDS| . (B-13)

The computer algorithm uses the equality, and then chooses

T = L/w = 3|RDS|/w , (B-14)

or uses T = TR given by equation (B-9), whichever is larger. For most cases

of interest, equation (B-14) yields the larger value.

The minimum number of points N needed in the discrete time series for

the A(t) is then obtained by the use of equations (B-4) and (B-14),

N  T /∆t = 3RDS4βνs

w = 6π

ksRDS , (B-15)

or by the use of equations (B-4) and (B-9),

N  T /∆t = 10πamax2βνs √2w = 10

√2

ksamax , (B-16)

whichever is larger. Note that usually RDS  amax.

Since N is used to dimension arrays, the FORTRAN source code actually

defines a default value N = 6048 as a parameter. The user should verify

that this is large enough to satisfy equations (B-15) and (B-16), and if it is

not, change it and recompile. The values of T and L are calculated internally by putting (T = (N − 1)∆t, L = wT) .

B-3 Turbule Model Parameters

As shown by Goedecke and Auvermann,1 equations (67) and (68), for Kolmogorov spectra of isotropic homogeneous turbulence, the structure parameters 

C2

T , C2

v



are related to the turbule parameters by

C2

T = (3.78/μ) (φ)



∆T1 /a1/3

1

2

J8/3 , (B-17)

C2

v = (0.68/μ) (φ)



Ω1 a2/3

1

2

J14/3 ,

1

G. H. Goedecke and H. J. Auvermann, “Acoustic scattering by atmospheric turbules,”

J. Acoust. Soc. Am. 102 (1997), pp 759–771.

44
where successively smaller turbule sizes aα are related by

aα+1 /aα = e−μ ≡ N −1

d , (B-18)

where a1 = amax is the outer scale length, φ = n1 a1 = nα aα is the packing

fraction (the same for all sizes), and nα is the number density of turbules of

size α,(∆T1 , Ω1 ) are the (temperature fluctuations, angular velocity) of the

largest turbules, and here

Js ≡

 ∞

0

dx xse−x2= 1

2

Γ

s + 1

2



, (B-19)

where Γ is the gamma function. From tables of Γ functions,

J8/3 = 1

2

Γ

11

6



≈ 0.47 , (B-20)

J14/3 = 1

2

Γ

17

6



= 11

12

Γ

11

6



≈ 0.86 .

The (∆Tα , Ωα ) for size aα are scaled from (∆T1 , Ω1 ) as in equations (31)

and (32) from the main report. The parameter μ = ln Nd in equations

(B-14) and (B-15) is calculated after the inner scale length aNs , the smallest turbule size, and the number of sizes Ns in the cascade are chosen by

the user:

aNs /a1 = exp[−μ(Ns − 1)] → μ = −(Ns − 1)−1 ln(aNs /a1 ) . (B-21)

Default values are Ns = 6, a1 = 1.5 m, aNs = a6 = 4.6875 cm, which

yield μ = ln 2, Nd = 2, i. e., a halving of sizes with each step. Doubling

Ns while keeping (a1 , aNs ) unchanged should approximately double the

running time of the code. The effects of changing Ns, a1 , and/or aNs had

not been investigated at the time of this writing.

Values of C2

T , C2

v



may be input, in meter-kilogram-second units, and then

the appropriate values of (ΓT = ∆T1 /T0 , Γv = Ω1 a1 /c) are calculated. If the

values C2

T = C2

v = 0 are input, then the code sets default values of(∆T1 , Ω1 )

by putting

∆T1 /T0 = ΓT , Ω1 a1 /c = Γv , (B-22)

where (ΓT , Γv) may be input. The algorithm begins calculating with the

largest turbules (size a1 ). It places the first tube of turbules of square cross

section d1 = a1 /φ1/3 along the wind direction (x-axis) just on the y side

of the origin and centered along the x

-axis. The origin is chosen halfway

between the source and the detector. That is, the centerline of the tube is

y = d1 /2, z = RS(3) + d0 + d1 /2. (The base of this tube is calculated to be

d0 m above the height of the source or detector above ground, whichever is

greater.) The length of the tube is chosen by equation (B-13). Then at t = 0,

one turbule is placed randomly inside each cubical cell of side length d1 in

the tube, such that each turbule center can have (x

, y, z) displacements

45
from its cell center from 0 to ±d/4 with equal independent probabilities.

Then AT (t), Av(t)



are calculated and accumulated for all these turbules.

Any part of a turbule that leaves the tube (at x = L/2) during time T is

introduced at x = −L/2 + w∆t, the next time step, thus keeping the

frozen turbulence model in steady state.

Successive tubes for the largest turbules are chosen in such a way that tubes

chosen later are further from the origin, in general. For example, the second tube would have a centerline 

y =−d1 /2, z =RS(3) + d0 + d1 /2



. The

third and fourth tubes would be above the first and second, respectively,

with centerlines (y = ±d1 /2, z = RS(3) + d0 + 3d1 /2). The fifth and

sixth tubes, if used, would be one step further away in y

, with centerlines

y = ±3d1 /2, z = RS(3) + d0 + d1 /2. This selection continues until input

limits or default limits are reached.

For smaller turbules, successive tubes are chosen that fill all the largest

tubes used. This implies, for example, that if a2 = a1 /2, then there will

be four times as many tubes of size d2 = a2 /φ1/3 = d1 /2 as there are for

size d1 . It is impractical to calculate turbule by turbule for the smaller turbules, because there are so many of them. Instead, the time-shift algorithm

is employed for turbules with RT S, RDT   

ksd2 or k−1 s



, whichever is

larger, as described in the main report. The crucial values of RT S and RDT

for this determination are the values for which the maximum signal occurs.

46
Appendix C. Quick Doppler Formula

Refer to figure 1 of the main report (p 3) for the geometry. The turbule is

at RT (t) = b(t) = b0 + wt, where w = wind velocity. Pulses are emitted

from the source S located at position RS ≡ s at times t

S

1 , tS2 > tS1 with

t

S

2 − t

S

1 ≡ 2π

ωs , ωs = source angular frequency. Then the distances from the

source to the turbule at times (t

T

1 , tT2 ) are

RT S

1 = |b(t

T

1 ) − s| = |b0 − s + wt

T

1 | ,

RT S

2 = |b0 − s + wt

T

2 | = |RT S1 + w(t

T

2 − t

T

1 )|

= RT S

1 + w · R T S1 (t

T

2 − t

T

1 ) . (C-1)

Thus these pulses get to the turbule at times (t

T

1 , tT2 ) given by

t

T

1 = t

S

1 +

RT S

1

c(1 + β · R T S

1 )

≈ t

S

1 +

RT S

1

c (1 − β · R T S

1 ) ,

t

T

2 = t

S

2 +

RT S

2

c(1 + β · R T S

2 )

≈ t

S

2 +

RT S

2

c (1 − β · R T S

2 )

≈ t

S

2 + 1

c



RT S

1 + w · R T S1 (t

T

2 − t

T

1 )



(1 − β · R T S

1 ) , (C-2)

so

t

T

2 − t

T

1 = t

S

2 − t

S

1 + β · R T S1 (t

T

2 − t

T

1 ) ≈ (t

S

2 − t

S

1 )(1 + β · R T S1 ) . (C-3)

The distances from the turbule to a detector located at RD = d at times

(t

T

1 , tT2 ) are

RDT

1 = |d − b0 − wt

T

1 | ,

RDT

2 = |d − b0 − wt

T

1 − w(t

T

2 − t

T

1 )| ≈ RDT1



w · R DT

1

 

t

T

2 − t

T

1

 , (C-4)

so the times (t

D

1 , tD2 ) at which the scattered pulses are received by the detector D are

t

D

1 ≈ t

T

1 +

1

c

RDT

1 ,

t

D

2 ≈ t

T

2 + 1

cRDT

2 ≈ t

T

2 +

1

c

RDT

1 −



β · R DT

1

 

t

T

2 − t

T

1

 . (C-5)

47
Thus to order β,

t

D

2 − t

D

1 = 

t

T

2 − t

T

1



(1 − β · R DT

1 )=(t

S

2 − t

S

1 )[1 − β · (R DT

1 − R T S

1 )] . (C-6)

Therefore the frequency observed at the detector is

ωobs = 2π

tD

2 − tD1

≈ ωs(1 + β · (R DT

1 − R T S1 )) . (C-7)

This is ωobs at D at time t

D

1 = t

T

1 + 1

c |d − b0 − wt

T

1 | ≈ t

S

1 + 1

c



RT S

1 + RDT1



.

Both R DT

1 and R T S2 should be evaluated at time t

T

1 in this formula for ωobs .

To first order in β, this is the same result obtained in the main report directly

from the expressions for the waves pT

as(t), pvas(t)



at the detector.

48
Appendix D. Computer Codes

This appendix contains descriptions of four FORTRAN codes. They are

written in structured FORTRAN 77 standard form. The Microsoft International Mathematical and Statistical Library is used for random number

generation and calculation of fast Fourier transforms. The Microsoft portlib

routine TIMER() was used to time the codes. The codes are input.for used

with sloweddy.for and input1.for used with slowone.for.

The input.for and input1.for codes write parameters to files. The parameters in the files can be edited by a text editor without having to rerun the

input codes.

The description and code (with a list of symbols) for program input follow.

49
*******************************************************************************

* Program Input *

* Version 1.0: 11-21-98 *

*******************************************************************************

* This program enters the environment parameters, and the turbule *

* parameters needed for the Slow Eddy program. These parameters *

* are written to the file "input.txt." *

*******************************************************************************

PROGRAM INPUT

IMPLICIT NONE

INTEGER NMAX

PARAMETER(NMAX=6400)

CHARACTER CONFIG

INTEGER IT,NS,NT

REAL ALPHA,AP,A1,C,D(3),GAMT,GAMV,MW,NUS,PF

REAL S(3),YC1(20),ZB,ZC1(20)

OPEN(1,FILE=’input.txt’)

c The following inputs have to do with the environment for the turbules.

PRINT *

PRINT *, ’Enter the sound speed in meters per second:’

READ *, C

PRINT *

PRINT *, ’Enter the sound frequency of the source in Hertz:’

READ *, NUS

PRINT *

PRINT *, ’Enter the reference distance’

PRINT *, ’from the source in meters:’

PRINT *, ’(If not sure, enter 1)’

READ *, AP

PRINT *

PRINT *, ’Enter the x,y,z coordinates’

PRINT *, ’of the source in meters:’

READ *, S

PRINT *

PRINT *, ’Enter the x,y,z coordinates’

PRINT *, ’of the detector in meters:’

READ *, D

PRINT *

PRINT *, ’Enter the boundary layer height in meters:’

PRINT *, ’(higher than either source or detector)’

READ *, ZB

IF((ZB.LE.S(3)).OR.(ZB.LE.D(3))) THEN

PRINT *, ’The turbule(s) will pass through source or detector.’

PRINT *, ’Please run the input program again, and enter a’

PRINT *, ’boundary layer height greater than the source or’

PRINT *, ’detector height.’

STOP

ENDIF

50
PRINT *

PRINT *, ’Enter the wind speed in meters per second:’

READ *, MW

PRINT *

PRINT *, ’Enter the wind direction in degrees’

PRINT *, ’measured counter clockwise from the x-axis:’

READ *, ALPHA

c The following inputs have to do with the turbules themselves.

PRINT *

PRINT *, ’Enter the size of the largest turbule in meters:’

READ *, A1

PRINT *

PRINT *, ’Enter the packing fraction:’

READ *, PF

PRINT *

PRINT *, ’Enter the temperature scale parameter:’

READ *, GAMT

PRINT *

PRINT *, ’Enter the velocity scale parameter:’

READ *, GAMV

PRINT *

PRINT *, ’Diagrams for configurations are given below:’

PRINT *, ’(hit "enter" to see more)’

READ *

PRINT *

PRINT *, ’ Configuration choices ’

PRINT *

PRINT *

PRINT *, ’ ---- ’

PRINT *, ’a: | | ’

PRINT *, ’ ---- ’

PRINT *

PRINT *

PRINT *, ’ ---- ---- ’

PRINT *, ’b: | | | | ’

PRINT *, ’ ---- ---- ’

PRINT *

PRINT *, ’(hit "enter" to see more)’

READ *

PRINT *

PRINT *, ’ ---- ---- ’

PRINT *, ’ |||| ’

PRINT *, ’ ---- ---- ’

PRINT *, ’ ---- ---- ’

PRINT *, ’c: | | | | ’

PRINT *, ’ ---- ---- ’

PRINT *

PRINT *

PRINT *, ’ ---- ---- ---- ---- ’

PRINT *, ’d: | | | | | | | | ’

PRINT *, ’ ---- ---- ---- ---- ’

PRINT *

PRINT *, ’(hit "enter" to see more)’

51
READ *

PRINT *

PRINT *, ’ ---- ---- ’

PRINT *, ’ |||| ’

PRINT *, ’ ---- ---- ’

PRINT *, ’ ---- ---- ---- ---- ’

PRINT *, ’e: | | | | | | | | ’

PRINT *, ’ ---- ---- ---- ---- ’

PRINT *

PRINT *

PRINT *, ’ ---- ---- ---- ---- ’

PRINT *, ’ |||||||| ’

PRINT *, ’ ---- ---- ---- ---- ’

PRINT *, ’ ---- ---- ---- ---- ’

PRINT *, ’f: | | | | | | | | ’

PRINT *, ’ ---- ---- ---- ---- ’

PRINT *

PRINT *

PRINT *, ’g: User defined geometry:’

PRINT *

PRINT *

PRINT *, ’Enter the letter (a-g) corresponding to the’

PRINT *, ’the appropriate configuration desired.’

PRINT *, ’(make sure letter is lower case)’

READ *, CONFIG

PRINT *, ’Configuration chosen was ’,CONFIG

PRINT *

PRINT *, ’Enter the total number of turbule sizes:’

READ *, NS

c If configuration g is chosen, the user will need to enter

c the number a large tubes, and the y,z coordinates to the

c center lines of these tubes.

IF(CONFIG.EQ.’g’) THEN

PRINT *

PRINT *, ’Enter the number of large tubes:’

PRINT *, ’(maximum number allowed is 20)’

READ *, NT

DO 5 IT=1,NT

PRINT *

PRINT *, ’Enter the Y and Z coordinates of the center’

PRINT *, ’of large tube number’,IT

PRINT *, ’in dimensionless units:’

READ *, YC1(IT),ZC1(IT)

5 CONTINUE

ENDIF

c All parameters entered above are written to the file ’input.txt’

c used in the Slow Eddy FORTRAN program.

WRITE(1,*) C

WRITE(1,*) NUS

WRITE(1,*) AP

WRITE(1,*) S

WRITE(1,*) D

WRITE(1,*) ZB

WRITE(1,*) MW

WRITE(1,*) ALPHA

52
WRITE(1,*) A1

WRITE(1,*) PF

WRITE(1,*) GAMT

WRITE(1,*) GAMV

WRITE(1,*) CONFIG

WRITE(1,*) NS

WRITE(1,*) NT

DO 20 IT=1,20

WRITE(1,*) YC1(IT),ZC1(IT)

20 CONTINUE

CLOSE(1)

END

*-----------------------------------------------------------------------------*

*-------------------------------- List of Symbols --------------------------*

*-----------------------------------------------------------------------------*

*

* ALPHA : real; angle of wind relative to non-rotated x-axis.

* AP : real; ’a prime’, reference distance from the source.

* A1 : real; size of largest turbule.

* C : real; speed of sound.

* CONFIG : character; abbreviation for ’configuration’ used to choose

* predefined tube configurations.

* D(3) : real; position vector to the detector.

* GAMT : real; ’capital gamma sub T’.

* GAMV : real; ’capital gamma sub v’.

* IT : integer; index for largest tube, with range from 1 to NT.

* MW : real; ’magnitude of wind’.

* NMAX : integer; maximum number of time (and frequency) data points

* NS : integer; number of sizes of turbules.

* NT : integer; number of largest tubes.

* NUS : real; ’Greek letter nu sub s’, the frequency

* of the source.

* PF : real; ’packing fraction’.

* S(3) : real; position vector to the source.

* YC1(20) : real; y coordinate to center line of largest tube, dimensionless.

* ZB : real; height above ground where scattering volume begins.

* ZC1(20) : real; z coordinate to center line of largest tube, dimensionless.

Descriptions of Program Slow Eddy and its eight subroutines and a list of symbols follows:

*******************************************************************************

* Program Slow Eddy *

* Version 1.0: 11-21-98 *

*******************************************************************************

* This program calculates the pressure amplitudes as a function of time, *

* and their Fourier transforms as a function of frequency. This program *

* uses input from the file ’input.txt’ that is generated by running the *

* ’input.for’ Fortran program. The Slow Eddy program is designed to populate *

* the scattering volume with tubes of turbules. There are six predefined *

* configurations (a-f) and one configuration choice specified by the user (g).*

* The tubes are populated with NS different sizes of turbules, keeping the *

* packing fraction constant. Each smaller size is obtained by halving the *

* previous size. A separate program is provided to run one turbule by itself *

* or one tube of turbules, all the same size. It is called ’slowone.for’. *

* Subroutines: choise, direct, dot, multi, onetrb, output, rotate, tau. *

*******************************************************************************

53
PROGRAM MAIN

SUBROUTINE CHOICE(CONFIG)

c This subroutine sets parameters for predefined configurations

c chosen by the user in the ’input.txt’ file. It is assumed that

c configuration a has been run before configuration b, etc.

SUBROUTINE DIRECT(AMPTT,AMPVT)

c This subroutine calculates ’AMPTT’ and ’AMPVT’ one at a time and

c accumulates the results directly (i.e., no time-shifting).

c The ’MSIMSL’ below stands for ’Microsoft International Mathematical

c and Statistical Library’. MSIMSL contains a random number generator

c that the ’DIRECT’ subroutine uses.

SUBROUTINE DOT(V1,V2,V1V2)

c This subroutine calculates the dot product of two, real, three

c component vectors ’V1’ and ’V2’, and names this dot product ’V1V2’.

SUBROUTINE MULTI(AMPT,AMPV)

c The ’MULTI’ subroutine is a key subroutine in this program. It is named

c ’MULTI’ because it is what drives the calculations for multiple tubes

c with multiple turbules.

SUBROUTINE ONETRB(B0,AT,BV,K)

c The ’ONETRB’ subroutine calculates ’AT’, BV’, and ’K’

c for one turbule, as the turbule travels the total length of

c the tube. The ’ONETRB’ subroutine is used by both the ’DIRECT’

c subroutine, and the time-shifting ’TAU’ subroutine.

c The ’MSIMSL’ subroutine library is need for the random number generator.

SUBROUTINE OUTPUT(AMPT,AMPV)

c This subroutine takes the complex amplitudes ’AMPT’ and ’AMPV’

c as functions of time, and calculates the magnitude, and the real

c and imaginary parts. It also performs the FFT on them. Then the

c subroutine writes the data to files called ’time.txt’ and ’fft.txt’

SUBROUTINE ROTATE(ALPHA,MW)

c The ’ROTATE’ subroutine rotates the coordinate system so that the wind

c direction is along the x-axis. The origin is also moved to the point

c located half way between the source and detector.

SUBROUTINE TAU(AMPTT,AMPVT)

c The ’TAU’ subroutine employs the time-shifting procedure.

*-----------------------------------------------------------------------------*

*-------------------------------- List of Symbols --------------------------*

*-----------------------------------------------------------------------------*

* A : real; characteristic size of turbule.

* AAT(NMAX) : complex; temperature amplitude accumulation variable used

* in the ’CHOICE’ subroutine.

* AAV(NMAX) : complex; velocity amplitude accumulation variable used

* in the ’CHOICE’ subroutine.

* AF : real; ’angular frequency’ equal to 2*Pi*nu_s.

* ALPHA : real; angle of wind relative to non-rotated x-axis.

* AMPT(NMAX) : complex; amplitude due to temperature fluctuations.

* AMPTT(NMAX): complex; temperature amplitude accumulation variable for a

* tube used in the ’DIRECT’ and ’TAU’ subroutines.

* AMPV(NMAX) : complex; amplitude due to velocity fluctuations.

* AMPVT(NMAX): complex; velocity amplitude accumulation variable for a

54
* tube used in the ’DIRECT’ and ’TAU’ subroutines.

* AT(NMAX) : complex; temperature amplitude for one turbule used

* in the ’ONETRB’ subroutine.

* AT0(NMAX) : complex; temperature amplitude for one turbule used

* in the ’TAU’ subroutine in the far field.

* AV(NMAX) : complex; velocity amplitude for one turbule used

* in the ’ONETRB’ subroutine.

* AV0 : complex; velocity amplitude for one turbule used

* in the ’TAU’ subroutine in the far field.

* AP : real; ’a prime’, reference distance from the source.

* A1 : real; size of largest turbule.

* BCY : real; y coordinate of centerline of largest tube.

* BCZ : real; z coordinate of centerline of largest tube.

* BRK : complex; ’bracket’, a defined variable.

* BT : real; the location of the turbule as a function of time.

* BTA(3) : real; the Greek letter ’beta’. One of the

* IMSL library subroutines used ’beta’, so ’BTA’ was used

* as an alternate abbreviation.

* BTARS : real; the dot product of ’BTA’ and ’RS’.

* BTAR0 : real; the dot product of ’BTA’ and ’R0’.

* BV(NMAX,3) : complex; ’B sub v’, a defined variable.

* B0(3) : real; initial location of turbule.

* C : real; speed of sound.

* CALPHA : real; ’cosine of alpha’.

* COEF : real; ’coefficient’, a defined variable.

* CONFIG : character; abbreviation for ’configuration’ used to choose

* predefined tube configurations.

* CPHI : real; ’cosine of phi’.

* CROSS(3) : real; the cross product of ’RDTH’ and ’K’.

* CTEMP : complex; a temporary complex variable.

* CTHETA : real; ’cosine of theta’.

* D(3) : real; position vector to the detector.

* DA : real; distance between center of cells for size A.

* DCM(2) : real; position vector to detector relative to ’center of mass’

* position of source and detector.

* DI : integer; time-shifting increment index.

* DNS : real; distance between center of cells for smallest size.

* DS : real; ’D sub s’, a defined variable.

* DT : real; time step increment.

* DTMAX : real; maximum time step increment.

* DW(2) : real; x and y components of position vector to detector

* in the rotated ’wind’ frame.

* D1 : real; distance between center of cells for largest size A1.

* ETA(3) : real; Greek letter ’eta’, random displacement vector.

* FAMPT(NMAX): complex; Fourier transform of ’AMPT’.

* FAMPV(NMAX): complex; Fourier transform of ’AMPV’.

* FRQ(NMAX) : real; abbreviation for ’frequency’, the independent variable

* of the Fourier transform.

* GAMMA2 : real; abbreviation for ’lower case gamma, squared’,

* a defined variable.

* GAMT : real; ’capital gamma sub T’.

* GAMV : real; ’capital gamma sub v’.

* GAUSS : real; ’gaussian’, a variable that has a form that

* resembles a gaussian function.

* I : integer; time index, with range from 1 to N.

* IAMPT : real; ’imaginary part of AMPT’.

* IAMPV : real; ’imaginary part of AMPV’.

* II : complex; the imaginary number ’i’, the square root of -1.

* IM : integer; index for turbule in a tube, with range from 1 to M.

55
* INS : integer; number of time points in smallest cell.

* IS : integer; index for sizes with range from 1 to NS.

* IT : integer; index for largest tube, with range from 1 to NT.

* I1 : integer; number of time points in largest cell.

* J : integer; index for vector components, with range 1 to 3.

* J2 : integer; running index for y component of a vector.

* J3 : integer; running index for z component of a vector.

* K(NMAX,3) : real; a defined variable.

* KETA : real; the dot product of ’K’ and ’ETA’.

* KS : real; wave number for the source ’k sub s’.

* L : real; length of tubes.

* LMIN : real; minimum length that tubes can have.

* M : integer; number of turbules in a tube.

* MAMPT : real; ’magnitude of AMPT’.

* MAMPV : real; ’magnitude of AMPV’.

* MBTA : real; ’magnitude of BTA’.

* MFAMPT : real; ’magnitude of Fourier transform of AMPT’.

* MFAMPV : real; ’magnitude of Fourier transform of AMPV’.

* MK : real; ’magnitude of K’.

* MRDT : real; ’magnitude of RDT’.

* MRS : real; ’magnitude of RS’.

* MRTS : real; ’magnitude of RTS’.

* MR0 : real; ’magnitude of R0’.

* MW : real; ’magnitude of W’.

* M1 : integer; number of turbules in largest tube.

* N : integer; number of time (and frequency) data points used.

* NMAX : integer; maximum number of time (and frequency) data points

* that can be used.

* NMIN : integer; minimum number of time (and frequency) data points

* that can be used.

* NS : integer; number of sizes of turbules.

* NST : integer; number of tubes in largest cell after each halving.

* NT : integer; number of largest tubes.

* NUS : real; ’Greek letter nu sub s’, the frequency

* of the source.

* OMEGAH(3) : real; the unit vector ’omega hat’. Omega

* is rotation velocity vector of a turbule.

* OMEGA1 : real; scaling parameter for the velocity turbules.

* PF : real; ’packing fraction’.

* PHASE : real; a defined variable in a complex exponent.

* PI : real; variable for 3.14159...

* RAMPT : real; ’real part of AMPT’.

* RAMPV : real; ’real part of AMPV’.

* RDT(3) : real; vector from turbule to detector.

* RDTH(3) : real; ’RDT hat’, unit vector in direction of ’RDT’.

* RN(4) : real; a random number array.

* RNA(3) : real; a random number array for angles.

* RNT(2) : real; a random number array for temperature sign of + or -.

* RPHASE : complex; abbreviation for the random phase, a defined variable.

* RS(3) : real; a defined variable.

* RTS(3) : real; vector from source to turbule.

* RTSH(3) : real; abbreviation for ’RTS hat’, a unit vector in the

* direction of ’RTS’.

* R0(3) : real; a defined variable.

* R0H(3) : real; abbreviation for ’R0 hat’, a unit vector in the

* direction of ’R0’.

* R0RDTH : real; the dot product of ’R0H’ and ’RDTH’.

* S(3) : real; position vector to the source.

* SALPHA : real; ’sine of alpha’.

56
* SCM(2) : real; position vector to source relative to ’center of mass’

* position of source and detector.

* SIGMA2 : real; ’sigma squared’, a defined variable.

* SPHI : real; ’sine of phi’.

* STHETA : real; ’sine of theta’.

* SW(2) : real; x and y components of position vector to source

* in the rotated ’wind’ frame.

* T : real; time variable.

* TEMP : real; a temporary real variable.

* TIMER : real; gives value of total time for the program to run.

* TMAX : real; largest value of time variable.

* TS : real; time for calculations for one size to be completed.

* TSF : real; ’temperature scale factor’, a defined variable.

* TWOPI : real; ’2*pi’.

* VSF : real; ’velocity scale factor’, a defined variable.

* V1(3) : real; arbitrary vector used in the ’DOT’ subroutine.

* V2(3) : real; arbitrary vector used in the ’DOT’ subroutine.

* V1V2 : real; result of the dot product of ’V1’ and ’V2’ in the

* ’DOT’ subroutine

* W(3) : real; wind velocity vector.

* XDS : real; x component of the vector from the source

* to the detector.

* YC : real; y coordinate to center line of tube for size A.

* YC1(20) : real; y coordinate to center line of tube for size A1.

* ZB : real; height above ground where scattering volume begins.

* ZC : real; z coordinate to center line of tube for size A.

* ZC1(20) : real; z coordinate to center line of tube for size A1.

* ZI(3) : real; Greek letter ’zi’, random variable

* between +1/4 and -1/4.

Descriptions of Program Input1 and Program Slow One follow:

*******************************************************************************

* Program Input1 *

* Version 1.0: 11-21-98 *

*******************************************************************************

* This program enters the environment parameters, and the turbule *

* parameters needed for the slowone.for program. These parameters *

* are written to the file "input1.txt." *

*******************************************************************************

*******************************************************************************

* Program Slow One *

* Version 1.0: 11-21-98 *

*******************************************************************************

* This program calculates the pressure amplitudes as a function of time, *

* and their Fourier transforms as a function of frequency. This program *

* uses input from the file ’input1.txt’ that is generated by running the *

* ’input1.for’ Fortran program. The Slow One program is designed to run *

* one tube of turbules of one size, or just one turbule by itself. *

* Subroutines: choice, direct, dot, multi, onetrb, output, rotate, tau. *

**********************************************************************************

57
Subroutines and symbols for Program Input1/Program Slow One are almost identical with those of Program Input/Program Slow Eddy.

Flow charts for Program Slow Eddy and its eight subroutines are given in

figures D-1, D-2, D-3, D-4, and D-5.

Figure D-1. Flow chart

for main program.

End

Call

CHOICE

Start

Slow Eddy

Read

Input

File

Call

ROTATE

Call

DOT

58
Figure D-2. Flow chart

for subroutine CHOICE.

 - -

CONFIG ≠ a

CONFIG = a

Star

subroutine

CHOICE

Branch

on

CONFIG

CONFIG = b

CONFIG ≠ b Branch

on

CONFIG

CONFIG ≠ f

CONFIG = g

Call

MULTI

Call

OUTPUT

Call

MULTI

Call

OUTPUT

CONFIG = f

Branch

on

CONFIG

Call

MULTI

Call

MULTI

Call

OUTPUT

Return

 - -

Call

OUTPUT

59
Figure D-3. Flow chart

for subroutine MULTI.

No

Yes

Turbule size loop

Tube loop

No

Yes

No

Yes

Start

subroutine

MULTI

Is A

>

0.5 m?

Call

DIRECT

Call

TAU

Is IS

=

NS ? 

Is IT

=

NT ?

Return

Yes

60
Figure D-4. Flow chart

for subroutines

ROTATE, DIRECT, and

TAU.

Start

subroutine

ROTATE

Set wind

vector

Return

No

Yes

Start

subroutine

DIRECT

Call

ONETRB

Return

Call

DOT

Is IM

=

M ?

Turbules in a tube

loop

No

Turbules in a tube

loop

Start

subroutine TAU

Call

ONETRB

Return

Call

DOT

Is IM

=

M ?

Yes

61
Figure D-5. Flow chart

for subroutines

ONETRB, OUTPUT,

and DOT.

No

Start

subroutine

ONETRB

Calculate

temperature

and velocity

amplitudes

(Six calls to

DOT)

Time loop

Yes

Return

Is I =

N ?

Start

subroutine

OUTPUT

Write

out

time

data

Call Fourier

transform

routine

Write

out

spectral

data

Return

Start

subroutine

DOT

Multiply

components

and add

Return

Not shown in figure D-2 for subroutine CHOICE is the intermediate save

operations in branches for which the variable CONFIG is b through e. These

are different in each branch, with typical file names given in table D-1.

Table D-1. Intermediate

save operations. CONFIG File name(s)

b cf12.txt

c cf34.txt and cf1234.txt

d cf56.txt

e cf123456.txt*

*This branch does not call subroutine

MULTI.

62
63

Distribution

Admnstr

Defns Techl Info Ctr

Attn DTIC-OCP

8725 John J Kingman Rd Ste 0944

FT Belvoir VA 22060-6218

Mil Asst for Env Sci Ofc of the Undersec of

Defns for Rsrch & Engrg R&AT E LS

Pentagon Rm 3D129

Washington DC 20301-3080

Ofc of the Secy of Defns

Attn ODDRE (R&AT)

The Pentagon

Washington DC 20301-7100

Ofc of the Secy of Defns

Attn OUSD(A&T)/ODDR&E(R) R J Trew

3080 Defense Pentagon

Washington DC 20301-7100

AMCOM MRDEC

Attn AMSMI-RD W C McCorkle

Redstone Arsenal AL 35898-5240

ARL Chemical Biology Nuc Effects Div

Attn AMSRL-SL-CO

Aberdeen Proving Ground MD 21005-5423

US Army Dugway Proving Ground

Attn STEDP 3

Attn STEDP-MT-DA-L-3

Attn STEDP-MT-M Bowers

Dugway UT 84022-5000

US Army Missile Cmnd

Attn AMSMI-RD-AS-SS R Alongi

Redstone Arsenal AL 35898-5253

US Army OEC

Attn CSTE EFS

4501 Ford Ave Park Center IV

Alexandria VA 22302-1458

Dir for MANPRINT

Ofc of the Deputy Chief of Staff for Prsnnl

Attn J Hiller

The Pentagon Rm 2C733

Washington DC 20301-0300

Hdqtrs Dept of the Army

Assist Secy of the Army for RD&A

Attn SARD-DOV H Fallin Rm 3E411

The Pentagon

Washington DC 20301-0103

Natl Security Agency

Attn W21 Longbothum

9800 Savage Rd

FT George G Meade MD 20755-6000

Redstone Scientific Info Ctr

Attn AMSMI-RD-CS-R

Bldg 4484

Redstone Arsenal AL 35898

SMC/CZA

2435 Vela Way Ste 1613

El Segundo CA 90245-5500

TECOM

Attn AMSTE-CL

Aberdeen Proving Ground MD 21005-5057

US Army ARDEC

Attn AMSTA-AR-TD M Fisette

Bldg 1

Picatinny Arsenal NJ 07806-5000

US Army CECOM RDEC

Night Vsn & Elect Sensors Dirctrt

Attn AMSEL-RD-NV-OD F Milton

10221 Burbeck Rd Ste 430

FT Belvoir VA 22060-5806

US Army CRREL

Attn CRREL-GP R Detsch

72 Lyme Rd

Hanover NH 03755-1290

US Army Info Sys Engrg Cmnd

Attn AMSEL-IE-TD F Jenia

FT Huachuca AZ 85613-5300

US Army Natick RDEC Acting Techl Dir

Attn SBCN-T P Brandler

Natick MA 01760-5002
Distribution (cont’d)

64

US Army Simulation, Train, & Instrmntn

Cmnd

Attn AMSTI-CG M Macedonia

Attn J Stahl

12350 Research Parkway

Orlando FL 32826-3726

US Army Soldier & Biol Chem Cmnd Dir of

Rsrch & Techlgy Dirctrt

Attn SMCCR-RS I G Resnick

Aberdeen Proving Ground MD 21010-5423

US Army Tank-Automtv Cmnd Rsrch, Dev, &

Engrg Ctr

Attn AMSTA-TR J Chapin

Warren MI 48397-5000

US Army TRADOC Anlys Cmnd—WSMR

Attn ATRC-WSS-R

White Sands Missile Range NM 88002

US Army Train & Doctrine Cmnd

Battle Lab Integration & Techl Dirctrt

Attn ATCD-B

FT Monroe VA 23651-5850

US Military Academy

Mathematical Sci Ctr of Excellence

Attn MDN-A LTC M D Phillips

Dept of Mathematical Sci Thayer Hall

West Point NY 10996-1786

Nav Air War Cen Wpn Div

Attn CMD 420000D C0245 A Shlanta

1 Admin Cir

China Lake CA 93555-6001

Nav Surface Warfare Ctr

Attn Code B07 J Pennella

17320 Dahlgren Rd Bldg 1470 Rm 1101

Dahlgren VA 22448-5100

Naval Surface Weapons Ctr

Attn Code G63

Dahlgren VA 22448-5000

AF Rsrch Lab Phillips Lab Atmospheric Sci

Div Geophysics Dirctrt

Hanscom AFB MA 01731-5000

Phillips Laboratory

Attn AFRL-VSBE Chisholm

29 Randolph Rd

Hanscom AFB MA 01731-3010

DARPA

Attn S Welby

3701 N Fairfax Dr

Arlington VA 22203-1714

NASA Marshal Spc Flt Ctr Atmos Sci Div

Attn Code ED 41 1

Huntsville AL 35812

Armstrong Atlantic State University

Attn R C Wood (2 copies)

11935 Abercorn Street

Savannah GA 31419-1997

New Mexico State Univ Dept of Physics

Attn G Goedecke (2 copies)

Las Cruces NM 88003-8001

Univ of Mississippi NCPA

Attn H E Bass

University MS 38577

Dept of Commerce Ctr

Mountain Administration

Attn Spprt Ctr Library R51

325 S Broadway

Boulder CO 80303

Hicks & Associates Inc

Attn G Singley III

1710 Goodrich Dr Ste 1300

McLean VA 22102

Natl Ctr for Atmospheric Research

Attn NCAR Library Serials

PO Box 3000

Boulder CO 80307-3000

NCSU

Attn J Davis

PO Box 8208

Raleigh NC 27650-8208
65

Distribution (cont’d)

Director

US Army Rsrch Ofc

Attn AMSRL-RO-D JCI Chang

Attn AMSRL-RO-EN W D Bach

PO Box 12211

Research Triangle Park NC 27709

US Army Rsrch Lab

Attn AMSRL-DD J M Miller

Attn AMSRL-CI-AI-R Mail & Records Mgmt

Attn AMSRL-CI-AP Techl Pub (3 copies)

US Army Rsrch Lab (cont’d)

Attn AMSRL-CI-LL Techl Lib (3 copies)

Attn AMSRL-IS-EE D K Wilson

Attn AMSRL-IS-EE H J Auvermann

(10 copies)

Attn AMSRL-SE-EE Z G Sztankay

Attn AMSRL-SE-SA N Srour

Adelphi MD 20783-1197
1. AGENCY USE ONLY

8. PERFORMING ORGANIZATION

 REPORT NUMBER

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

12a. DISTRIBUTION/AVAILABILITY STATEMENT

10. SPONSORING/MONITORING

 AGENCY REPORT NUMBER

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

6. AUTHOR(S)

REPORT DOCUMENTATION PAGE

2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

11. SUPPLEMENTARY NOTES

14. SUBJECT TERMS

13. ABSTRACT (Maximum 200 words)

Form Approved

OMB No. 0704-0188

(Leave blank)

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources,

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this

collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson

Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503.

12b. DISTRIBUTION CODE

15. NUMBER OF PAGES

16. PRICE CODE

17. SECURITY CLASSIFICATION

 OF REPORT

18. SECURITY CLASSIFICATION

 OF THIS PAGE

19. SECURITY CLASSIFICATION

 OF ABSTRACT 20. LIMITATION OF ABSTRACT

NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102

Doppler Broadening of Acoustic Waves Scattered by

Turbulence Flowing With a Horizontal Wind

July 2000 Final, Oct 1997 to Aug 1999

This report investigates the scattering of a monochromatic acoustic wave by an assumed

unchanging atmospheric turbulence distribution that moves with a uniform constant horizontal

wind. The source and detector are at rest at different positions in a ground-fixed frame.

Theoretical expressions are derived for the scattered acoustic pressure at the detector. Because

distances and scattering angles change with time as any turbulent fluctuation (turbule) flows

through the scattering volume, the detector signal has a time-dependent Doppler frequency shift.

A simple formula for the consequent frequency broadening is derived. A computer code is

developed that calculates the scattered signal and its fast Fourier transform scattered by a single

turbule, or by a collection of turbules of many different scale lengths. The code uses a new

approximate but quite accurate time-shift algorithm that reduces the calculation time

substantially. Several numerical results from this code are presented for a model scenario that is

similar to a recent experimental scenario, with good agreement among spectral widths and

shapes.

Acoustic scattering, atmospheric turbulence, frequency spectrum

Unclassified

ARL-TR-2131

0FEJ60

61110253A11

B53A

61102A

ARL PR:

AMS code:

DA PR:

PE:

Approved for public release;

distribution unlimited.

78

Unclassified Unclassified UL

2800 Powder Mill Road

Adelphi, MD 20783-1197

U.S. Army Research Laboratory

Attn: AMSRL-IS-EP hauverma@arl.mil

U.S. Army Research Laboratory

67

George H. Goedecke (New Mexico State University), Roy

C. Wood (Armstrong Atlantic State University), Harry J.

Auvermann (ARL)

email:

2800 Powder Mill Road

Adelphi, MD 20783-1197
DEPARTMENT OF THE ARMY

U.S. Army Research Laboratory2800 Powder Mill RoadAdelphi, MD 20783-1197

An Equal Opportunity Employer