A method for approximating high frequency sound radiation—The plane projection Rayleigh integral Open Access

The simulation of sound radiated by vibrating arbitrary objects, suspended in air or other fluids such as water, is a common task in technical acoustics. In order to enable rapid and cost-effective development, it is crucial to predict the sound radiated by vibrating surfaces at an early stage. High frequencies, resulting in small wavelengths compared to object sizes or radiation distances, pose challenges, leading to significant computational effort. The typical field of applications includes airborne ultrasonics or medical ultrasound diagnostics. However, these challenges also present opportunities to use approximations that reduce the problem's complexity. We propose a method for simulating sound radiation at high frequencies that is low in implementational complexity, efficient in computational effort, and proven to provide high accuracy for our test case of a radiating sphere at high frequencies. This method is referred to as the plane projection Rayleigh integral (PPRI).

The most commonly used methods for simulating sound radiation are the boundary element method (BEM),¹ the finite element method (FEM),² and the source superposition method.³ Each of these methods has its own advantages and disadvantages. A common challenge among them is the rapid increase in computational effort at high frequencies. This issue arises due to the higher discretization density needed to achieve the desired accuracy. Efforts are being made to address these challenges by improving the methods themselves, such as through enhanced versions of BEM.^4–10 Additionally, alternative methods are being explored that use approximations at high frequencies. These methods simplify the underlying acoustic equations so that simulation requires only the numerical computation of a surface integral. This category of approximation methods includes the visible element Rayleigh integral (VERI) and the high frequency boundary element method (HFBEM).¹¹ 

The proposed PPRI also falls into this category. The methodology consists of two steps (Fig. 1). First, the surface velocity distribution $v→¯n,s$ of a harmonically vibrating three-dimensional object is projected onto a virtual plane leading to $v→¯n,p$⁠. Second, the Rayleigh integral is used to calculate the sound radiated from this projection plane as an approximation of the sound radiated by the three-dimensional object. It is important not to confuse the method with the projection plane method, which addresses singularities in boundary integrals by analytically transforming the integral onto a tangential plane.¹² 

Another method that demonstrates how a virtual plane can represent a more complex three-dimensional surface by approximating sound radiation using the Rayleigh integral is the boundary element Rayleigh integral method (BERIM).¹³ The BERIM is successfully used to simulate radiated sound from an open cavity by virtually enclosing the cavity with a plane. The interior problem is solved using BEM, and the exterior problem is then addressed by applying the Rayleigh integral to the plane covering the opening.

The objective of this paper is to provide an initial assessment of the method in comparison to the VERI and HFBEM, to ascertain its operational boundaries, and to determine whether it combines the attributes of these methods in practical applications as intended. The methods are evaluated by calculating the sound field of a breathing and oscillating sphere across a range of radius and frequency combinations at various field points, primarily varying the distance from the surface. Analytical solutions for sound radiation in both cases are used as a reference for error calculations. The amplitude and phase errors of the sound pressure are considered separately. Median error values over the spatial field point distribution and the associated standard deviation are used as metrics for general precision. Additionally, a custom metric is introduced to assess the dependence on the distance of the field points and convergence of error values.

First, the underlying acoustic equations are introduced, followed by the high frequency approximations used in the VERI and HFBEM. The PPRI is introduced thereafter. Alongside the reference solution, the implementation of the methods for the radiating sphere is presented. The parameterization of the test case is outlined, and the metrics for evaluation are introduced. Finally, the numerical results are presented and discussed.

Although the method is evaluated for airborne sound, the method is applicable to all fluid domains and the results are generally valid.

The calculation requires knowledge of both the surface velocity distribution $v→¯(x→s)$ and the sound pressure $p¯(x→s)$ immediately above the surface. In many cases, the surface velocity distribution is known. However, determining the surface sound pressure is a time-consuming process. The primary factor influencing the calculation time of the BEM is solving the integral equation [Eq. (1)] for the surface sound pressure.¹ 

The equations from Sec. II represent the analytical solutions for their respective problems.^16–19 Additionally, approximation methods to reduce complexity and computational effort are based on these integral equations. These approximations require the consideration of high frequency sound radiation and are necessary due to the significant increase in computational effort at those frequencies, as a consequence of higher discretization densities.

The initial and most straightforward approximation is to utilize the Rayleigh integral [Eq. (2)] itself to calculate the sound field produced by a non-planar object.^20–23 In many cases, objects have surface areas that are nearly planar. Even curved areas are approximated by the Rayleigh integral. This approximation is more effective at higher frequencies, due to the fact that the areas appear more planar in relation to the shorter wavelengths, i.e., when the radius of curvature is large compared to the wavelength. The simulation of the sound field is consequently a weighted sum of the normal surface velocity at discrete surface points.

When attempting to approximate sound pressure radiated by objects in all spatial directions, the Rayleigh integral alone is inadequate. This is due to the fact that the integral itself does not distinguish between surface areas facing towards or away from the field point considered. In order to address this issue, the VERI is proposed in the literature as a solution.¹¹ The idea is to include only those surface areas that are “visible” to the field point in the Rayleigh integral. The determination of visible surface points is accomplished through the scalar product of the surface normal, represented by the vector $n→,$ and the direction from the surface point to the field point, represented by the vector $d→$ (Fig. 3). For a visible point, the scalar product will be positive; for an invisible point, it will be negative. This way of determining visible surface points generally only applies to convex surfaces. The border is formed by perpendicular vectors and a scalar product of zero. Consequently, the relevant surface area for the sound pressure approximation is reduced significantly and changes with each field point considered. This is shown by the purple marked area denoted as $Sv$ in Fig. 3.

The existing literature indicates that additional requirements must be met to achieve optimal approximations using the HFBEM.²⁷ Primarily, the dimensions of the surface areas undergoing vibration in phase must exceed the wavelength. So, the ratio of the spatial change in surface velocity to the wavelength limits the validity of the method. Second, the radiator must possess overall convexity to prevent backscattering of radiated sound by the object itself.

In contrast to the VERI, the surface points are not weighted binary (visible or not) but continuously by the scalar product $n→⋅d→$⁠. The weights range from zero, indicating a completely averted surface point, to 2, indicating a surface point completely facing the field point considered. This is due to the subtraction of the scalar product from 1 in the integral equation [Eq. (5)]. The points indexed 1 and 2 are assigned weights of 1 and zero by the VERI and weights of 1.93 and 0.24 by the HFBEM (Fig. 3).

For the derivation of the method, consider a baffled circular piston transducer. The Rayleigh integral [Eq. (2)] is the analytical solution for this case. The radius of surface curvature, being infinite for plane surfaces, is reduced gradually. The sound pressure radiated by the slightly curved surface is still approximated with high accuracy by applying the Rayleigh integral to the original plane, as long as the maximum deflection due to the curvature is much smaller than the wavelength of radiated sound. This approach effectively neglects the occurring phase shift. If the radius is decreased further, the phase shift can no longer be disregarded. For approximating the radiated sound pressure by using the Rayleigh integral on the original plane, the phase shift needs to be included, similarly to the example of approximating the sound field generated by a conical loudspeaker.¹⁵ The surface curvature also has an influence on the direction of sound radiation. Again, this effect is negligible for sufficiently large radii of surface curvature but needs to be considered when the curvature radius is reduced. Furthermore, the directivity of sound radiation increases with radiation frequency. Thus, the surface points are weighted according to the direction of sound radiation in addition to the phase shift, increasing the accuracy of this approach at high frequencies. These steps lead to the proposed PPRI. Because the development is not restricted to constant radii of surface curvature, the method is applicable for approximating the sound radiation of arbitrary objects at high frequencies.

A crucial point of the method is that the position and orientation of the projection plane varies with each field point considered for the sound pressure calculation. The normal direction of the plane coincides with the direction from the object's center of mass or volume to the field point. Many options for the positioning along this axis may be considered, with close proximity to the real surface being favored for smaller errors. For the sake of getting a first general evaluation of the method, only one option is used. The plane is positioned tangentially to the surface at the foremost surface point in the normal direction. Therefore, the original and projected surface velocities are the same at the point of contact.

In some respects, the PPRI methodology integrates elements of both the HFBEM and VERI approaches. In particular, some surface points are completely disregarded in the sound pressure approximation, in line with the concept of visible elements in the VERI. The Rayleigh integral is also used to approximate the sound pressure of a non-planar radiator. In addition, other surface points are continuously weighted according to the scalar product between the normal directions of the projection plane and the original surface, a characteristic feature of the HFBEM, where the scalar product is formed between the normal direction of the surface and the direction to the field point considered.

The PPRI therefore aims to combine the lower computational effort of the VERI with the higher precision of the HFBEM. This is achieved by ignoring surface points that have little influence on the sound pressure at the field point but taking into account the other points with individual weights. Considering arbitrary surfaces and surface velocity distributions, the projection can be computed using established methods from computer graphics, when a virtual discretized model of the object is generated. This promises to be fast and efficient.²⁸ Together with the smaller and simpler surface for numerical integration with the Rayleigh integral, this method is expected to have advantages in computational time compared to the HFBEM and the VERI.

The harmonically radiating sphere with an angular frequency of $ω$ is used as a three-dimensional test object for the evaluation of the PPRI in comparison to the VERI and the HFBEM. It is the most basic sound radiator with a full set of analytical solutions. The analytical solution of the sound radiation for specific surface velocity distributions on the sphere is compared to all three methods under test. The Rayleigh integral itself is the analytical solution for planar radiators, and in combination with rounded surfaces that resemble or are close to a spherical shape, many relevant geometries are covered in some way by this investigation. In accordance with the previous derivation of the PPRI, the use of the sphere represents a worst case for the evaluation of the method. Here, reducing the radius of surface curvature to the radius of the piston transducer will result in the sphere, when keeping the curvature radius constant for all surface points. Thus, the errors calculated in this evaluation indicate an upper bound for the method in this consideration.

In order to implement the methods for the reference objects, the respective surface integrals have to be transformed into spherical coordinates and then evaluated numerically over the surface of the sphere. Even in the case of the VERI, the integral is evaluated over the entire surface of the sphere. Here, the selection of the visible elements can be done in the integrated function by multiplying the integrand by the step function of the scalar product between the normal surface direction and the direction of the surface point to the field point. Consequently, the integrand is equal to zero for all surface points that are not visible to the field point, resulting in a negative scalar product.

The matlab function integral2 is used for numerically evaluating the surface integrals. The function uses both Gauß and Kronrod quadratures to assess the integration error section-wise and adapt the discretization of the section accordingly to yield efficient distribution of evaluation points for high precision and low computational effort.^30,31

The realization of the PPRI is presented in more detail here. Instead of directly applying an integral to the surface of the sphere, the method adds the step of first projecting the normal surface velocity from the sphere $v¯n,s(R,θ,φ)$ onto a virtual tangential projection plane, $v¯n,p(a,β),$ and then applying the Rayleigh integral there. Generally, the projection mapping surface points from the original object to corresponding points on the virtual plane can be split into three parts: first, finding corresponding points on both surfaces; second, projecting the original surface velocity vector onto the normal direction of the virtual plane by using the scalar product; third, shifting the phase of the velocity according to the distance the sound wave theoretically travels from the object to the virtual plane. These three steps are applicable to all arbitrary discretized objects.

By using the sphere as a reference object with surface velocity distributions defined by a continuous function, the projection can be done continuously as well. That means a function for the normal surface velocity on the virtual plane, depending on its local coordinates, is derived here. Thereafter, the Rayleigh integral can be evaluated numerically on the virtual plane.

In order to assess the performance and limitations of the method, the radius and frequency parameters are varied, as shown in Table I. Radii started from 1 cm and the frequency in the hearing range at 10 kHz, where the main requirement for the HFBEM is not met [Eq. (4)]: $R/λ≈0.3$ for $c=343 m/s,$ which is significantly smaller than a factor of 10 for the condition of R being much larger than $λ$⁠. The largest radius used is 70 cm, which in combination with the highest used frequency of 80 kHz strongly satisfies the HFBEM's main requirement of $R/λ≈163.$

For all combinations of these parameters, the complex sound pressure is calculated for discrete field points using all three methods. The results are compared with the analytical reference solutions. Different positions and distances from the spherical surface are used to evaluate the dependence of the method on these variables. Due to the different symmetries in the sound fields of the breathing and oscillating spheres, different field point distributions are used for these two cases.

In the case of the breathing sphere, a single axis is sufficient to evaluate the entire field, due to the point symmetry of the sound field. This axis extends from the near field of the sphere to a maximum distance of 5.5 m from its center. Along this axis, 6000 field points are equally spaced to allow adequate evaluation of the sound pressure phase, with at least four points per wavelength.

For the oscillating sphere, a quarter circle area in a plane coinciding with the axis of oscillation is required to evaluate the entire sound field. Here, 900 field points are evenly spaced in the area of the quarter circle to evaluate the sound pressure amplitude. The points are positioned on circular arcs with increasing radii from above the surface up to 5.5 m in constant steps. They are equidistant on the circular arcs of all radii. Again, the sound pressure amplitude is also evaluated on an axis discretized by 6000 field points. This time two axes are used to evaluate the orientation aspect. One axis is in line with the oscillation axis of the sphere, and the other is at a 45° angle to it.

The error of the approximation methods with respect to the analytical reference solution is the basis for the evaluation in this paper. For the sound pressure amplitude, the relative error is calculated, and for the sound pressure phase, the absolute error is used. Three metrics are used for both phase and amplitude errors to reduce the amount of resulting data for evaluation.

The sound pressure values are used for the evaluation. Compared to a possible evaluation with the quantity sound power, this is the more demanding case. By integrating sound pressure values to get the radiated sound power, random numerical errors are reduced due to averaging effects and sound pressure errors in different field points might even cancel each other out. Therefore, the method's accuracy only increases when considering sound power calculations.

The overall precision of the methods is quantified by the median error of the spatial field point distribution. The median is used instead of the mean value because large error values in field point regions where the analytical reference solution is close to zero shift a medium error value disproportionately, leading to false conclusions. In addition, the standard deviation from this median is used as a second metric to qualify and quantify the variation of error values across the spatial distribution of field points.

As another result, the PPRI is examined in a more detailed and frequency-orientated manner. For a single test case of the breathing sphere with a constant radius, the method is compared to a version of itself without the use of the amplitude factor [Eq. (15)].

A way to visualize the numerical results is to plot the sound pressure amplitude over the radial distance from the sphere's center (Fig. 6). It illustrates the behavior of the methods, which is further quantified by the metrics used to evaluate the performances.

The VERI method results in a particularly different curve compared to the reference solution, HFBEM, and PPRI [Fig. 6(a)]. In addition to the exponential decrease in sound pressure amplitude, an oscillatory component is present. This phenomenon may be explained by the fact that the area designated as “visible” changes on the order of magnitude of the wavelength, depending on the distance of the evaluated field point from the surface of the sphere. This variation significantly impacts the interference pattern of the radiated sound from the spherical cap. There will be constructive interference at one field point and destructive interference at another.

The other two methods match the reference solution better, as seen in the enlarged detail in Fig. 6(a). The PPRI appears to be even closer to the reference solution than the HFBEM. This can be further examined through the error curves of these methods [Fig. 6(b)]. Both methods increase their prediction accuracy with increasing distance from the surface. Similar to the VERI, the PPRI shows an oscillatory component, but of significantly smaller magnitude. This may also be explained by the interference of sound waves due to the limited projection plane. These aspects are particularly quantified by the introduced variation factor V.

Our results showed virtually no difference between the error metrics of the two axes used to assess the field of the oscillating sphere. Therefore, the performance of the methods does not depend on the positioning of the field point at a given distance around the sphere. In this case, only the axis of oscillation will be used for evaluation.

The median of the relative amplitude and absolute phase error combined with the standard deviation are the metrics used to evaluate the overall precision of the methods (Fig. 7). The previously introduced lower limit for the Helmholtz number and a 1% error threshold are used to put the results into perspective. For the phase, the threshold is applied to $2π$⁠. The ordinates of the amplitude and phase plots are scaled accordingly to allow easy comparison.

The VERI shows large error values and comparatively large standard deviations [Figs. 7(a) and 7(b)]. For the amplitude, the maximum error values almost reach 100% and hardly fall below the 1% threshold. The method is able to approximate the sound pressure phase with lower error values. There is no clear trend for the dependence of the median values on the Helmholtz number, but the standard deviation seems to increase with increasing Helmholtz number. A large standard deviation is in accordance with the oscillatory component of the trend shown earlier [Fig. 6(a)]. Nevertheless, the spread of the median values around a constant is large. Furthermore, the method shows a clear dependence on the vibration mode of the surface. In fact, the accuracy in both amplitude and phase increases significantly when the oscillating sphere is considered rather than the breathing.

In contrast, there is a clear monotone decrease for the HFBEM error values [Figs. 7(c) and 7(d)]. This was to be expected given the main requirement for the method's deployment, which is more strongly fulfilled at larger Helmholtz numbers. The spread around a linear decline in this double logarithmic representation is quite large, more so for the phase error than for the amplitude error. The data indicate that the upper limit of the spread is linearly decreasing. It may be the case that the performance of the method is more dependent on the radius of the sphere than on the vibration frequency. In that case, the method will compute lower error values across the range of Helmholtz numbers if the largest sphere with a radius of 70 cm is considered at the different frequencies. A preliminary examination of the data before compressing them by means of the Helmholtz number suggests that this is a plausible hypothesis. Again, the amplitude error values are slightly larger than the phase error values. In case of the amplitude error, the values fall below the 1% threshold at the theoretically derived lower limit of the Helmholtz number of 20. The method shows a negligible difference between the curves for the breathing mode and oscillating mode. The standard deviation is so small compared to the error values that the error bars are not visible, which concludes a rather constant error progression along the radial distance from the sphere surface.

The PPRI shows the same decreasing trend in error values with increasing Helmholtz numbers [Figs. 7(e) and 7(f)]. The spread around a linear decline is significantly smaller compared to the HFBEM, and the error values decrease at a larger rate. This time the sound pressure amplitude errors are slightly lower than the phase error values. In this case, a lower Helmholtz number limit of about 5 can be readout visually. A small dependence of the method's performance on the vibration mode of the object can be observed for the sound pressure amplitude. The error values decrease at a larger rate for the oscillating sphere than the breathing mode. The standard deviation, on the other hand, is much larger compared to the HFBEM, but still significantly below those obtained by the VERI. The standard deviation increases with increasing Helmholtz number. This can also be attributed to a kind of oscillatory component in the error progressions, as seen above [Fig. 6(b)], originating from interference effects in the area closer to the projection plane, therefore depending on the ratios of frequency, radius, and distance.

The variation factor can be evaluated in a similar way (Fig. 8). Here, all plots are scaled equally and the range of values is given, allowing a straightforward comparison between methods in this respect. The variation factor further quantifies the aspects already considered by the standard deviation and provides a statement on the convergence towards a constant error value.

The VERI results show the highest values of the three methods [Figs. 8(a) and 8(b)], which increase as the Helmholtz number increases. This shows that the oscillatory component of the course increases at larger Helmholtz numbers and converges to a constant value at greater distances. This supports the given explanation for the behavior of the VERI, since the “visible” area changes at a decreasing rate as the distance of the field points increases. This rate of change depends mainly on the ratio of the distance to the radius of the sphere and therefore on the Helmholtz number, since the range of distances considered is fixed. Again, the VERI yields lower values when considering an oscillating sphere rather than a breathing sphere.

The HFBEM also shows this increasing tendency, but with the smallest value range of all methods [Figs. 8(c) and 8(d)]. Therefore, the error curves of the method can be considered as rather constant. In contrast to the VERI, this method yields smaller values for the breathing mode than for the oscillating mode.

The PPRI ranks between the other methods, but the values decrease slightly with increasing Helmholtz number [Figs. 8(e) and 8(f)]. The range overlaps slightly with the VERI's range at smaller Helmholtz numbers and continues to decrease, whereas the VERI values start to increase strongly. This means that the small oscillatory component [Fig. 6(b)] seen above decreases with increasing Helmholtz number and converges at smaller distances to the sphere surface. Similar to the VERI, the method results in smaller values for the oscillatory mode.

In all cases, the amplitude and phase values calculated behave similarly. Therefore, both aspects of the sound pressure show the same dependence on the distance of the field point from the object surface.

Almost all results show a better performance (precision and dependence on distance) of the methods when applied to an oscillating sphere. This gives confidence that the breathing sphere was not a low-demand reference object for the evaluation of the PPRI's performance.

Another result is a more frequency-oriented evaluation of the PPRI and the consideration of computational time (Fig. 9). The previous results yield a general understanding of the method's performance, depending on radius, frequency, and distance. For a more detailed view, the breathing sphere is investigated with a constant radius of 0.1715 m but varying frequency. The error of the sound pressure amplitude is calculated in a single field point [Fig. 9(a)]. This radius is chosen because of its relation to the speed of sound. At a frequency of 1 kHz, the wavelength equals the diameter of the sphere, resulting in a Helmoltz number of 1. Thus, the Helmholtz number is computed by dividing the frequency by 1 kHz, enabling straightforward consideration of frequencies and Helmholtz numbers at the same time. In addition to the PPRI, also a version without the amplitude factor is considered, where only the phase shift is taken into account.

The results show that an upper bound for the error of the PPRI in both versions exists. Above a certain frequency, this upper bound decreases gradually with increasing Helmholtz number and frequency. This is in accordance with the previous observation that the PPRI is only applicable for higher frequencies and improves its accuracy with increasing frequency. The necessity of the amplitude factor becomes evident in this consideration, because it accelerates the convergence of the error to zero considerably. Without the amplitude factor, the 1% threshold is undercut only in the mid-ultrasonic frequency range (about 50 kHz), beyond the derived Helmholtz number threshold of 20 for the HFBEM. With the amplitude factor, the error falls below 1% at frequencies above 5 kHz in this case. This Helmholtz number threshold of 5 for the PPRI is also derived from the previous, more general evaluations.

The focus of the implementation of the methods is on convenience rather than computational time. For example, the continuous projection on the level of functions in the PPRI is beneficial for the computational time, so the time increases when implementing the method for arbitrary discretized objects. Nevertheless, this consideration gives a first impression of what can be expected in terms of computational effort. The computational effort of the methods is compared by using the time needed for calculating the sound pressure at 80 kHz in a single field point [Fig. 9(b)]. The data were computed in matlab on an Apple M1 chip with 8 Gbytes of RAM. The HFBEM proves to be the most time-consuming method, with about 0.104 s in this specific case, followed by the VERI, which is about twice as fast, with 0.06 s, and PPRI, which is 5 times as fast, with about 0.019 s.

The computational effort increases linearly with the number of nodes in all methods. The advantage of the PPRI over the HFBEM in this regard is that many fewer nodes on the considered surfaces are required to achieve the same discretization density. Furthermore, the numerically integrated functions are less complex and require fewer computations to evaluate.

An approximation method for the simulation of high frequency sound radiation was successfully introduced. The approach offers benefits compared to established methods in both computational effort and precision.

In terms of overall precision, the PPRI results in significantly smaller error values than the VERI. The error values are of similar magnitude to those computed by the HFBEM and show a declining trend as well. However, the PPRI demonstrates a more rapid and predictable decline, quantifiable as a lower Helmholtz number limit. When considering a 1% error threshold, the theoretically derived Helmholtz number limit for the HFBEM of 20 is recognizable. For the PPRI, a limit as low as 5 can be determined. Consequently, the method is applicable to even smaller radii of curvature or performs better at the same radii as the HFBEM.

With regard to the dependency of the method's performance on the distance of the field points from the surface of the object, the PPRI is inferior to the HFBEM. However, the dependence both observed [Fig. 6(b)] and quantified (Fig. 8) is so minor that it does not constitute a significant issue. This is in contrast to the pronounced oscillations observed for the VERI [Fig. 6(a)].

In terms of computational efficiency, the PPRI exhibits notable advantages, demonstrating a significant increase in speed compared to the HFBEM in the calculations conducted for this assessment. This suggests that the method is achieving its intended objectives, combining the lower computational effort of the VERI with the higher precision of the HFBEM, even surpassing both methods in those regards.

Though the method is only evaluated with airborne sound, its application is not restricted to that case and can be applied to vibrating objects in other fluids as well.

Following this initial assessment of the PPRI through the utilization of an analytical reference solution, further research requirements can be addressed. Presently, the application of the PPRI is presumed to be limited to convex surfaces. Surfaces with great complexity and overlapping areas due to concavity are probably an issue. As the literature suggests, these are applications where the HFBEM struggles as well.²⁷ Further research will be concerned with those more complex and asymmetric geometries. In this context, it is appropriate to compare the methods to the regular BEM as a reference and to examine the computational complexity and effort more closely. Furthermore, optimization of the positioning of the projection plane will be a component of future research.

This work has been partially funded by the technology transfer program “TTP Leichtbau” of the German Federal Ministry for Economic Affairs and Climate Action (Bundesministerium für Wirtschaft und Klimaschutz; Reference No. 03LB3029).

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.