c
TU Delft, The Netherlands, 2006
NUMERICAL SOLUTION OF NONLINEAR ACOUSTIC WAVE
PROBLEMS EMPLOYING A GREEN’S FUNCTION APPROACH
Jacob Huijssen∗ and Martin D. Verweij∗
∗Delft University of Technology, Laboratory of Electromagnetic Research
Mekelweg 4, 2628 CD Delft, The Netherlands e-mail: j.huijssen@tudelft.nl
Key words: Nonlinear acoustics, medical ultrasound, Neumann iterative method, Green’s func-tion method, convolufunc-tion integral, Burgers’ equafunc-tion
Abstract. The design of phased array transducers for medical diagnostic ultrasound asks for
an understanding of the nonlinear propagation of acoustic wavefields. Most existing numerical models are based on the linearized model equations, but in the recent decades several numerical models have been developed that incorporate weak nonlinear propagation as well.
In this study, we present an approach that enables the computation of large-scale, three-dimensional, weakly nonlinear wavefields in the time domain. It is based on the Neumann iterative method, which enables the successive use of the solution of a linear wave problem, where the nonlinearity is treated as a contrast source. The linear wave problem is formulated as a convolution integral with a free space Green’s function kernel. If the Green’s function is adequately regularized, this method enables us to discretize the spatiotemporal domain with step sizes up to the limit given by the Nyquist-Shannon sampling theorem for a bandlimited signal. The remaining spatiotemporal derivatives are handled with high-order finite difference schemes.
The proposed method is validated with a one-dimensional nonlinear wave problem. We evaluate the wavefield, including harmonic frequencies up to the fifth harmonic. The results are compared with a solution of the lossless Burgers’ equation. It is observed that already after a small number of iterations the computed result shows perfect agreement with those of the Burgers’ equation.
In this paper, we make plausible that the proposed method can be straightforwardly extended to more complex problems. The volume integral equation approach opens the road to solving time domain nonlinear wave problems in three dimensions. Furthermore, the method can be adapted to media with attenuation and inhomogeneity.
1 INTRODUCTION
score of models has been developed that predict the continuous-wave or pulse-time wavefields that are generated by phased array transducers with arbitrarily steered, focused and apodized excitations1,2. These models operate under a linear approximation of the basic acoustic equa-tions, thus neglecting the nonlinear character of acoustic propagation due to the convective term and due to excitation-dependent medium behavior. These effects are indeed of second order and under small-signal conditions they can be safely neglected. However, with the sound pressure levels that are common in echography, distortion of the pulse shape due to nonlinear propagation is clearly observable3.
Moreover, a new imaging modality called Tissue Harmonic Imaging (THI) has found its way into the medical diagnostic ultrasound practice that specifically benefits from this effect4. The nonlinear distortion of the transmit pulse during propagation through a tissue-like medium results in the generation of higher harmonic frequencies. With respect to traditional fundamental imaging (FI), the selective imaging of the second harmonic frequency, which is done in THI, benefits from an improved beam shape and thus in an increase in image quality. In a wide area of medical applications, THI has shown to give significantly improved imaging results when compared with FI5,6. Recently, it was suggested that selective imaging of even higher harmonics will result in further improvement7.
As THI is explicitly based on nonlinear wave propagation, the existing linear models can-not account for it. In the recent years, several approaches have been developed to account for the nonlinear distortion. A number of models is based on a forward-propagating wave approximation8,9, other models have attempted to solve the full-wave problem based on a sec-ond order nonlinear partial differential equation10,11or a set of first order basic equations12. The applied numerical techniques vary from finite difference, finite element, and pseudospectral to angular spectrum methods. The nonlinear distortion is accounted for in a number of ways:
• a plane-wave nonlinear propagator based on the Burgers’ equation8,9 • an iterative update scheme10
• a direct discrete solution of the nonlinear basic equations11,12 • an incrementally linear update scheme13
In order to model the nonlinear acoustic field of a phased array transducer propagating in tissue, we would need a model that can handle a large-scale, 3-dimensional configuration in the time domain. To reduce the vast computational complexity of this, most of the previ-ously mentioned models have been derived for problems with reduced dimensionality, e.g. for continuous-wave problems or for problems exhibiting a cylindrical symmetry. However, the challenge remains to develop a full-wave, nonlinear wave propagation model that can handle a large scale 3-dimensional configuration in the time domain, with an acceptable cost in terms of memory and computation time. The most recent contributions to address this challenge are presented in Refs. 14-16, all of which use the plane-wave nonlinear propagator to incorporate the nonlinearity.
to the linear wave problem. By correcting the linear wave problem by means of a Neumann iterative method, we obtain the solution to the nonlinear wave problem with any desired degree of accuracy. An equivalent approach was used in Ref. 10. For the linear wave problem, we use the Green’s function method. If the Green’s function is adequately regularized, we can ob-tain accurate wavefield results with a discretization up to the Nyquist criterion for the smallest wavelength of interest, i.e. with two points per wavelength.
In the subsequent sections, the Neumann iterative method is outlined and the application of the Green’s function method is discussed. It is shown how the Green’s function can be regularized to allow an almost optimal discretization. Next, the method is demonstrated with a planar nonlinear lossless wave propagation problem. We compare the results with an implicit solution of the Burgers’ equation using an excitation that is representative for medical diagnostic applications. It is observed that already after a small number of iterations the model gives accurate results. In the closing section we discuss future extensions of the method.
2 PROBLEM FORMULATION AND NEUMANN ITERATIVE METHOD
The nonlinear acoustic wave propagation can be described by well known second order par-tial differenpar-tial equations (PDE’s) like the Westervelt equation, the KZK equation and the Burg-ers’ equation17. Within the same order of approximation, an equivalent but more basic formu-lation is the set of first order PDE’s18
∇p + ρDtv = f , (1)
∇· v + κDtp = q, (2)
wherep is the acoustic pressure, v is the fluid velocity, and f and q denote the external force and volume source densities. Dt = ∂t+ v · ∇ is the total time derivative, and ρ and κ are the field-dependent mass density and compressibility, given as [19, 20]
ρ = ρ0[1 + κ0p], (3)
κ = κ0[1 + κ0(1 − 2β)p]. (4)
Here,β is the nonlinearity parameter and ρ0andκ0 are low-amplitude quantities. As a key step of our scheme, we rewrite the set of equations that result from Eqs. (1) to (4) by gathering the nonlinear terms in a contrast source term on the right hand side. The resulting equations can be cast in the general matrix form
DF =DL
x + MLDLt F = S + SN(F ), (5)
where DLx, DLt denote the linearized differential operators in space and time, and ML repre-sents the linear medium behavior. In the present case, the field variables, external sources and nonlinear contrast source terms are written as
We assume that the contrast source term SN(F ) will be small with respect to the other terms. If we can obtain a linear operator D−1that solves F in terms of S and SN, i.e.
F = D−1S + SN(F ) ,
(7) then we can write an iterative scheme to approximate the field solution of the nonlinear problem as
F0 = D−1[S] , (8)
Fj = D−1S + SN(Fj−1) . (9)
This represents our Neumann iterative method.
3 GREEN’S FUNCTION METHOD AND EFFICIENT DISCRETIZATION
One method to solve a linear PDE, i.e. to evaluate the linear operator D−1 in Eq. (7), is the Green’s function method. With this method, the PDE can be translated to an integral representation of the form
F = G ∗x,tS, (10)
where G denotes the Green’s function for the specific wave problem, S accounts for the rel-evant source terms and ∗x,t denotes the convolution with respect to the spatial and temporal dimensions.
We study the efficient approximation and evaluation of the convolution integral. For clarity, let us restrict ourselves in this paragraph to a single-variable, scalar convolution integral, given as
h(t) = Z
R
K(t − τ )f (τ )dτ, (11)
where the kernel functionK(t) can be non-singular as well as singular. We approximate f (t) on a basis of functionsφn(t). This gives an approximation for h(t) as
h∗(t) = Z R K(t − τ )X n fnφn(τ )dτ. (12)
The set of functionsφn(t) is defined by a single function φ(t) with φn(t) = φ(t − n∆t), where ∆t is a certain step size. The function φ(t) is normalized such that R
Rφ(τ )dτ = ∆t. The expansion coefficientsfnare determined by
lim ∆t↓0 Z R fnφn(τ )dτ = lim ∆t↓0 Z ∆t 2 −∆t 2 f (τ − n∆t)dτ, (13)
which for regular functionsf (t) and φ(t) is realised by taking fn = f (n∆t). With this, we can write Eq. (12) as
h∗(t) = ∆tX n
where[K](t) is the regularized form of K(t), [K](t) = 1 ∆t Z R K(t − τ )φ(τ )dτ. (15)
When f (t) is of limited support, say t ∈ [−∆t/2, N ∆t − ∆t/2], it follows that fn = 0 for n < 0 and n ≥ N . And when we restrict our interest to finding h(t) on the points t = m∆t for 0 ≤ m < N , we end up with a finite convolution sum approximation of Eq. (11),
h∗ m = ∆t N −1 X n=0 [K]m−nfn, 0 ≤ m < N, (16) whereh∗
m = h∗(m∆t) and [K]m−n = [K](m∆t − n∆t). This convolution sum can efficiently be evaluated with algorithms related to the Fast Fourier Transform (FFT)21,22.
The choice of a specific expansion functionφ(t), together with the step size ∆t, determines the accuracy of the approximation. For instance, withφ(t) = δ( t
∆t) we end up with a sampling approach, and Eq. (16) becomes the midpoint rule approximation of Eq. (11). The accuracy of the sum now depends on∆t and the smoothness of K(t) and f (t) in terms of their magnitude and derivatives. However, for the modeling of nonlinear acoustic propagation we would like to focus on the spectral content of the wavefield. We are interested in the accurate reproduction of the field up to a certain angular frequencyωmax. If we would like to obtain an approxima-tion that is most efficient in terms of the step size, then the absolute minimum is given by the Nyquist-Shannon theorem for bandlimited signals, that prescribes that at least two points per wavelength, for the smallest wavelength of interest, should be taken. An expansion function that is appropriate for bandlimited input functions is the sine cardinal function
φ(t) = sinc(πt/∆t) = (
sin(πt/∆t)
πt/∆t , t 6= 0
1, t = 0. (17)
Iff (t) is a function whose spectrum is bandlimited to |ω| < ωmaxand if∆t ≤ π/ωmax, then the result of Eq. (16) is exact, i.e. h∗
m ≡ h(m∆t). If f (t) is not bandlimited, then we can exactly reproduce the spectrum ofh(t) for |ω| < π/∆t with Eq. (16) when the coefficients fnare taken from a form of f (t) for which the frequency content for |ω| ≥ π/∆t has been removed. In signal processing terms this is equal to filteringf (t) with an ideal low-pass filter with a cutoff angular frequencyω = π/∆t before it is sampled.
1. Obtain the regularized form ofK(t).
2. Window it to a region[−(N − 12)∆t, (N − 12)∆t] (which is equal to multiplication with a box function).
3. Sample the result.
4. Obtain the DFT coefficients by means of an FFT.
But as an alternative and approximate procedure, we can perform the following steps: 1′. Obtain the continuous transform domain equivalent ˆK(ω) of K(t).
2′. Perform the windowing operation in the transform domain.
3′. Obtain the regularized form in the transform domain (which for the sinc expansion func-tion is equal to multiplicafunc-tion with a box funcfunc-tion).
4′. Sample the result to obtain the desired DFT coefficients.
Both procedures are equivalent except for the interchange in order of the windowing and reg-ularization operation. This gives an error that is small as long as the window size is large with respect to the step size, i.e. N ≫ 1. In the alternative procedure, the difficulty may lie in performing the windowing operation in the Fourier domain. This will be worked out for a one-dimensional wave problem in Section 5.
4 PARTIAL DERIVATIVES IN THE CONTRAST SOURCE TERM
In the determination of the contrast source terms SN(F ) in Eq. (6), a number of temporal and spatial derivatives of the wavefield quantities are needed. When the wavefield is available on a coarse grid, we can revert to spectral difference schemes or high-order finite difference (FD) schemes to obtain these with sufficient accuracy23. For a non-periodic function on an equidistant grid, spectral differences will result in Gibb’s phenomena at the grid boundaries. To avoid this, a high-order finite difference scheme can be applied with single-sided stamps at the boundaries. The order of these FD stamps has to remain low since the application of high-order single-sided FD schemes gives rise to numerical instability. In Ref. 23 an algorithm is described that obtains the weights for FD schemes of arbitrary order and skewness.
5 APPLICATION TO A PLANAR WAVE PROPAGATION PROBLEM
In order to illustrate and evaluate the presented method, we will apply it to a nonlinear one-dimensional space-time domain wave propagation problem in a homogeneous, lossless medium. The set of equations described by Eqs. (1)–(6) reduces to two scalar equations, given as
∂xp + ρ0∂tv = f − ρ0v∂xv − ρ0κ0pDtv, (18) ∂xv + κ0∂tp = q − κ0v∂xp − κ20(1 − 2β)pDtp. (19) We define the source to be a volume impulse source atx = 0, excited with a harmonic signal with a gaussian envelope, i.e. S = [ 0 q(x, t) ]T with
Here,Q0, td, twandf0are the source amplitude, pulse delay, pulse width and center frequency, respectively. For the linear propagation problem, the pressure level at the source is related to Q0viap = 12ρ0c0Q0.
It can be shown that the Green’s function for this problem is G(x, t) = κ0∂t −∂x −∂x ρ0∂t g(x, t), g(x, t) = c0 2H t − |x| c0 , (21)
whereg(x, t) is the one-dimensional free space Green’s function24, andH(t) is the Heaviside step function.
We would like to limit the spectral content ofq(x, t) in the spatial and temporal dimensions to frequencies belowk = π/∆x and ω = π/∆t, respectively. Therefore, the source is filtered in thex-dimension using an ideal filter, resulting in
q′(x, t) = sinc(πx/∆x) q(t). (22)
In the t-dimension, the source function is not filtered, as it is already sufficiently bandlimited for signals with a pulse length of six or more pulses and step sizes ∆t = ∆x/c0 < 4/f0. Subsequently,q(x, t) is windowed and discretized, giving
q(m∆x, n∆t) = q(n∆t)/∆x, m = 0
0, otherwise. (23)
Regarding the regularization of G, we will focus on the function g(x, t). Evaluation of Eq. (15) in both x and t gives an expression for which no closed form is available. Instead, we use the alternative procedure outlined in Section 3 to obtain the DFT coefficients ofg(x, t). For this specific Green’s function, the windowing operation in x and t can be expressed and simplified as
gX,T(x, t) = g(x, t)[H(x + X) − H(x − X)][H(t + T ) − H(t − T )]
= g(x, t)[H(t) − H(t − T )], (24)
whereT = X/c0 = (N − 12)∆t. We study the windowing operation in the spatial Fourier and temporal Laplace domain. The transform ofg(x, t) is obtained as
ˆ g(k, s) = Z ∞ −∞ Z ∞ −∞ g(x, t) exp(−ikx − st)dxdt = 1 k2+ c−2 0 s2 , (25)
wherek is the spatial Fourier parameter and s = s0+ iω is the temporal Laplace parameter. If s0 ↓ 0 the temporal Fourier transform of g(x, t) is obtained, with the singularities s = ±ic0|k|. The transform ofgX,T can be written as
Re(s′) Im(s′) jc0|k| −jc0|k| s s0
Figure 1: Poles and integration contour used for the evaluation of Eq. (26).
where 0 < sc < s0. We evaluate Eq. (26) by means of complex contour integration. The integration path is closed by a semicircle in the left halfspace, see Fig. 1. As the integrand is O(s′−3) for |s′| → ∞, the integral over the semicircle goes to zero when the radius goes to infinity, and by the residues of the poles Eq. (26) results in
ˆ gX,T(k, s) = ˆg(k, s) 1 − s exp(−sT )sin(c0kT ) c0k − exp(−sT ) cos(c0kT ) . (27)
The third step in the alternative procedure, i.e. the regularization ofgˆX,T, is a trivial task in the transform domain, and we discretize the result with step sizes ∆k = π/X, ∆ω = π/T . To avoid the singular points in ˆg, we keep a small positive s0 in the transform parameter. In the FFT’s, this is accounted for by multiplying the function withexp(−s0t) before the forward FFT and byexp(s0t) after the inverse FFT. These terms are kept close to 1 by setting s0 = 0.1/T .
6 RESULTS
The algorithm is validated with an implicit solution of the lossless Burgers’ equation17. For the medium we take water, with parametersρ0 = 998 kg m−3, c0 = 1480 m s−1, κ0 = 1/ρ0c20 and β = 3.5. For the source excitation we use parameters that represent typical values for diagnostic ultrasound,f0 = 1 MHz, Q0 = 1 s−1,td = 6/f0andtw = 1.5/f0, which results in a peak source pressure of740 kPa.
0 1 2 3 4 5 6 7 −80 −60 −40 −20 0 f(MHz) dB 0 0.02 0.04 0.06 0.08 0.1 −80 −60 −40 −20 0 x(m) dB n=1 n=2 n=3 n=4 n=5 n=6 n=7 F0 H2 H3 H4 H5
Figure 2: Left: Power spectrum of the signal pressure p at x = 0.1 m as obtained by the Burgers’ equation solution (thick line) and by our algorithm aftern = 1 to n = 7 iterations with ∆t = ∆x/c0 = 1/14f0. We
observe that each iteration gives a better estimate for increasingly higher harmonics. Right: development of the fundamental (F0) and four of the higher harmonics (H2 to H5) during propagation.
each iteration gives a better estimate for increasingly higher harmonics. Table 1 shows the relative square errors for the peak values of the fundamental and several harmonic components atx = 0.1 m for n = 1 to n = 7 iterations. From the boldface entries in the table we observe that in this situation we needn+2 iterations to reproduce the n-th harmonic within a 2% relative error.
Square relative error (%)
n F0 H2 H3 H4 H5 1 0.184 97.0 100 100 100 2 0.184 6.38 86.3 100 100 3 0.02 6.37 12.2 65.8 91.2 4 0.06 0.11 8.89 17.7 34.1 5 0.00 0.24 0.57 8.31 18.9 6 0.00 0.04 0.42 0.88 7.49 7 0.00 0.06 0.11 0.17 1.89
Table 1: Relative square error against number of iterationsn for the fundamental (F0) and four of the higher harmonics (H2 to H5) atx = 0.1 m .
In Fig. 3 we investigate the stability of the algorithm close to the theoretical shock formation distance x = 2c¯ 2
2 4 6 8 10 −3 −2 −1 0 1 2 3 t (µs) p (MPa) 0 1 2 3 4 5 6 7 −80 −60 −40 −20 0 f (MHz) dB
Figure 3: Waveform (left) and power spectrum (right) of the signal pressurep close to shock formation distance
(x = 0.98¯x) as obtained by the Burgers’ equation solution (dashed), and by our algorithm (solid), with fmax= 7f0
andn = 8 iterations. We only observe a small deviation in the waveform at the steeper slopes and overestimates
in the highest frequency components that are accounted for.
figure, the waveform and power spectrum are plotted at x = 0.98¯x. It appears that as long as X < ¯x the iterative scheme is stable and results in a waveform that is as accurate as is permitted by the number of included harmonics. The highest harmonic components are slightly overestimated due to the influence of the now significant harmonic components at frequencies larger thanfmax. This means that as we come close tox, we need to decrease the step size to¯ obtain accurate harmonics. When we have ∆t = ∆x/c0 = 1/18f0, then the square relative error for the fifth harmonic (H5) after 10 iterations is 1.48%.
7 DISCUSSION AND CONCLUSION
In this paper, we have introduced a nonlinear contrast source method based on the Neumann iterative method combined with a Green’s function approach as an algorithm to compute non-linear wave propagation. We have shown that with the Green’s function approach the problem can be efficiently discretized. The method was applied on a planar nonlinear wave propagation problem. The results already showed to be accurate after a small number of iterations, and the stability of the iterative scheme near shock wave formation distance was also shown to be good. The experiments that were performed for this paper were done by a Matlab program on a standard desktop PC and required computation times in the order of minutes.
all spatial and temporal dimensions, which will be vital for the computation of 3D, large-scale, full-wave, nonlinear wave propagation problems in the time domain.
Two drawbacks of the method are that firstly the method requires the storage of the field in all dimensions, and secondly that the used operators, a multi-dimensional convolution sum and a number of high-order finite differences, are expensive in terms of computation time. Both issues are countered by the efficient discretization. In the further development of the method, these issues will become our main challenge.
In the future, we will extend the method to three dimensions and to media with attenuation. To enable the computation of 3D time domain nonlinear wave problems, both memory and com-putation time issues seem to require the employment of high performance parallel computing systems.
ACKNOWLEDGEMENTS
The authors would like to acknowledge the financial support of the ICT Delft Research Centre, Delft, The Netherlands.
REFERENCES
[1] J.A. Jensen and N.B. Svendsen, ”Calculation of pressure fields from arbitrarily shaped, apodized and excited ultrasound transducers”, IEEE Trans. Ultrason. Ferroelec. Freq. Control, 39(2), 262-267, (1992).
[2] M. Tabei, T.D. Mast and R.C. Waag, ”A k-space method for coupled first-order acoustic propagation equa-tions”, J. Acoust. Soc. Am., 111(1), 53-63, (2002).
[3] G. Muir and E.L. Carstensen, ”Prediction of nonlinear acoustic effects at biomedical frequencies and inten-sities”, Ultras. Med. Biol., 6, 341-344, (1980).
[4] V.F. Humphrey, ”Nonlinear propagation in ultrasonic field: measurements, modelling and harmonic imag-ing”, Ultrasonics, 38, 267-272, (2000).
[5] R.P. Ward, K.A. Collins et al., ”Harmonic imaging for endocardial visualization and myocardial contrast echocardiography during transesophageal echocardiography”, J. Am. Soc. Echocardiogr., 17(1), 10-14, (2004).
[6] K.S. Sodhi, R. Sidhu et al., ”Role of tissue harmonic imaging in focal hepatic lesions: comparison with conventional sonography”, J. Gastroenterol. Hepatol., 20(10), 1488-93, (2005).
[7] A. Bouakaz and N. de Jong, ”Native tissue imaging at superharmonic frequencies”, IEEE Trans. Ultrason. Ferroelec. Freq. Control 50(5), 496-506, (2003).
[8] T.J. Christopher and K.J. Parker, ”New approaches to nonlinear diffractive field propagation”, J. Acoust. Soc. Am., 90(1), 488-499, (1991).
[9] Y.-S. Lee and M.F. Hamilton, ”Time-domain modeling of pulsed finite-amplitude sound beams”, J. Acoust. Soc. Am., 97(2), 906-917, (1995).
[10] J. Hoffelner, H. Landes, M. Kaltenbacher and R. Lerch, ”Finite Element simulation of nonlinear wave propagation in thermoviscous fluids including dissipation”, IEEE Trans. Ultrason. Ferroelec. Freq. Control, 48(3), 779-786, (2001).
[11] I.M. Hallaj and R.O. Cleveland, ”FDTD simulation of finite-amplitude pressure and temperature fields for biomedical ultrasound”, J. Acoust. Soc. Am., 105(5), L7-L12, (1999).
[12] S. Ginter, ”Numerical simulation of ultrasound-thermotherapy combining nonlinear wave propagation with broadband soft-tissue absorption”, Ultrasonics, 37, 693-696, (2000).
[13] Wojcik, G., J. Mould, et al., ”Nonlinear modeling of therapeutic ultrasound”, Proc. 1995 IEEE Ultrasonics Symposium, 1617-1621, (1995).
[15] X. Yang and R.O. Cleveland, ”Time domain simulation of nonlinear acoustic beams generated by rectangular pistons with application to harmonic imaging”, J. Acoust. Soc. Am., 117(1), 113-123, (2005).
[16] R.J. Zemp, J. Tavakkoli and R.S.C. Cobbold, ”Modeling of nonlinear ultrasound propagation in tissue from array transducers”, J. Acoust. Soc. Am., 113, 139-152, (2003).
[17] M.F. Hamilton and D.T. Blackstock, (Eds.), Nonlinear Acoustics, Ac. Press, San Diego, (1998). [18] A.T. De Hoop, Handbook of radiation and scattering of waves, Ac. Press, San Diego, (1995).
[19] J. Huijssen and M.D. Verweij, ”Nonlinear constitutive equations derived for fluids obeying an ideal gas, a Tait-Kirkwood or a B/A type equation of state”, Proc. of 17th ISNA, in press, (2006).
[20] M.D. Verweij and J. Huijssen, ”Nonlinear and dissipative constitutive equations for coupled first-order acoustic field equations that are consistent with the generalized Westervelt equation”, Proc. of 17th ISNA, in press, (2006).
[21] R. Tolimieri, M. An and C. Lu, Algorithms for discrete fourier transform and convolution, Springer-Verlag, New York, (1989).
[22] W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, Numerical Recipes in Fortran 77, 2nd ed., Cambridge University Press, New York, (1996).
