Interference of Echo-Signals from Two Buried Spherical Targets

: A numerically efﬁcient technique is presented for computing the backscattered ﬁelds from two spherical targets embedded in an underwater sediment. The bottom is assumed to be a homogeneous liquid attenuating half-space. The transmitter/receiver is located in a homogeneous water half-space. The distances between the transmitter/receiver and objects of interest are supposed to be large compared to the acoustic wavelengths in water and seabed. In simulations, the spherical scatterers of the same radius are assumed to be acoustically rigid. The interactions between two spheres are not taken into account because of the strong attenuation in the bottom. The scattering from one sphere in a wide frequency range is determined using the Hackman and Sammelmann’s general approach. The arising scattering coefﬁcients of the sphere are evaluated using the steepest descent method. The obtained asymptotic expressions for the scattering coefﬁcients essentially allowed to decrease a number of summands in the formula for the form-function of the backscattered acoustic ﬁeld.

In the present paper, the interference of echo-signals from two buried spherical scatterers is investigated. The bottom is assumed to be a homogeneous liquid attenuating half-space. The point source emitting the spherical incident wave of an angular frequency ω is placed at point M of a homogeneous water half-space. Modeling is performed in a wide frequency range of 40-60 kHz. The distance between scatterers and source/receiver is about 100 m.
Generally speaking, the two spheres may be of different sizes. However, for simplicity, later on in the discussions, the two spheres will be assumed to have the same radius.
The spherical scatterers of the same radius a are assumed to be acoustically rigid (i.e., homogeneous Neumann boundary conditions on the outer surface). The distance between their centers is d. The geometry of the problems is depicted in Figure 1.
The scattering from one sphere in a wide frequency range is determined using the Hackman and Sammelmann's general approach [19][20][21][22]. The arising scattering coefficients are evaluated using the steepest descent method [23,24].
An alternative and potentially attractive technique for computing the scattering coefficients of the sphere is the method of complex images [25][26][27][28][29][30]. The complex image solution is based on the approximation of the reflection and transmission coefficients via a -the radii of the spheres; -the distance between the centers of the spheres; -the distance from the centers of the spheres to the water/ground interface; , , -the horizontal and vertical distances between the centers of the spheres and the point source ; and -the angles corresponding to the mirror reflection path for the left sphere.
The scattering from one sphere in a wide frequency range is determined using the Hackman and Sammelmann's general approach [19][20][21][22]. The arising scattering coefficients are evaluated using the steepest descent method [23,24].
An alternative and potentially attractive technique for computing the scattering coefficients of the sphere is the method of complex images [25][26][27][28][29][30]. The complex image solution is based on the approximation of the reflection and transmission coefficients via a discrete sum of image point sources, with complex source point coordinates. The two main advantages of this method are its validity in the near field and the existence of a simple recursion relation in frequency [27] for computing the source parameters. The fitting of the transmission coefficient, which is needed in those configurations where the source is located in one of two media and the scatterer is located in the other medium, is affected by convergence problems which have been reported in the literature [28]. For this reason, the present work is focused on the steepest descent approximation of the scattering coefficients.
The goal of this paper is to investigate in detail the structure of the interference of the echo-signals from two or more spherical targets submerged in an underwater sediment and evaluate the influence of the Franz type surface wave and water/sediment interface reflections on the form-function.
Computing the far field scattered by objects inside an underwater sediment can be also performed via the Helmholtz-Kirchhoff integral (see, for instance, [31,32]).
In [6] it is shown that in the case of two spherical scatterers in water it is possible to neglect the multiple scattering between scatterers if 8 . When targets are buried in the attenuating bottom, the condition 8 is more than sufficient to neglect the multiple scattering between spheres. All computations and plots in this paper were performed using the computer algebra system Wolfram Mathematica. For special functions, such as spherical Hankel functions or spherical Bessel functions, the built-in functionality is applied. This paper is organized as follows. In Section 2, we propose a numerically efficient technique for computing the echo-signal from the only buried spherical target. In Section The goal of this paper is to investigate in detail the structure of the interference of the echo-signals from two or more spherical targets submerged in an underwater sediment and evaluate the influence of the Franz type surface wave and water/sediment interface reflections on the form-function.
Computing the far field scattered by objects inside an underwater sediment can be also performed via the Helmholtz-Kirchhoff integral (see, for instance, [31,32]).
In [6] it is shown that in the case of two spherical scatterers in water it is possible to neglect the multiple scattering between scatterers if d > 8a. When targets are buried in the attenuating bottom, the condition d > 8a is more than sufficient to neglect the multiple scattering between spheres.
All computations and plots in this paper were performed using the computer algebra system Wolfram Mathematica. For special functions, such as spherical Hankel functions or spherical Bessel functions, the built-in functionality is applied. This paper is organized as follows. In Section 2, we propose a numerically efficient technique for computing the echo-signal from the only buried spherical target. In Section 3, modeling of the scattered form-function of the only buried spherical target is carried out and discussed. In Section 4, we present the results of the numerical calculations of the scattering form-function for two and three spherical scatterers. Section 5 is the summary of relevant points. Finally, Section 6 concludes the paper.

Echo-Signal from the Only Buried Sphere
Throughout the entire paper, the attention is restricted to harmonic oscillations of constant frequency f . The complex time dependence e iωt , with t representing time, i = √ −1 and ω = 2π f , is factored out of equations.
For definiteness, let us consider the scatterer with a center at point O (see Figure 1). To evaluate the echo-signal from the buried spherical scatterer, let us use the method proposed in [19][20][21][22], where the acoustic potential of the echo-signal is represented in the form, In Equation (1), k = ω/c is the wave number in water, and T l are the elements of the free-field T-matrix for the acoustical scattering by the considered target. These elements are found using the separation of variables and applying the appropriate boundary conditions to the spherical harmonics. For the acoustically rigid sphere of a radius a embedded in an underwater sediment, where, k b = ω/c b is the wave number of a longitudinal wave propagating in the bottom, h l (x) is the spherical Hankel function of the 1st kind, j l (x) is the spherical Bessel function, and the prime marking the spherical functions in Equation (2) indicates a derivative with respect to the whole argument. The scattering coefficients of a sphere A ml (r) in Equation (1) are as follows: Here, ε 0 = 1, ε m = 2 for m ≥ 1; J m is the cylindrical Bessel function of the m-th order, q and h(q) = k 2 − q 2 are the horizontal and vertical components of the incident vector in water, respectively, and Π m l (x) are the normalized associated Legendre functions of the order l and rank m (see, for example, [33]) W(q) is the plane-wave transmission coefficient for propagation from water to bottom [24] Here, The densities of water and bottom are denoted as ρ and ρ b , respectively. The following inequalities are assumed to be satisfied in the complex q-plane: Imh(q) ≥ 0 and Imh b (q) ≥ 0.
In this paper, we will use the single-scatter approximation when C ml (r) = A ml (r). The full multiple scattering solution modeling the backscattered field from a thin air-filled spherical elastic shell in water close to the seabed or to the air-water interface was studied in [34].
If the echo-signal from the scatterer is evaluated taking into account the full multiple scattering, the coefficients C ml (r) in Equation (1) are found using a linear system of algebraic equations (see [19][20][21][22]34]). In fact, to find C ml , it is necessary at each fixed frequency to solve a system of l max + 1 equations with l max + 1 unknowns at m = 0, etc., and finally, a system consisting of one equation with one unknown C l max ,l max at m = l max . Coefficients of these systems are integrals with slowly decreasing and rapidly oscillating integrands. This makes computing the coefficients C ml (r) rather time consuming, and therefore in the present paper we will consider the single-scatter approximation.
The truncation level l max in Equation (1) is set by a rule suggested by Kargl and Marston [35], l max = ka + 4.05(ka) where, [x] is the integer part of x. For a = 0.3 m, c = 1500 m/s and f = 60 kHz, Equation (6) gives l max = 95. Thus, computing the backscattered field (1), it will be necessary to sum up more than 4500 summands. The integral representation of the scattering coefficients (3) is valid for arbitrary frequencies and distances between the source and the target. At frequencies of interest and distances y from 90 m up to 100 m, considered in this paper for illustration, the integrand in (3) is rapidly oscillating and slowly decreasing that makes the straightforward calculation of scattering coefficients rather time consuming. To speed up the computation of integrals (3), we will evaluate them using the steepest descent method (see, for example, [23,24]).
Let us change the variable in integral (3), putting q = ksin θ.
where, n = k b /k = c/c b is the refractive index at the water/bottom interface and χ = ρ b /ρ is the ratio of the densities of the bottom and water half-spaces.

Using properties of the Hankel functions H
(1) m and the associated Legendre functions P m l , the obtained integral (under the assumption that kysin θ 1) can be written in the form, where, the integration contour is going along the straight line Reθ = −π/2 from i∞ to 0, then along the real axis from θ = −π/2 to θ = π/2 and finally along the straight line Reθ = π/2 from 0 to −i∞; a ml is given using Equation (4), Since, by assumption, kysin θ 1, it is expedient to analyze the integral (8) using the saddle point method (see, for example, [23,24]). The saddle point is found from equation, which can be rewritten in the form, The incident angle θ is connected with the refraction angle θ b by the equality (the Snell's law), From Figure 1 it follows that, Acoustics 2023, 5

513
Thus, the geometrical interpretation of Equation (11) is the validity of equality, Let us denote the root of Equation (11) via θ * and the corresponding value of or, using Equation (11), where, LM and OL are the lengths of intervals forming the ray going from M to the origin of coordinates O, respectively.
Thus, the estimation of integral (8) by the saddle point method gives, Here, The angle θ * is found from Equation (11), and cos θ * b = 1 − sin 2 θ * /n 2 , χ = ρ b /ρ is the ratio of densities of the bottom and water half-spaces.
Taking account of the second approximation in the saddle point method will add to the expression (12) one more summand of the order 1/(ky) 2 . This summand will be only essential at small angles θ * . In this paper, we will consider parameters b = 0.5 m, z − b = 50 m, and y = 100 m, when Equation (11) gives θ * = 1.097 rad (θ * = 62.9 • ).
Our analysis of A ml (r) is complete if angle θ * is less than the angle of total internal reflection α = arcsinn. If θ * > α, it is necessary to take into account the two-valuedness of the function U(θ) (see Equation (7)) which enters into integral (8). Due to this twovaluedness, in the expression (12) for A ml (r), an additional term corresponding to the branch cut contribution appears.
Taking into account the representation (12) for the scattering coefficients A ml (r), Equation (1) for the acoustic potential of the echo-signal from the sphere will take the form, The summation theorem for the associated Legendre functions (see [33]) gives, Thus, The obtained formula for the acoustic potential of the echo-signal needs summation only with respect to l. The function B(r) is expressed via elementary functions. This formula is valid for any spherical scatterer: acoustically rigid, elastic, or elastic shells, or to use the proper elements of the T-matrix (for example, the T-matrix elements for a spherical elastic shell filled with air can be found in Appendix A of [36]).

Modeling of the Scattering Form-Function of the Only Buried Spherical Target
To study the frequency dependence of an echo signal, let us consider the form-function of acoustic scattering. In a free water space this function is determined as, where, Φ(ka) is the acoustic potential of the echo-signal at the reception point and Φ inc = exp(ikr)/(4πr) is the potential of the incident wave at the origin when there is no scatterer (see [21]). For the radiated spherical wave and the scattering geometry considered in this paper (see Figure 1) at n < 1, Φ inc is as follows (see [23,24] and the notations used in Equation (3)), Evaluating this integral as performed to obtain Equation (14), we get, Finally, for the form-function of the acoustic scattering we get the expression, In this paper, a medium model with the following parameters is considered: For the considered scattering geometry of the problem θ * = 1.097, α = arcsinn = 1.141. Thus θ * < α.
In Figure 2, the graph of the form-function (17) of the backscattered field from a spherical target of a radius a = 0.3 m immersed in a nonattenuating sandy sediment is shown. Since the sphere is considered as a rigid body, there is no wave propagating inside the sphere (no elastic effects) and the scattering problem reduces to a simple reflection problem in the exterior domain for the sphere. In this case, the echo-signal consists of three components: the mirror reflection (see the dashed lines in Figure 3a,b), the first Franz creeping wave which is exited at the sphere boundary between the illuminated and shadowed parts of the sphere (see the solid line in Figure 3a), and the wave which is shown by a solid line in Figure 3b. Additionally, the second Franz wave envelopes the spherical target one more time. As a result, its amplitude essentially decreases and the following Franz waves are not visible on the plot of the scattering form-function. The interference produces oscillations with a period equal to 1.336.  Similar computations were conducted in the case of the absorbing bottom. The bottom attenuation was taken into account by considering as a complex value with Im 4.46 m/s. As before, the form-function will be obtained using the formula (17), where is a complex value / 0.909 0.002 . In Figure 4, the graph of the form-function (17)    Similar computations were conducted in the case of the absorbing bottom. The bottom attenuation was taken into account by considering as a complex value with Im 4.46 m/s. As before, the form-function will be obtained using the formula (17), where is a complex value / 0.909 0.002 . In Figure 4, the graph of the form-function (17)   Similar computations were conducted in the case of the absorbing bottom. The bottom attenuation was taken into account by considering c b as a complex value with Imc b = −4.46 m/s. As before, the form-function F(ka) will be obtained using the Formula (17), where n is a complex value n = c/c b = 0.909 + 0.002i.
In Figure 4, the graph of the form-function (17) of the backscattered field from a spherical target of a radius a = 0.3 m immersed in an attenuating sandy sediment is shown. The decrease observed in the scattering form-function as a function of ka is connected with the factor exp{−k · Im f (θ * )} = exp{−(ka)q}, where, q = OL/a · Imn and OL = LQ/cos θ * b = 0.5/0.204 = 2.453 m. An approximate value of the form-function F(ka) obtained using (17) and an exact value calculated using (1)-(4) and (15) for attenuating sandy sediment were compared at several frequencies from the interval 40 ≤ f ≤ 60 kHz. At all these points, the difference between them was of the order of 10 −5 . Thus, the agreement of the exact and asymptotic formulae is excellent. An approximate value of the form-function obtained using (17) and an exact value calculated using (1)-(4) and (15) for attenuating sandy sediment were compared at several frequencies from the interval 40 60 kHz. At all these points, the difference between them was of the order of 10 . Thus, the agreement of the exact and asymptotic formulae is excellent.

Modeling of Interference of Echo-Signals from Two Spherical Scatterers
The form-function for two spherical scatterers of the same radius is naturally defined as (see [6]), Here, Φ is the acoustic potential of the sphere with a center at the point (see Fig

Modeling of Interference of Echo-Signals from Two Spherical Scatterers
The form-function for two spherical scatterers of the same radius a is naturally defined as (see [6]), Here,  An approximate value of the form-function obtained using (17) and an exact value calculated using (1)-(4) and (15) for attenuating sandy sediment were compared at several frequencies from the interval 40 60 kHz. At all these points, the difference between them was of the order of 10 . Thus, the agreement of the exact and asymptotic formulae is excellent.

Modeling of Interference of Echo-Signals from Two Spherical Scatterers
The form-function for two spherical scatterers of the same radius is naturally defined as (see [6]), Here, Φ is the acoustic potential of the sphere with a center at the point (see Figure 1), and Φ is the potential of the incident wave at the origin when there is no scatterer.
Let us consider two cases:  In Figure 6, the form-function (18) for two spherical targets immersed in an attenuating sediment is shown for d = 10 m (y = 100 m,ŷ = 90 m). The dashed line is the sum of the scattering form-function for the only sphere at y = 100 m (see Figure 4) and the scattering form-function for the only sphere atŷ = 90 m (see the dashed line in Figure 5). The dotted line is the absolute value of the difference between these two form-functions. The graph of the of the scattering form-function is an oscillating curve with a small period of oscillations. Therefore, we plot these graphs for 55 ≤ ka ≤ 60. In Figure 6, the form-function (18) for two spherical targets immersed in an attenuating sediment is shown for 10 m ( 100 m, 90 m). The dashed line is the sum of the scattering form-function for the only sphere at 100 m (see Figure 4) and the scattering form-function for the only sphere at 90 m (see the dashed line in Figure   5). The dotted line is the absolute value of the difference between these two form-functions. The graph of the of the scattering form-function is an oscillating curve with a small period of oscillations. Therefore, we plot these graphs for 55 60. In Figure 7, the form-function (18) for two spherical targets buried in an absorbing sandy sediment is shown for 5 m ( 100 m, 95 m). The dashed line is the sum of the scattering form-function for the only sphere at 100 m (see Figure 4) and the scattering form-function for the only sphere at 95 m (see solid line in Figure 5).
The dotted line is the absolute value of the difference between these two form-functions. In Figure 7, the form-function (18) for two spherical targets buried in an absorbing sandy sediment is shown for d = 5 m (y = 100 m, ∼ y = 95 m). The dashed line is the sum of the scattering form-function for the only sphere at y = 100 m (see Figure 4) and the scattering form-function for the only sphere at ∼ y = 95 m (see solid line in Figure 5). The dotted line is the absolute value of the difference between these two form-functions. In Figure 6, the form-function (18) for two spherical targets immersed in an attenuating sediment is shown for 10 m ( 100 m, 90 m). The dashed line is the sum of the scattering form-function for the only sphere at 100 m (see Figure 4) and the scattering form-function for the only sphere at 90 m (see the dashed line in Figure   5). The dotted line is the absolute value of the difference between these two form-functions. The graph of the of the scattering form-function is an oscillating curve with a small period of oscillations. Therefore, we plot these graphs for 55 60. line is the modulus of the difference between these two curves.
In Figure 7, the form-function (18) for two spherical targets buried in an absorbing sandy sediment is shown for 5 m ( 100 m, 95 m). The dashed line is the sum of the scattering form-function for the only sphere at 100 m (see Figure 4) and the scattering form-function for the only sphere at 95 m (see solid line in Figure 5).
The dotted line is the absolute value of the difference between these two form-functions.  line is the modulus of the difference between these two curves.
In Figure 8, the form-function for two spherical targets buried in an underwater sediment is shown for 5 m ( 95 m, 90 m). Based on the numerical calculations of the backscattered form-function for two spheres with different separation between two spheres and different distances , and , following conclusions may be appropriate: (1) Due to coherent addition of the two signals backscattered from two targets, the upper envelope of the two-sphere scattering form-function is the same as a sum of the formfunction curves for each of the two spheres. The lower envelope of the two-sphere scattering form-function is the same as the modulus of the difference between these form-function curves for each of the two spheres. (2) The oscillation period of the two-sphere scattering form-function expressed via is calculated as, For 100 m, this period increases about two times (from 0.2141 to 0.4257) when increases from 5 m to 10 m. In case does not change, but the bigger for the pair of spheres changes from 100 m to 95 m, the period of oscillation changes from 0.4251 to 0.4305.
(3) The influence of the sphere reflection and water/sediment interface reflection (see Figure 3a,b, respectively) is not visible on the plots of the form-function because this contribution is too weak to be detected on the curves.
The scattering form-function for three and more spherical scatterers of the same radius can be defined and calculated similarly to (18). For example, for three targets at 100 m, 95 m, and 90 m, the scattering form-function will be of the form shown in Figure 9. Based on the numerical calculations of the backscattered form-function for two spheres with different separation between two spheres and different distances y, ∼ y andŷ, following conclusions may be appropriate: (1) Due to coherent addition of the two signals backscattered from two targets, the upper envelope of the two-sphere scattering form-function is the same as a sum of the formfunction curves for each of the two spheres. The lower envelope of the two-sphere scattering form-function is the same as the modulus of the difference between these form-function curves for each of the two spheres. (2) The oscillation period of the two-sphere scattering form-function expressed via ka is calculated as, 2πa For y = 100 m, this period increases about two times (from 0.2141 to 0.4257) when d increases from 5 m to 10 m. In case d does not change, but the bigger y for the pair of spheres changes from y = 100 m to y = 95 m, the period of oscillation changes from 0.4251 to 0.4305.
(3) The influence of the sphere reflection and water/sediment interface reflection (see Figure 3a,b, respectively) is not visible on the plots of the form-function because this contribution is too weak to be detected on the curves.
The scattering form-function for three and more spherical scatterers of the same radius a can be defined and calculated similarly to (18). For example, for three targets at y = 100 m, ∼ y = 95 m, andŷ = 90 m, the scattering form-function will be of the form shown in Figure 9. When the distance between the centers of two spherical scatterers increases, the values of the form-function decrease. When the distance increases to some extent, the three-sphere form-function will be close to the form-function of the nearer spherical target from the source/receiver.

Discussion
An efficient method for computing the backscattered field from two or more targets immersed in an absorbing underwater sediment has been presented. In simulations, the targets are assumed to be spherical of radius and acoustically rigid. The bottom is assumed to be a homogeneous liquid attenuating half-space. The transmitter/receiver is located in a homogeneous water half-space. The distances between the transmitter/receiver and objects of interest are supposed to be large compared to the acoustic wavelengths in water and sediment.
For calculation of the echo-signal from the only buried sphere, we have followed the Hackman and Sammelmann's general approach. The arising scattering coefficients of the sphere were evaluated using the steepest descent method. The use of the obtained asymptotic expansions also allowed to essentially decrease the number of summands in the formula for the form-function of the backscattered acoustic field. The computational examples demonstrate efficiency of the presented technique.
The obtained asymptotic formula for the form-function can be used for the elastic target or the spherical elastic shell as well. In this case, it will only be necessary to replace elements of the free-field T-matrix.
In the case of the only acoustically rigid spherical scatterer, the echo-signal consists of three components: the mirror reflection, the Franz creeping wave which is exited at the sphere boundary between the illuminated and shadowed parts of the sphere, and the wave reflected from the water/non-absorbing sand interface. Taking into account the attenuation in the sediment, provides results depending essentially on the geometry of the problem (on depth of the scatterer and the length of the ray path propagating in a seabed).
The two spherical scattering form-function is an oscillation curve with a small period of oscillations. The oscillation period of the two sphere form-function expressed via is obtained using (19). The upper envelope is the same as a sum of the form-function curves The plot of the scattering form-function for these three spherical targets immersed in an attenuating sandy sediment becomes more complicated, i.e., quasiperiodic. The value maxF(ka) increases in comparison with the analogous values for the plots shown in Figures 6 and 7.
When the distance d between the centers of two spherical scatterers increases, the values of the form-function decrease. When the distance d increases to some extent, the three-sphere form-function will be close to the form-function of the nearer spherical target from the source/receiver.

Discussion
An efficient method for computing the backscattered field from two or more targets immersed in an absorbing underwater sediment has been presented. In simulations, the targets are assumed to be spherical of radius a and acoustically rigid. The bottom is assumed to be a homogeneous liquid attenuating half-space. The transmitter/receiver is located in a homogeneous water half-space. The distances between the transmitter/receiver and objects of interest are supposed to be large compared to the acoustic wavelengths in water and sediment.
For calculation of the echo-signal from the only buried sphere, we have followed the Hackman and Sammelmann's general approach. The arising scattering coefficients of the sphere were evaluated using the steepest descent method. The use of the obtained asymptotic expansions also allowed to essentially decrease the number of summands in the formula for the form-function of the backscattered acoustic field. The computational examples demonstrate efficiency of the presented technique.
The obtained asymptotic formula for the form-function can be used for the elastic target or the spherical elastic shell as well. In this case, it will only be necessary to replace elements of the free-field T-matrix.
In the case of the only acoustically rigid spherical scatterer, the echo-signal consists of three components: the mirror reflection, the Franz creeping wave which is exited at the sphere boundary between the illuminated and shadowed parts of the sphere, and the wave reflected from the water/non-absorbing sand interface. Taking into account the attenuation in the sediment, provides results depending essentially on the geometry of the problem (on depth of the scatterer and the length of the ray path propagating in a seabed).
The two spherical scattering form-function is an oscillation curve with a small period of oscillations. The oscillation period of the two sphere form-function expressed via ka is obtained using (19). The upper envelope is the same as a sum of the form-function curves from each of the two spheres. The lower envelope of the two-sphere form-function is the same as the modulus of the difference between these form-function curves for each of the two spheres.
At d = 10 m and 5 m, the periods calculated using this formula and the oscillation periods of the curves in Figures 5 and 6 coincide with the accuracy of 10 −6 .
Simulations held for three spherical targets showed more complicated (quasiperiodic) behavior of the plot of the scattering form-function, especially if the distance between two neighboring scatterers is not constant.

Conclusions
This paper has presented an efficient method for computing the backscattered field from two or more spherical targets immersed in an absorbing underwater sediment. The scattering from one sphere in a wide frequency range is determined using the Hackman and Sammelmann's general approach.
This technique is based on the evaluation of the arising scattering coefficients of the sphere via the steepest descent method. The efficiency comes from the availability to significantly reduce the number of summands in the formula for the form-function, as well as to substitute the integrals, which have slowly decreasing and rapidly oscillating integrands, with the expression that only involves elementary functions.
An extension to include more general background media, such as stepwise stratified media, should be fairly straightforward by adapting the method proposed by Hackman and Sammelmann.