Fundamental method and its hardware implementation for the generic prediction and analysis of multiple scattering of waves in particulate composites

[0001] Attached hereto there is a CD-R, hereby incorporated by reference, as an embodiment for the present invention. By following industry-standard installation procedures in a computer operating under a member of the MS-Windows family of operating systems and a Pentium II or higher equivalent processor, full disclosure for the current invention is enhanced, and its operability can be confirmed with actual experimental data provided in the disk. A commercial product based on this invention, to be called SCATTERER will be marketed in the near future. The copyright owner has no objection to legal use and reproduction of this patent document and attached CD-R, as they appear in the Patent and Trademark Office records, reserving otherwise all rights.

[0002] This invention relates to a method and its hardware implementation for the fundamental prediction, analysis, and parametric studies of the interaction of multiply scattered waves of any nature (acoustic, electromagnetic, elastic) with particulate composites (solid or fluid particulate phases included in a solid or fluid continuous phase). Any embodiment for the present invention constitutes a powerful research and development tool in a vast variety of engineering and scientific fields dealing with wave propagation in composites. The METAMODEL Prediction Engine 21, a crucial component of this invention (FIG. 1), is suitable for integration into existing acoustic or optical particle size analyzers to thereby extending substantially their concentration ranges.

[0003] Several systems to assist in the modeling, analysis, and parametric studies of physical systems are known in the prior art and available in the market. Some of them are general mathematical tools with thousands of functions in calculus, linear algebra, ordinary differential equations, partial differential equations, Fourier analysis, dynamic systems, etc. An incomplete list includes MATLAB (The MathWorks, Inc), Maple V (Waterloo Maple Software), SCIENTIFIC WORKPLACE (TCI Software Research), MathCad (MathSoft, Inc.), and Mathematica (Wolfram Research, Inc.). These systems constitute, in a hierarchy of research and development tools, a next level up from e.g. the well-known IMSL FORTRAN/C library of scientific routines (Visual Numerics, Inc.). By becoming a little more specific in scope, they offer a higher-level language to the user, reducing the need for proficiency in computer programming and increasing the efficiency with which the scientist/engineer solves his/her problems. Some other systems are even more specific like FEMLAB (COMSOL, Inc.) or MaX-1 (A visual Electromagnetics Platform, J. Wiley) or RealTime (A Two-Dimensional Electromagnetic Field Simulator, CRC Press). The last two systems are overly specific simulators, as they cover only certain aspects of electromagnetic waves with no flexibility offered to modeling their multiple scattering in particulates. The first one, FEMLAB, using MATLAB as a software platform, offers to the user a less specific high-level language oriented towards the solution of partial differential equations using the finite element technique for problems in physics and engineering. It allows the user to create a suitable mathematical model for the physical problem at hand from available building blocks (partial differential equations, functional relationships, etc.) in an ample variety of disciplines including acoustics, fluid dynamics, structural mechanics, heat transfer, chemical processing, and electromagnetics. After a short training, the user only needs to be concerned with physics and engineering, and is required no major proficiency in computer programming and numerical intricacies. However, no pre-built modeling blocks for multiple scattering of waves in particulate composites is available either.

[0004] All these systems in the prior-art are either exceedingly general or overly specific to efficiently solve the multiply interleaved wave equations required to accurately predict the multiple scattering of waves in particulates. In order to take advantage of those existing tools, it would inevitably require an expert on both the physics and the mathematics of the problem whoin additionwould have to command the novel know-how to be disclosed in this specification, plus unacceptable computer time and human effort. A system tackling all those intricate physico-mathematical aspects of the problem, and making them transparent to the regular physicist/engineer has been clearly in need for many years. The current invention provides such a tool.

[0005] Two distinct approaches have historically developed in trying to mathematically model scattering of waves with particulates. One is microscopic (fundamental) and called the analytical wave-scattering theory; the other is macroscopic (phenomenological) and called the transport theory or radiative transfer theory (22 A. Ishimaru, Wave Propagation and Scattering in Random Media, Academic Press, 1978; 48 F. Alba et al, Ultrasound Spectroscopy: A Sound Approach to Sizing of Concentrated Particulates, Handbook on Ultrasonic and Dielectric Characterization Techniques for Suspended Particulates, The American Ceramic Society, 1998). The former starts with the wave equation, rigorously predicts the scattering of a wave with a single particle in the medium, and statistically describes the interaction between all the scattered fields around all the particles, to arrive to expected values for the attributes of the composite in terms of the attributes of its constituents. It is a mathematically rigorous fundamental approach and, as a consequence, all the multiple scattering, diffraction, and interference effects are contemplated with the broad validity and accuracy range inherent to a first-principle technique. This is the physico-mathematical background of the current invention. In the latter (macroscopic) approach, instead, the basic wave equation is not the starting point. It deals directly with the transport of energy through the particulate composite, arriving to a basic differential equation (the equation of transfer). A basic assumption is that there is no correlation between the coexisting scattered fields and, therefore, their intensities can be added directly (as opposed to vectorially adding the fields). A crucial magnitude is the composite depth (particle concentration by number times total cross-section of a single particle times path length). The solution of the equation of transfer depends upon this depth and the so-called phase function. The theory is heuristic as opposed to fundamental and, ergo, it lacks the mathematical rigor and ample validity range of a first-principle approach. As it would be expected, this approach is simpler and more tractable from the numerical point of view. It has flourished during the last 100 years, in great part, due to the lack of sufficient cost-effective computational power to implement the fundamental approach. Nowadays, duplication of computer power every year or so at lower costs, has made it possible to consider full numerical implementation of the fundamental approach using non-expensive computer hardware.

[0006] Scientific literature on the analytical multiple-scattering theory is abundant but dispersed throughout a myriad of publications. Main authors who have significantly contributed in the last fifty years are L. L. Foldy, M. Lax, P. C. Waterman and R. Truell, V. Twersky, P. Lloyd and M. V. Berry, J. G. Fikioris and P. C. Waterman, A. K. Mal and S. K. Bose, N. C. Mathur and K. C. Yeh, V. N. Bringi V. V. Varadan, V. K. Varadan, Y. Ma, L. M. Schwartz, A. J. Devaney, L. Tsang and J. A. Kong, D. D. Phanord and N. E. Berger, de Daran, F. et al, etc. Work of these researchers and others are listed in this patent. The great majority of these publications are dedicated to specific problems in specific fields of applications. A very few, however, recognizing the common underlying structural physico-mathematical core in all these fields, have attempted to formulate generic equations for the prediction of either acoustical, electromagnetic, or elastic waves in particulate composites. The publication 37 Multiple Scattering Theory for Acoustic, Electromagnetic, and Elastic Waves in Discrete Random Media by V. V. Varadan, V. K. Varadan and Y. Ma (Multiple Scattering of Waves in Random Media and Random Rough Surfaces, The Pennsylvania State University, 1985) belongs to this kind of general work and constitutes direct theoretical prior-art to this invention as far as the Prediction Engine 21 is concerned (FIG. 1). Even though the approach described in this prior-art and several others has a common theoretical background with the Prediction Engine 21, the latter constitutes a step forward in generality and accuracy, as the present specification will thoroughly demonstrate.

[0007] There is finally a third approach in dealing with the multiple-scattering phenomenon that is more limited in scope: it consists basically in attempting to suppress it by conceiving ad-hoc hardware designs and/or software schemes through which the multiple-scattered signal is either significantly attenuated in favor of the single-scattering component, or cross-correlation techniques are employed to filter one component from the other. U.S. Pat. No. 6,100,976, No. 5,956,139, and No. 4,380,817 belong to this group. This approach is mentioned here for the sake of completeness even though, since no attempt to mathematically modeling multiple-scattering phenomena is made, these patents do not truly constitute prior art to the present invention.

[0008] Two important technical fields where the present invention meets a long-term need are those of particle sizing using light scattering and acoustic spectroscopy. This is simply because, in order to indirectly measure a certain constituent attribute (e.g. particle size distribution) with a wave-based instrument, it is strictly necessary to be able to predictas a function of this constituent attributethat composite attribute that is directly being measured (e.g. scattering, attenuation, intensity). During the last 30 years or so, measuring of particle size in suspensions, emulsions, and aerosols using the scattering of light have developed to its mature state with a plethora of commercial instruments available in the market. Nonetheless, in practically all of those instruments a sine qua non condition for accurate performance is that the particle concentration should be reduced before analysis to extremely low values (typically under 0.1% v/v)when actual industrial concentrations are several orders of magnitude higher. One of the main reasons for such an operating constraint is that the mathematical model (e.g. Mie theory) internally employed by these light-scattering instruments is a single-scattering model, i.e. a model that can predict accurately the composite attributes only when the particles are well-separated between each other, so that no interaction between contiguous electromagnetic fields occurs. Multiple scattering is simply not modeled.

[0009] A light-scattering instrument for aerosols (host medium is gaseous) that can work at higher concentrations (light transmissions down to 2% as opposed to typical values of about 60% to 80% in other instruments) is the one described by Harvill, T. L. and Holve, D. J. in U.S. Pat. No. 5,619,324, issued Apr. 8, 1997, in which correction for multiple scattering phenomena is conducted through a suitable adaptation of the theory described by E. D. Hirleman in 40 Modeling of Multiple Scattering Effects in Fraunhofer Diffraction Particle Size Analysis, Particle Characterization 5, 57-65 (1988), and 41 A General Solution to the Inverse Near-Forward Scattering Particle Sizing Problem in Multiple Scattering Environments: Theory, Proceedings of the 2^nd International Congress on Optical Particle Sizing, Mar. 5-8 1990, pp. 159-168. With such a correction, this instrument extended the concentration range by a factor of 5-10 (only up to about a few percent by volume). Hirleman's model was developed for the specific case of small-angle scattering with the assumption that scattered fields add incoherently (on an intensity rather than amplitude basis). A black box approach is employed through defining a multiple-scattering redistribution matrix that represents how optical energy incident in a certain direction is redistributed by the whole composite into another direction. This matrix not only depends on the nature of medium and particles with their size distribution and concentration, but also on the geometry of the instrument as well as the optical depth of the particulate system. In plain words, the matrix can only be determined (ad hoc) once the instrument (hardware) at hand is known. Hirleman uses also a successive orders of scattering approach to express the multiple-scattering matrix in terms of the single-scattering one. Both matrices are global characterizations of the particulate system. No fundamental modeling of the interaction of the wave with a single particle and interaction between scattered waves around particles is undertaken. By further assuming that a scattering event of order n can be represented by the n-power of the single-scattering matrix applied to the incident signal and summing all contributions, Hirleman arrives to an analytical relation between the single and multiple-scattering matrices. The black-box approach plus the dependence of the matrix upon geometry and extension of the medium gives the model a phenomenological character as opposed to a fundamental (first principles) approach. The narrow applicability of a phenomenological approach is therefore in effect. In U.S. Pat. No. 5,818,583, issued Oct. 6, 1998, E. Sevick-Muraca et al describe a method that applies the so-called diffusion model of light propagation to predict some aspects of multiple scattering of electromagnetic waves with suspensions (host medium is a fluid). When the particle concentration is larger than about 1%, a further approximation to the Radiative Transfer Theory can be made leading to the Diffusion Approximation Theory (22 A. Ishimaru, 1978). In this diffusion model, it is assumed that the diffusion intensity is scattered in all directions with an almost uniform angular distribution. In what the inventors describe as photon migration techniques, the light intensity of the source is varied in time and the phase shift of the received signal for a number of different wavelengths is directly measured. The isotropic scattering coefficient is then calculated with the diffusion model and, using this coefficient, the particle size and concentration is determined inverting an integral equation that relates the particle size distribution with the scattering coefficient. Due to the diffusion regime assumption, the model cannot be used at low concentration, because the presence of substantial multiple scattering is a necessity for proper performance. Another limitation is that the referred integral equation relating the scattering coefficient with the size distribution assumes no interaction between particles occurs (the scattering coefficient is proportional to the particle concentration). This equation can only be valid when particles are far apart. As a result, the model will only be accurate when multiple scattering is strong enough for the diffusion approximation to hold but yet the particles being far apart enough so no interaction occurs between them so equation 3 in column 9, line 55 of this prior-art patent can be valid. These are opposite requirements in terms of particle concentration. Consequent with this fact, the experimental data shown in the patent correspond to concentrations around 1%when actual industrial concentrations are much higher than that. This is specifically discussed by the authors of the patent in 46 Photon-Migration Measurement of Latex Size Distribution in Concentrated Suspensions, AIChE Journal, March 1997, Vol. 43, No. 3. In addition, due to the assumption of almost uniform angular distribution of scattering, validity of the model will break down as the ratio between the size of the particles and the wavelength increases. This is because the larger the particle size as compared to the wavelength, the stronger the directivity of the scattering in the forward direction is.

[0010] In U.S. Pat. No. 5,121,629, issued Jun. 16, 1992, I describe a method and apparatus for measuring particle size distribution that employs a fundamental mathematical model for scattering of ultrasound waves in suspensions or emulsions (host medium is a fluid). The modeling aspects of this prior patent actually constitute the most direct prior-art patent to the current invention as far as its application to predicting the attenuation of acoustic waves in suspensions is concerned. Modeling details in this prior patent start in column 8, line 29 up to column 12, line 46. As in the METAMODEL Prediction Engine 21 of the current invention, modeling in this prior patent employs the fundamental wave scattering theory. It starts with the first-principle description of the interaction of a wave with a single particle (the particle/medium signature), and proceeds to statistically formulate the field relations within an ensemble of particles of different sizes. It ends with an analytical equation that relates the propagation constant of the composite with the propagation constant of the host medium, the physical properties of both medium and particles, as well as their size distribution and concentration. There are, however, three essential differences with the current invention when applied to acoustics. Two of them are associated with the particle/medium signature for which this prior-art patent uses the Epstein-Carhart-Allegra-Hawley (ECAH) theory while the present invention adds to this ECAH theory the following: a) modeling for possible viscoelastic behavior of both host medium and particle and b) modeling for phenomena due to surface tension which is important in gaseous particles. The third paramount difference has to do with the overlapping of contiguous fields: the prior patent arrives to equation 8 (column 10, line 60), which is the statistical relation obtained by Waterman and Truell in 9 Multiple Scattering of Waves, J. Math. Phys. Vol. 2, No.4, July-August, 1961, as a closed-form dispersion equation for the particulate composite. Such a simple closed-form equation was obtained under the assumption of statistically independent point scatterers. In plain words, the expression point scatterers refers to the abstract concept of a particle with no physical dimensions but with the scattering properties of its actual size. It was applied in order to simplify the mathematics and its numerical coding into computer software. As a result, model accuracy deteriorates as the size of the particles and/or frequency increase. Regarding the assumption of statistical independence between scatterers positions, it eventually breaks untrue as concentration increases because of the inevitable spatial correlation between particles. Furthermore, P. Lloyd and M. V. Berry (17 Wave propagation through an assembly of spheres IV. Relations between different multiple scattering theories, Proc. Phys. Soc., 1967, Vol. 91), proved there was an error in the referred Waterman and Truell's formula arising from an artificial stratification of the medium into thin slabs. All the referred drawbacks, though not damaging in many cases as demonstrated with the good results disclosed in this prior patent, would instead definitely contribute to destroy the accuracy of a generic model (valid for waves of any nature) if the broadest range of materials, particle size, concentration, and wave frequencies is aimed at. In fact, it will be demonstrated when the Prediction Engine 21 is described in detail that, if scatterers are treated as what they are: finite-size inclusions in the host medium, and the spatial correlation between them is not neglected, it is not possible any longer to arrive at a simple closed-form expression for the dispersion equation. Mathematics and computer programming complexities increase considerablybut prediction accuracy and range of applicability improve drastically.

[0011] In U.S. Pat. No. 6,109,098, issued Aug. 29, 2000, A. Dukhin and P. Goetz describe a device that combines acoustic and electroacoustic spectrometry to measure particle size distribution and zeta potential for concentrated dispersed systems. Focusing only on the mathematical modeling part of this prior patent, they use a macroscopic phenomenological approach to modeling the interaction between the sound wave and the concentrated suspension (host medium is a fluid). Specifically, they use the so-called coupled phase and cell models applying the long-wave regime condition, i.e. the particle size has to be much smaller than the sound wavelength. As a result, their invention delivers useful results only for particles with sizes smaller than 10 m. It is worth noting that the authors categorically reject beforehand the very fundamental approach adopted in the present invention when discussing U.S. Pat. No. 5,121,629 by saying in column 4, line 55: . . . He (Alba) suggests using EKAH (they mean ECAH) theory in combination with a multiple scattering approach This is not adequate because even the multiple scattering approach requires a single particle acoustic field which is known only for a single particle in infinite media . . . FIG. 3 through FIG. 11 shown in the current specification will irrefutably demonstrate that the Prediction Engine 21 of the present invention has accomplished what U.S. Pat. No. 6,109,098 says it was unfeasible.

[0012] In U.S. Pat. No. 6,119,510, issued Sep. 19, 2000, M. Carasso. et al describe an improved process for determining the characteristics of dispersed particles. The authors contemplate the use of their process with acoustic as well as light waves. Focusing again only in the mathematical modeling, they describe a computationally improved Allegra-Hawley Model for acoustics and the original Mie Model for light waves. As already stated in depth, those models are single-scattering models that will inevitably fail as the concentration increases (typically around 10% for acoustics and 1% for optics). No attempt whatsoever is made to model multiple-scattering phenomenanot even by employing equation 8 (column 10, line 60) of U.S. Pat. No 5,121,629 that is the closest established prior-art per their own admission. The authors simply say: Although this predicted spectrum is an approximation (owing to non-linearities induced by multiple-scattering phenomena) it provides desirable results in this invention. The meaning of the word desirable is nonetheless never tested because the acoustic experimental data shown in this prior patent correspond to concentrations of 1% and 5% for which multiple scattering is negligible. Similarly, the optical turbidity data correspond to a concentration of only 1%.

[0013] In summary, there are no research and development tools available for the generic prediction, analysis, and parametric studies of the interaction between generic waves (acoustic, electromagnetic, elastic) and particulate composites. The current invention meets such a long-term need by providing a powerful research and development tool that makes all theoretical physics and numerical intricacies transparent to the human operator. A fundamental mathematical model is a key component of such a tool. The lack of a complete fundamental mathematical model for the prediction of multiple scattering of waves in particulates has, in addition, limited the concentration range of currently existing wave-based instruments. All such instruments which operation is based on the interaction of any type of wave traveling through suspensions, emulsions or solids containing inclusions, could benefit considerably from the prediction engine of this invention. Because the METAMODEL Prediction Engine 21 is a separate module connected to the other modules of this invention through industry-standard interfaces, it can be directly integrated to thosealready existinginstruments extending their concentration ranges several orders of magnitude. Further objects and advantages of my invention will become apparent from a consideration of the ensuing description.

[0014] This invention comprises of three crucial components (FIG. 1): a) a fundamental mathematical model, with its corresponding numerical algorithms to predict the composite physical attributes from the physical attributes of the constituent phases. I call this component the METAMODEL Prediction Engine 21, b) a set of optimization algorithms to numerically invert the aforesaid mathematical model so as to estimate a subset of the attributes of the constituent phases and composite from the knowledge of the rest of the attributes of the constituent phases and composite. I call this component the METAMODEL Inversion Engine 22; and c) an ergonomic graphical interface between the prediction and inversion engines and the final user so he/she can easily set up complex simulation experiments to assist in parametric studies, analysis, and synthesis of real systems in which waves of any nature interact with composites. I call this component the SCATTERER Graphical Interface 23. The Prediction Engine 21 uses a fundamental approach that considers the composite as a random medium and consequently describes it statistically as indicated in block 24 of FIG. 2. The particle/host medium signature is deterministically modeled from the knowledge of the physical properties of particle and host medium and the physics associated with the nature of the problem (acoustical, electromagnetic, elastic) as expressed in block 25. Particle spatial and attribute correlations are modeled as indicated in block 26. Stochastic field equations, assuming a particle is fixed in space with its physical attributes known, are established as represented in block 27. Further averaging over the physical attributes of the fixed particle is performed as expressed in block 28, arriving to the dispersion equations that are solved to determine the propagation constant of the random medium as indicated in block 29. All physical attributes of the composite are thus determined from the composite propagation constant and the physical attributes of the constituent phases as depicted in block 30. The graphical interface 23 is implemented following de facto industry standards for graphical user interfaces. Hierarchies of windows, icons, options and command buttons, checkbox controls, context menus, image controls, tool-tips, on-line help and the like are used to produce an ergonomic efficient interaction between this invention and the human operator. A general-purpose computer 19 and data-carrier media 20 are needed for the full operability of the invention.

[0015] FIG. 1 schematically depicts main components of the invention with their structural relationship.

[0016] FIG. 2 schematically displays main modeling steps in the METAMODEL prediction Engine.

[0017] FIG. 3 validates the prediction engine for an extremely fine dispersion of silica particles (19 nm) used in the semiconductor industry for polishing, at a particle concentration of C_v11.4%

[0018] FIG. 4 validates the prediction engine for an extremely fine dispersion of silica particles (19 nm) used in the semiconductor industry for polishing at C_v17.1%

[0019] FIG. 5 validates the prediction engine for a dispersion of standard silica particles (100 nm30 nm) at C_v10.0%

[0020] FIG. 6 validates the prediction engine for a dispersion of standard silica particles (100 nm30 nm) at C_v30.0%

[0021] FIG. 7 validates the prediction engine for a dispersion of standard silica particles (300 nm30 nm) at C_v16.0%

[0022] FIG. 8 validates the prediction engine for a dispersion of standard silica particles (300 nm30 nm) at C_v24.0%

[0023] FIG. 9 validates the prediction engine for a dispersion of standard silica particles (450 nm30 nm) at C_v10.0%

[0024] FIG. 10 validates the prediction engine for a dispersion of standard silica particles (450 nm30 nm) at C_v30.0%

[0025] FIG. 11 validates the prediction engine for a broad distribution of soy oil droplets in water (0.06 m to 2 n) at C_v10.6%

[0026] Detailed Description/Operation of the METAMODEL Engines 21 and 22

[0027] An essential part of this invention is a generic fundamental mathematical model to predict the interaction between a generic wave (acoustic, electromagnetic, elastic) and a dense generic particulate system (fluid or solid particles included in fluid or solid media). Two numerical engines operate around this mathematical model. The Prediction Engine 21 starts with the operating attribute vector o, the attribute vector m of the host medium, and the attribute vector a of the particles to predict the attribute vector c of the composite (FIG. 1). Vector o contains operating parameters like frequency, temperature, scattering angle, etc. Vector m contains all the physical properties for the suspending medium. Vector a comprises of all the physical properties for the material of the particles plus the particle size distribution and concentration. Some of these properties can be frequency-dependent (e.g. for non-Newtonian fluid or viscoelastic solids, or the refractive index in the optical case). Properties may also be temperature-dependent. The components of vector c are the physical attributes of the composite (attenuation, velocity, intensity, phase shift, elastic constants, etc.). These attributes are in general functions of frequency, temperature, scattering angle, etc. As a clarifying example, for the acoustic case treated in U.S. Pat. No. 5,121,629 the knowledge of the seven disclosed physical properties for both the suspending medium and suspended particles plus the particle size distribution and particle concentration, made it possible to uniquely predict the attenuation coefficient of the ultrasound wave as a function of frequency. The Inversion Engine 22 handles the inverse problem, i.e. knowing some of the attributes of the composite and constituents, it estimates the values of the missing attributes of the constituents and composite (FIG. 1). The symbol P. stands for a subset of all properties (the partial knowledge set). The symbol {overscore (P)}. stands for the missing subset (set complement). As a concrete example, most of the referred patents describe particle size analyzers, viz. instruments that estimate the particle size distribution (a constituent attribute) from the measurement of a composite attribute (attenuation spectrum for U.S. Pat. No. 5,121,629 and U.S. Pat. No. 6,109,098, light intensity as a function of angle for U.S. Pat. No. 5,619,324, and phase shift as a function of wavelength for U.S. Pat. No. 5,818,583) and the knowledge of the rest of the attributes. All the inherent numerical intricacies associated with the inversion of integral equations are well known by those skilled in the art. Non-linear optimization algorithms using Powell or related techniques in multi-dimensional space (W. H. Press et al, Numerical Recipes, The Art of Scientific Computing, Cambridge Press, 1988) have been successfully employed for the specific acoustic problem described in U.S. Pat. No. 5,121,629. The present invention, using similar but more powerful and elaborate numerical techniques, provides for a generalized direct/invert algorithmic scheme allowing the human operator to manipulate the mathematical model in any desired direction specifying those parameters that are known as well as those that are not known and are to be estimated. As an example, the engineer/scientist could set up the current invention to use measured attenuation spectral data, particle size, and concentration as inputs to estimate the elastic constants of the particles.

[0028] In the attached embodiment for this invention, both numerical engines 21 and 22 are programmed in Compaq Visual FORTRAN as a Dynamically Linked Library (DLL) called by the Graphical Interface 23 that is programmed in Microsoft Visual Basic for Windows. Any other language or combination of languages interacting with a different liking technology, e.g. Component Object Model (COM) server technology, could obviously be used to embody the current invention.

[0029] In this generic modeling problem I start considering a random spatial distribution of N particles (scatterers) with arbitrary physical attributes (thermodynamical, size, shape, orientation, acoustic, elastic, electromagnetic, transport, etc.) suspended or included in an otherwise continuous medium with an arbitrary concentration. The particular position and attributes of each scatterer are unknown but the ensemble of them can be statistically described through a probability density function p(r₁, . . . ,r_i, . . . ,r_N;a₁, . . . ,a_i, . . . ,a_N). The symbols r_i and a_i represent the position and attribute vectors for particle i respectively.

[0030] In plain words: p(r₁, . . . ,r_i, . . . ,r_N;a₁, . . . ,a_i, . . . ,a_N)dv₁ . . . dv_i . . . dv_Nda₁ . . . da_i . . . da_N is the joint probability of finding one scatterer in the volume dv₁ (centered at r₁) with attribute vector in the neighborhood of a₁ (aa₁), another scatterer in dv₂ with attribute vector in the neighborhood of a₂, etc. The symbol dv_i stands for the volume differential element in R³ with R being the set of real numbers. Calling nA the number of attributes for each particle, note that the symbol da_i stands not for the differential of the vector a_i, but for the hyper-volume differential element in R^nA. One of the components of a_i is the size of particle i that will be denoted as d_i. The simplest practical case commonly considered in the prior multiple-scattering literature is that in which attributes of all scatterers are identical. A few authors considered as well the case in which there is a distribution of sizes throughout the particle population.

[0031] A generic wave incites the above random medium. As it is well known from the prior art, assuming the particle concentration is low enough, the local exciting wave at each scattering center can be considered to be simply the global incident wave (Born's approximation). This is the well-known single-scattering regime where the behavior of the ensemble is simply proportional to the behavior of the individual particles and the number of them per unit volume. For the acoustical case, this incident wave is normally a longitudinal wave (compressional), i.e., its displacement vector is parallel to the direction of propagation of the wave. At each scattering center this longitudinal wave converts into multiple modes inside as well as outside the particle: two compressional waves (refracted and reflected), two thermal longitudinal waves (inside and outside), and two viscous-inertial transversal waves (inside and outside). At low concentrations these new modes extinguish in the close neighborhood of each particle before interacting with other particles. As the concentration increases, these local waves around the particles overlap and multiple-scattering effects become significant. The detected wave at the receiver is the remaining propagating compressional mode. For the electromagnetic case, the incident wave is transversal, i.e., the plane of vibration of the electrical and magnetic fields is perpendicular to the direction of propagation of the wave. There are normally no conversion modes at each scatterer. The elastodynamic case is the most complex of all, as its vectorial nature is present in both the incident wave as well as in the scattering process that occurs at each particle. The incident wave can have both longitudinal and transversal components and each one of them convertsat each scattering centerinto longitudinal and transversal modes.

[0032] The incident field is assumed to be that of a monochromatic wave without loss of generality. The factor e^-j is consequently implicit in all equations. Except for the incident field, all fields in the particulate system depend upon the spatial/attribute configuration of scatterers r₁, . . . ,r_i, . . . ,r_N;a₁, . . . ,a_i, . . . ,a_N. This lengthy notation will be maintained throughout the present description for the sake of an enabling disclosure. I call the incident field at observation point r by I(r), and the field that excites particle i (located at point r_i with attribute vector a_i) at observation point r by E(r/r_i;a_i:r₁, . . . ,r_i, . . . ,r_N;a₁, . . . ,a_i, . . . ,a_N). Similarly, calling the field scattered by particle j (located at r_j with attribute a_j) when particle i with attribute a_i is located at r_i by NS(r/r_j;r_i;a_j,a_i:r₁, . . . ,r_i, . . . ,r_j, . . . ,r_N; a₁, . . . ,a_i, . . . ,a_j, . . . ,a_N), the following important generic equation is valid: $\begin{array}{cc}\underset{\_}{E}\left(\underset{\_}{r}/{\underset{\_}{r}}_i;{\underset{\_}{a}}_i:\left\{{\underset{\_}{r}}_1,,{\underset{\_}{r}}_i,,{\underset{\_}{r}}_N;{\underset{\_}{a}}_1,,{\underset{\_}{a}}_i,,{\underset{\_}{a}}_N\right\}\right)=\underset{\_}{I}\left(\underset{\_}{r}\right)+\underset{j=1,ji}{\overset{N}{}}\underset{\_}{S}\left(\underset{\_}{r}/{\underset{\_}{r}}_j{\underset{\_}{r}}_i;{\underset{\_}{a}}_j{\underset{\_}{a}}_i:\left\{{\underset{\_}{r}}_1,,{\underset{\_}{r}}_i,,{\underset{\_}{r}}_j,,{\underset{\_}{r}}_N;{\underset{\_}{a}}_1,,{\underset{\_}{a}}_i,,{\underset{\_}{a}}_j,,{\underset{\_}{a}}_N\right\}\right)& \left(1\right)\end{array}$

[0033] Following Foldy's 2 approach to handle the statistical nature of the physical system under consideration, I define the expected (average) value of a field at the observation point r as the configurational average over all possible spatial/attribute distributions of the scatterers: $\begin{array}{cc}\underset{\_}{F}\left(\underset{\_}{r}\right)={}_V^{}{}_V^{}{}_D^{}{}_D^{}\underset{\_}{F}(\underset{\_}{r}:\left\{{\underset{\_}{r}}_1,,{\underset{\_}{r}}_i,,{\underset{\_}{r}}_N;{\underset{\_}{a}}_1,,{\underset{\_}{a}}_i,{\underset{\_}{a}}_N\right\})p\left({\underset{\_}{r}}_1,,{\underset{\_}{r}}_i,{\underset{\_}{r}}_N;{\underset{\_}{a}}_1,,{\underset{\_}{a}}_i,,,{\underset{\_}{a}}_N\right)v_1{v}_i{v}_N{\underset{\_}{a}}_1{\underset{\_}{a}}_{{}_i}{\underset{\_}{a}}_N\underset{\_}{F}\left(\underset{\_}{I},\underset{\_}{E},\underset{\_}{S}\right\}& \left(2\right)\end{array}$

[0034] The domain of integration V is the volume accessible to the scatterers, i.e., a subset of R³. The domain D of the attributes is a subset of R^nA. Because of the random nature of the system all fields can be thus expressed as:

[0035] Symbols <F(r)> and F₁ (r) represent the coherent and incoherent fields respectively. The incoherent field is the truly random part of the total field. Through the introduction of conditional probabilities, the expected value of a field observed at r when a particle is assumed to be at r_i with attribute vector a_i (first order conditional average) is: $\begin{array}{cc}\underset{\_}{F}\left(\underset{\_}{r}/{\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right)={}_{\mathrm{Vi}}^{}{}^{}{}_{\mathrm{Vi}}^{}{}_D^{}{}^{}{}_D^{}\underset{\_}{F}\left(\underset{\_}{r}/{\underset{\_}{r}}_i;{\underset{\_}{a}}_i:\left\{{\underset{\_}{r}}_1,,{\underset{\_}{r}}_i,,{\underset{\_}{r}}_N;{\underset{\_}{a}}_1,,{\underset{\_}{a}}_i,,{\underset{\_}{a}}_N\right\}\right)p\left({\underset{\_}{r}}_1,{,}^{}{\underset{\_}{r}}_N;{\underset{\_}{a}}_1,{,}^{},,{\underset{\_}{a}}_N/{\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right)v_1,{}^{}{v}_N{\underset{\_}{a}}_1{}^{}{\underset{\_}{a}}_N\underset{\_}{F}\left\{\underset{\_}{I},\underset{\_}{E},\underset{\_}{S}\right\}& \left(4\right)\end{array}$

[0036] The primes in the previous formula mean that r_i and a_i are not varied throughout the averaging process as clearly expressed by the conditional probability. The domain of integration Vi is V minus a sphere of diameter d_i. For the case of assuming particles i and j respectively at locations r_i and r_j with properties a_i and a_j (second order conditional average) it is obtained: $\begin{array}{cc}\underset{\_}{F}\left(\underset{\_}{r}/{\underset{\_}{r}}_i,{\underset{\_}{r}}_j;{\underset{\_}{a}}_i,{\underset{\_}{a}}_j\right)={}_{\mathrm{Vij}}^{}{}^{}{}^{}{}_{\mathrm{Vij}}^{}{}_D^{}{}^{}{}^{}{}_D^{}\underset{\_}{F}\left(\underset{\_}{r}/{\underset{\_}{r}}_i,{\underset{\_}{r}}_j;{\underset{\_}{a}}_i,{\underset{\_}{a}}_j\right):\left\{{\underset{\_}{r}}_1,,{\underset{\_}{r}}_i,,{\underset{\_}{r}}_j,,{\underset{\_}{r}}_N;{\underset{\_}{a}}_1,,{\underset{\_}{a}}_i,,{\underset{\_}{a}}_j,,{\underset{\_}{a}}_N\right\})p\left({\underset{\_}{r}}_1{,}^{},{,}^{},{\underset{\_}{r}}_N;{\underset{\_}{a}}_1,{,}^{},{,}^{},,{\underset{\_}{a}}_N/{\underset{\_}{r}}_i,{\underset{\_}{r}}_j;{\underset{\_}{a}}_i,{\underset{\_}{a}}_j\right)v_1{}^{}{}^{}{v}_N{\underset{\_}{a}}_1{}^{}{}^{}{\underset{\_}{a}}_N\underset{\_}{F}\left\{\underset{\_}{I},\underset{\_}{E},\underset{\_}{S}\right\}& \left(5\right)\end{array}$

[0037] Without further explanation, statistics of higher order are straightforwardly defined. The global configurational average <F(r)> can then be expressed in terms of the first partial average as follows: $\begin{array}{cc}\underset{\_}{F}\left(\underset{\_}{r}\right)={}_V^{}{}_D^{}\underset{\_}{F}\left(\underset{\_}{r}/{\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right)p\left({\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right){\underset{\_}{a}}_i{v}_i\underset{\_}{F}\left\{\underset{\_}{I},\underset{\_}{E},\underset{\_}{S}\right\}& \left(6\right)\end{array}$

[0038] Similarly, the first partial average can be expressed in terms of the second as: $\begin{array}{cc}\underset{\_}{F}\left(\underset{\_}{r}/{\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right)={}_{\mathrm{Vi}}^{}{}_D^{}\underset{\_}{F}\left(\underset{\_}{r}/{\underset{\_}{r}}_i,{\underset{\_}{r}}_j;{\underset{\_}{a}}_i,{\underset{\_}{a}}_j\right)p\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j/{\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right){\underset{\_}{a}}_j{v}_j={}_{\mathrm{Vi}}^{}{}_D^{}\underset{\_}{F}\left(\underset{\_}{r}/{\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)p\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j/{\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right){\underset{\_}{a}}_j{v}_j\underset{\_}{F}\left\{\underset{\_}{I},\underset{\_}{E},\underset{\_}{S}\right\}& \left(7\right)\end{array}$

[0039] In general then, the partial average of order k can be expressed in terms of the partial average of order k1 and the corresponding conditional probability. Pursuing this recurrence throughout the whole population is obviously impractical. In order to break the recurrence, I used in the above second equality a generalization of the so-called quasi-crystalline approximation introduced first by Lax 3,4. It consists in assuming that <F(r/r_i,r_j;a_i,a_j(><F(r/r_j;a_j)>. Breaking the recurrence at a higher statistical level is also contemplated in this invention but it would complicate considerably more the mathematical modeling process as well as its disclosure in this specification with no significant practical advantage.

[0040] Taking now the configurational average of both terms in equation 1 and using equation 7, it is thus obtained: $\begin{array}{cc}\underset{\_}{E}\left(\underset{\_}{r}/{\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right)=\underset{\_}{I}\left(\underset{\_}{r}\right)+\underset{j=1,j1}{\overset{N}{}}{}_{\mathrm{Vi}}^{}{}_D^{}\underset{\_}{S}\left(\underset{\_}{r}/{\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)p\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j/{\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right){\underset{\_}{a}}_j{v}_j& \left(8\right)\end{array}$

[0041] Let us now define n(r_j,a_j)da_jdv_j as the number of scatterers in the neighborhood of point r_j and with attribute vector a in the neighborhood of a_j. Let us further define p(r_j;a_j)da_jdv_jn(r_j,a_j)da_jdv_j/N, i.e., the ratio between the number of particles with attribute in the neighborhood of a_k and location in the neighborhood of r_j, and the total number of particles in the ensemble. Assuming there is no global correlation between position and attribute, it is possible to simplify the joint probability as p(r_j;a_j)n(r_j)n(a_j)/N²n(r_j)f_n(a_j)/N. We have introduced the function f_n(a_j) as the density distribution by number for the attributes of particle j in the neighborhood of r_j. When the attribute is particle size, f_n(a_j) is the standard particle size distribution by number but local to the neighborhood of r_j. The non-segregation assumption also implies f_n (a_j)f_n (a_i)f_n (a) which is now the standard ensemble particle size distribution. Along similar intuitive lines, the conditional probability can be defined as:

[0042] The function g(r_j,a_j,r_i,a_i) is a measure of the local correlation between the positions and attributes of the particles. It basically depends on the nature and range of the interparticle forces. A very successful approximation commonly used in statistical mechanics (the so-called Pair Correlation Function) consists in assuming particles are impenetrable and that the correlation is only a function of their distance rr_jr_i and the sum of their radii D(d_jd_i)/2:

[0043] The function H(x) is the Heaviside step function which value is 0 for x<0 and 1 for x0. When CF(r)1 it is the simplest approximation to correlation: as long as particles do not interpenetrate, given that particle i is located at r_i, the chance of particle j to be located at r_j with attribute a_j is governed exclusively by p(r_j,a_j). This is called the hole correction in the prior art 13. In general, the function CF(r) is oscillatory approaching unity as the distance between the two particles increases simply because as two particles separate away, their spatial correlation disappear. There is an extended literature regarding the calculation of CF(r). The most successful and well-known form is the so-called Percus-Yevick approximation from the theory of classical fluids. Other approaches are the so-called hypernetted chain equation, the Virial expansion, the self-consistent approximation, etc. A powerful technique recently developed is based on Monte-Carlo simulation of a large number of scatterer spatial configurations. The current invention, in its graphical user interface, provides for most of these well-known forms for CF(r) as set-up options. Another powerful option provided by this invention is to allow for this oscillatory function to be characterized by three parameters (amplitude, frequency, and decay) and let the user experiment with them in trying to accurately modeling the unique spatial correlation characteristics of the problem at hand. The following relation is then valid: $\begin{array}{cc}\underset{\_}{E}\left(\underset{\_}{r}/{\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right)=\underset{\_}{I}\left(\underset{\_}{r}\right)+{}_{\mathrm{Vi}}^{}{}_D^{}\underset{\_}{S}\left(\underset{\_}{r}/{\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)n\left({\underset{\_}{r}}_j,{\underset{\_}{a}}_j\right)g\left({\underset{\_}{r}}_j,{\underset{\_}{a}}_j,{\underset{\_}{r}}_i,{\underset{\_}{a}}_i\right){\underset{\_}{a}}_j{v}_j& \left(11\right)\end{array}$

[0044] The field scattered by a single particle has been separately studied in the prior literature for all acoustic, electromagnetic, and elastic cases. The most general form of expressing the scattered field generated by scatterer j due to its exciting field is through the so-called transition matrix operator as follows:

[0045] In plain words, the transition matrix constitutes a behavioral signature for the basic individual interaction between the wave and a single particle suspended in the medium. I call it the particle/medium signature. Notice that both the exciting and scattered fields depend upon the spatial/attribute configuration of scatterers, while the transition operator only depends upon the host medium properties, the particle j properties, the operating parameters (frequency, temperature, etc) and the nature of the incident physical wave (acoustic, electromagnetic, elastic). It is a representation of the Green's function for a single scatterer in the host medium (single scattering problem). When the incident wave is electromagnetic, this invention employs the well-known Mie-Scattering model to define the transition operator. In the acoustic case the transition operator is obtained from: 1) 5 P. S Epstein and R. R. Carhart, The Absorption of Sound in Suspensions and Emulsions, I. Water Fog in Air, J. Acoust. Soc. Am., 25 3 553-565 (1953), 2) A correction to the thermal diffusion effect due to surface tension is included from the work of J. C. F. Chow, 14 Attenuation of Acoustic Waves in Dilute Emulsions and Suspensions, J. Acoust. Soc. Am., 36 12 2395-2401 (1964), 3) A more complete treatment for both solid and fluid particles is taken from J. R. Allegra and S. A. Hawley, 19 Attenuation of Sound in Suspensions and Emulsions: Theory and Experiments, J. Acoust. Soc. Am., 51, 1545-1564 (1972), and 4) Viscoelastic behavior of both particles and host medium is modeled along the lines described in C. Verdier and M. Piau, 47 Acoustic Wave Propagation in Two-phase viscoelastic fluids: The Case of Polymer emulsions, J. Acoust. Soc. Am. 1014, 1868-1876, April 1997 and J. Y. Kim et al, 44 Dispersive Wave Propagation in a Viscoelastic Matrix Reinforced by Elastic Fibers, J. Acoust. Soc. Am., 95 (3), March 1994. The transition operator for elastic waves is obtained through a combination of the treatment described for the acoustic case and the approach described in 1) N. G. Einspruch and R. Truell, 7 Scattering of a Plane Longitudinal Wave by a Spherical Fluid Obstacle in an Elastic Medium, J. Acoust. Soc. Am., 32 2, 214-220, February 1960, and 2) N. G. Einspruch et al, 8 Scattering of a Plane Transverse Wave by a Sherical Obstacle in an Elastic Medium, J. Appl. Phys., 31 5, 806-817, May 1960.

[0046] In terms of this basic signature and equation 11, the global behavior of the ensemble (multiple scattering) will be predicted. Taking the configurational average on both sides of equation 12 it is obtained:

[0047] Replacing equation 13 in equation 11 the following important relation is valid: $\begin{array}{cc}\underset{\_}{E}\left(\underset{\_}{r}/{\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right)=\underset{\_}{I}\left(\underset{\_}{r}\right)+{}_{\mathrm{Vi}}^{}{}_D^{}\underset{\_}{\underset{\_}{T}}\left(\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_j\right)\underset{\_}{E}\left(\underset{\_}{r}/{\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)n\left({\underset{\_}{r}}_j,{\underset{\_}{a}}_j\right)g\left({\underset{\_}{r}}_j,{\underset{\_}{a}}_j,{\underset{\_}{r}}_i,{\underset{\_}{a}}_i\right){\underset{\_}{a}}_j{v}_j& \left(14\right)\end{array}$

[0048] Given a fixed spatial/attribute distribution of scatterers, all fields (external and internal to the particles) are supposed to verify the Helmholtz's vector wave equation:

[0049] The symbol ² stands for the Laplacian of a vectorial field and k is the wave number for the particular field. The general solution to this equation may be expressed as a linear combination of three independent discrete sets of vector solutions L_nm (longitudinal), M_nm (transversal) and N_nm (transversal) with arbitrary constants f_nm1,f_nm2,f_nm3 as follows 1: $\begin{array}{cc}\underset{\_}{F}\left(\underset{\_}{r}\right)=\underset{n=0}{\overset{}{}}\underset{m=-n}{\overset{m=n}{}}\left[f_{\mathrm{nm1}}{L}_{\mathrm{nm}}\left(\underset{\_}{r}\right)+f_{\mathrm{nm2}}{M}_{\mathrm{nm}}\left(\underset{\_}{r}\right)+f_{\mathrm{nm3}}{N}_{\mathrm{nm}}\left(\underset{\_}{r}\right)\right]\underset{\_}{F}\left\{I,E,S\right\}& \left(16\right)\end{array}$

[0050] Purely irrotational fields could be expressed exclusively in terms of the L_nm eigenvectors; purely solenoidal fields could be expressed exclusively in terms of both M_nm and N_nm eigenvectors. All three types of basis vector functions are directed related to the discrete set of solutions (r)_nm for the corresponding scalar Helmholtz wave equation 1. For particles with shapes randomly deviating from the spherical shape, spherical coordinates are appropriate. For fiber-shaped particles, cylindrical coordinates are the suitable choice. For the sake of continuing the disclosure, using spherical coordinates (r,,), the associated discrete sets L_nm M_nm N_nm are:

[0051] The P_n^m(cos) set comprises of the well known associated Legendre Polynomials; and the z_n set is either the set of spherical Bessel functions j_n or the set of Hankel functions h_n. The Bessel set is chosen for those fields that are to be regular at the origin (incident field and the internal field of a particle); the Hankel functional set is used for those fields that are to verify the radiation condition at infinity (the scattered fields). Superscripts B or H will be added to L_nm,M_nm,N_nm according to the case.

[0052] Let the symbol Ndw stand for Number of different waves occurring during the interaction between the exciting wave and the particle. For instance, in the acoustic problem, as described in U.S. Pat. No. 5,121,629, Ndw would be equal to 3 with the value 1 corresponding to the compressional waves, value 2 corresponding to the longitudinal thermal waves, and the value 3 corresponding to the viscous transversal waves. Adding then a new subscript k1, Ndw, and defining k_k as the corresponding wave number, the following equalities are valid in spherical coordinates: $\begin{array}{cc}\begin{array}{c}{L}_{\mathrm{nmk}}\left(\underset{\_}{r}\right)=\left\{\frac{}{r}{z}_n\left(k_{k}r\right)P_n^m\left(\cos \right){}^m\right\}{\underset{\_}{e}}_r+\\ \left\{z_n\left(k_{k}r\right)\frac{}{rr}{P}_n^m\left(\cos \right){}^m\right\}{\underset{\_}{e}}_{}+\\ \left\{\frac{m}{r\sin }{z}_n\left(k_{k}r\right)P_n^m\left(\cos \right){}^m\right\}{\underset{\_}{e}}_{}\\ M_{\mathrm{nmk}}\left(\underset{\_}{r}\right)=\left\{\frac{m}{\sin }{z}_n\left(k_{k}r\right)P_n^m\left(\cos \right){}^m\right\}{\underset{\_}{e}}_{}+\\ \left\{z_n\left(k_{k}r\right)\frac{}{}{P}_n^m\left(\cos \right){}^m\right\}{\underset{\_}{e}}_{}\\ N_{\mathrm{nmk}}\left(\underset{\_}{r}\right)=\left\{\frac{n\left(n+1\right)}{k_{k}r}{z}_n\left(k_{k}r\right)P_n^m\left(\cos \right){}^m\right\}{\underset{\_}{e}}_r+\\ \left\{\frac{}{k_{k}rr}{\mathrm{rz}}_n\left(k_{k}r\right)\frac{}{}{P}_n^m\left(\cos \right){}^m\right\}{\underset{\_}{e}}_{}+\\ \left\{\frac{}{k_{k}rr}{\mathrm{rz}}_n\left(k_{k}r\right)\frac{m}{\sin }{P}_n^m\left(\cos \right){}^m\right\}{\underset{\_}{e}}_{}\end{array}& \left(18\right)\end{array}$

[0053] It is worth noting that equations 5 and 6 of prior art U.S. Pat. No. 5,121,629 are special versions of our generic equation 16 in spherical coordinates (upon use of equation 18) for the simpler acoustic case where the incident wave is purely longitudinal (a field of gradients). Let us now apply our equation 16 to the exciting field surrounding particle i located at position r_i as well as to the field scattered by particle j located at position r_j, and take configurational averages: $\begin{array}{cc}\begin{array}{c}\underset{\_}{E}\left(\underset{\_}{r}/{\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right)=\underset{n=0}{\overset{}{}}\underset{m=-n}{\overset{m=n}{}}\underset{k=1}{\overset{\mathrm{Ndw}}{}}[e_{\mathrm{nmk}1}\left({\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right)L_{\mathrm{nmk}}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)+\\ e_{\mathrm{nmk}2}\left({\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right)M_{\mathrm{nmk}}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)+\\ e_{\mathrm{nmk}3}\left({\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right)N_{\mathrm{nmk}}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)]\\ \underset{\_}{r}:\underset{\_}{r}-{\underset{\_}{r}}_{i}b\end{array}& \left(19\right)\\ \begin{array}{c}\underset{\_}{S}\left(\underset{\_}{r}/{\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)=\underset{n=0}{\overset{}{}}\underset{m=-n}{\overset{m=n}{}}\underset{k=1}{\overset{\mathrm{Ndw}}{}}[s_{\mathrm{nmk}1}\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)L_{\mathrm{nmk}}^H\left(\underset{\_}{r}-{\underset{\_}{r}}_j\right)+\\ s_{\mathrm{nmk}2}\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)M_{\mathrm{nmk}}^H\left(\underset{\_}{r}-{\underset{\_}{r}}_j\right)+\\ s_{\mathrm{nmk}3}\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)N_{\mathrm{nmk}}^H\left(\underset{\_}{r}-{\underset{\_}{r}}_j\right)]\\ \underset{\_}{r}:\underset{\_}{r}-{\underset{\_}{r}}_{j}b\end{array}& \left(20\right)\end{array}$

[0054] The symbol b stands for the distance of closest approach between particles i and j (b(d_jd_i)/2). Following the expansions in equations 19 and 20, applying a suitable set of boundary conditions around particle j, and adding the polarization index q1,2,3, the abstract nature of the operator T(o,m,a_j) materializes into a multi-dimensional infinite matrix that relates each scattering coefficient in equation 20 with all the coefficients in equation 19 as: $\begin{array}{cc}\begin{array}{c}{s}_{\mathrm{nmkq}}\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)=\underset{r=0}{\overset{}{}}\underset{s=-r}{\overset{s=r}{}}\underset{t=1}{\overset{\mathrm{Ndw}}{}}\underset{u=1}{\overset{3}{}}{T}_{\mathrm{rstu}}^{\mathrm{nmkq}}\left(\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_j\right)e_{\mathrm{rstu}}\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)\\ =\underset{\mathrm{rstu}}{}{T}_{\mathrm{rstu}}^{\mathrm{nmkq}}\left(\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_j\right)e_{\mathrm{rstu}}\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)\end{array}& \left(21\right)\end{array}$

[0055] In the electromagnetic case, due to the purely solenoidal character of the wave and the spherical symmetry involved, only T_n112ⁿ¹¹² (o,m,p) and T_n113ⁿ¹¹³ (o,m,p_j) are needed and they coincide with the well known Mie coefficients an and b, in the prior art 27. In the acoustic and elastic cases both longitudinal and transversal waves exist. Applying the single summation convention used in equation 21, equation 20 now becomes: $\begin{array}{cc}\begin{array}{c}\underset{\_}{S}\left(\underset{\_}{r}/{\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)=\underset{\mathrm{nmk}}{}[\left\{\underset{\mathrm{rst}}{}{T}_{\mathrm{rst}1}^{\mathrm{nmk}1}\left(\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_j\right)e_{\mathrm{rst}1}\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)\right\}\\ L_{\mathrm{nmk}}^H\left(\underset{\_}{r}-{\underset{\_}{r}}_j\right)+\underset{\mathrm{rst}}{}{T}_{\mathrm{rst}2}^{\mathrm{nmk}2}\left(\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_j\right)\\ e_{\mathrm{rst}2}\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)\}{M}_{\mathrm{nmk}}^H\left(\underset{\_}{r}-{\underset{\_}{r}}_j\right)+\\ \left\{\underset{\mathrm{rst}}{}{T}_{\mathrm{rst}3}^{\mathrm{nmk}3}\left(\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_j\right)e_{\mathrm{rst}3}\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)\right\}\\ N_{\mathrm{nmk}}^H\left(\underset{\_}{r}-{\underset{\_}{r}}_j\right)]\end{array}& \left(22\right)\end{array}$

[0056] And, consequently, equation 14 converts to: $\begin{array}{cc}\underset{\mathrm{nmk}}{}\left[e_{\mathrm{nmk}1}\left({\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right)L_{\mathrm{nmk}}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)+e_{\mathrm{nmk}2}\left({\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right)M_{\mathrm{nmk}}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)+e_{\mathrm{nmk}3}\left({\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right)N_{\mathrm{nmk}}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)\right]=\underset{\mathrm{nmk}}{}\left[i_{\mathrm{nmk}1}\left({\underset{\_}{r}}_i\right)L_{\mathrm{nmk}}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)+i_{\mathrm{nmk}2}\left({\underset{\_}{r}}_i\right)M_{\mathrm{nmk}}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)+i_{\mathrm{nmk}3}\left({\underset{\_}{r}}_i\right)N_{\mathrm{nmk}}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)\right]+{}_{\mathrm{Vi}}{}_{\mathrm{Di}}n\left({\underset{\_}{r}}_j,{\underset{\_}{a}}_j\right)g\left({\underset{\_}{r}}_j,{\underset{\_}{a}}_j,{\underset{\_}{r}}_i,{\underset{\_}{a}}_i\right)\underset{\mathrm{nmk}}{}\left[\left\{\underset{\mathrm{rst}1}{}{T}_{\mathrm{rst}1}^{\mathrm{nmk}1}\left(\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_j\right)e_{\mathrm{rst}1}\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)\right\}{L}_{\mathrm{nmk}}^H\left(\underset{\_}{r}-{\underset{\_}{r}}_j\right)+\left\{\underset{\mathrm{rst}}{}{T}_{\mathrm{rst}2}^{\mathrm{nmk}2}\left(\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_j\right)e_{\mathrm{rst}2}\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)\right\}{M}_{\mathrm{nmk}}^H\left(\underset{\_}{r}-{\underset{\_}{r}}_j\right)+\left\{\underset{\mathrm{rst}}{}{T}_{\mathrm{rst}3}^{\mathrm{nmk}3}\left(\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_j\right)e_{\mathrm{rst}3}\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)\right\}{N}_{\mathrm{nmk}}^H\left(\underset{\_}{r}-{\underset{\_}{r}}_j\right)\right]{\underset{\_}{a}}_j{v}_j& \left(23\right)\end{array}$

[0057] In the above equation the eigenvectors are expanded around the location of particle i in both the exciting field and the incident field. However, in the terms associated with the scattered fields, the eigenvectors are expanded around the location of particle j. It is worth noting as well that the latter eigenvectors are expressed in terms of Hankel functions while the former are expressed in terms of Bessel functions. In order to take subsequent advantage of the orthogonality of these vector functions, all of them need to be expressed in terms of the same type of functions and around a single location. For the purpose, I use the translational addition theorem for spherical vector wave functions well known in the prior literature which can be stated as: $\begin{array}{cc}\begin{array}{c}{L}_{\mathrm{nmk}}^H\left(\underset{\_}{r}-{\underset{\_}{r}}_j\right)=\underset{v=0}{\overset{}{}}\underset{=-v}{\overset{v}{}}{A}_v^{\mathrm{nmk}1}\left({\underset{\_}{r}}_j-{\underset{\_}{r}}_i\right)L_{vk}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_j\right)\\ =\underset{v}{}{A}_v^{\mathrm{nmk}1}\left({\underset{\_}{r}}_j-\underset{\_}{r_i}\right)L_{vk}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)\\ M_{\mathrm{nmk}}^H\left(\underset{\_}{r}-{\underset{\_}{r}}_j\right)=\underset{v}{}{A}_v^{\mathrm{nmk}2}\left({\underset{\_}{r}}_j-{\underset{\_}{r}}_i\right)M_{vk}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)+\\ A_v^{\mathrm{nmk}3}\left({\underset{\_}{r}}_j-{\underset{\_}{r}}_i\right)N_{vk}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)]\\ N_{\mathrm{nmk}}^H\left(\underset{\_}{r}-{\underset{\_}{r}}_j\right)=\underset{v}{}{A}_v^{\mathrm{nmk}3}\left({\underset{\_}{r}}_j-{\underset{\_}{r}}_i\right)M_{vk}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)+\\ A_v^{\mathrm{nmk}2}\left({\underset{\_}{r}}_j-{\underset{\_}{r}}_i\right)N_{vk}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)]\end{array}& \left(24\right)\end{array}$

[0058] Where the translational coefficients A_v^nmkq are given by: $\begin{array}{cc}{A}_v^{\mathrm{nmkq}}=\underset{p}{}{a}_q\left(n,m,v,,p\right)h_p\left(k_k{r}_{\mathrm{ij}}\right)P_p^{m-}\left(\cos {}_{\mathrm{ij}}\right)\exp (\left(m-\right){}_{\mathrm{ij}}& \left(25\right)\end{array}$

[0059] The index p varies over the range nv,nv2, . . . nv and the coefficients a_q(n,m,v,,p) are: $a_1\left(n,m,v,,p\right)={\left(-1\right)}^{}{}^{v+p-n}\left(2v+1\right)a\left(m,n-,vp\right)$ $a_2\left(n,m,v,,p\right)={\left(-1\right)}^{}{}^{v+p-n}a\left(n,v,p\right)a\left(m,n-,vp\right)$ $a_3\left(n,m,v,,p\right)={\left(-1\right)}^{}{}^{v+p-n}b\left(n,v,p\right)a\left(m,n-,vp,p-1\right)$ $a\left(n,v,p\right),b\left(n,v,p\right),a\left(m,n,-,vp\right)\mathrm{and}$ $a\left(m,n,-,vp,p-1\right)\mathrm{are}\text{:}$ $a\left(n,v,p\right)=\left[2v\left(v+1\right)\left(2v+1\right)+\left(v+1\right)\left(n-v+p+1\right)\left(n+v-p\right)-v\left(v-n+p+1\right)\left(n+v+p+2\right)\right]/\left[2v\left(v+1\right)\right]b\left(n,v,p\right)={\left[\left(n+v+p+1\right)\left(v-n+p\right)\left(n-v+p\right)\left(n+v-p+1\right)\right]}^{1/2}\left(2v+1\right)/\left[2v\left(v+1\right)\right]a\left(m,n,-,vp\right)={\left(-1\right)}^{m+}{\left(2p+1\right)\left[\frac{\left(n+m\right)!\left(v+\right)!\left(p-m-\right)!}{\left(n-m\right)!\left(v-\right)!\left(p+m+\right)!}\right]}^{1/2}\left(\frac{nvp}{m-\left(m+\right)}\right)\left(\frac{nvp}{000}\right)a\left(m,n,-,vp,q\right)={\left(-1\right)}^{m+}{\left(2p+1\right)\left[\frac{\left(n+m\right)!\left(v+\right)!\left(p-m-\right)!}{\left(n-m\right)!\left(v-\right)!\left(p+m+\right)!}\right]}^{1/2}\left(\frac{nvp}{m-\left(m+\right)}\right)\left(\frac{nvp}{000}\right)$

[0060] The symbol $\left(\frac{j_1{j}_2{j}_3}{m_1{m}_2-\left(m_1+m_2\right)}\right)$

[0061] is the well-known Wigner 3j symbol 6. Upon application of equations 24, equation 23 transforms to: $\begin{array}{cc}\underset{\mathrm{nmk}}{}\left[e_{\mathrm{nmk}1}\left({\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right)L_{\mathrm{nmk}}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)+e_{\mathrm{nmk}2}\left({\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right)M_{\mathrm{nmk}}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)+e_{\mathrm{nmk}3}\left({\underset{\_}{r}}_i;{\underset{\_}{a}}_i\right)N_{\mathrm{nmk}}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)\right]=\underset{\mathrm{nmk}}{}\left[i_{\mathrm{nmk}1}\left({\underset{\_}{r}}_i\right)L_{\mathrm{nmk}}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)+i_{\mathrm{nmk}2}\left({\underset{\_}{r}}_i\right)M_{\mathrm{nmk}}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)+i_{\mathrm{nmk}3}\left({\underset{\_}{r}}_i\right)N_{\mathrm{nmk}}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)\right]+{}_{\mathrm{Vi}}{}_{D}n\left({\underset{\_}{r}}_j,{\underset{\_}{a}}_j\right)g\left({\underset{\_}{r}}_j,{\underset{\_}{a}}_j,{\underset{\_}{r}}_i,{\underset{\_}{a}}_i\right)\underset{\mathrm{nmk}}{}\left[\left\{\underset{\mathrm{rst}}{}{T}_{\mathrm{rst}1}^{\mathrm{nmk}1}\left(\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_j\right)e_{\mathrm{rst}1}\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)\right\}\underset{v}{}{A}_v^{\mathrm{nmk}1}{L}_{vk}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)+\left\{\underset{\mathrm{rst}}{}{T}_{\mathrm{rst}2}^{\mathrm{nmk}2}\left(\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_j\right)e_{\mathrm{rst}2}\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)\right\}\underset{v}{}\left[A_v^{\mathrm{nmk}2}{M}_{vk}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)+A_v^{\mathrm{nmk}3}{N}_{vk}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)\right]+\left\{\underset{\mathrm{rst}}{}{T}_{\mathrm{rst}3}^{\mathrm{nmk}3}\left(\underset{\_}{o},\underset{\_}{m},{\underset{\_}{p}}_j\right)e_{\mathrm{rst}3}\left({\underset{\_}{r}}_j;{\underset{\_}{a}}_j\right)\right\}\underset{v}{}\left[A_v^{\mathrm{nmk}3}{M}_{vk}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)+A_v^{\mathrm{nmk}2}{N}_{vk}^B\left(\underset{\_}{r}-{\underset{\_}{r}}_i\right)\right]\right]{\underset{\_}{a}}_j{v}_j& \left(26\right)\end{array}$