Comparison of the Results of Modeling a Dispersed Medium by Wave and Ray Methods Sergey Pozdnyakov1, Sergey Ershov1, Nikolay Deryabin1 and Alexey Voloboy1 1 Keldysh Institute of Applied Mathematics RAS, Miusskaya sq. 4. Moscow, 125047, Russia Abstract A large number of works on dispersed medium modeling use either pure ray optics or light transport equation in which propagation of light obeys geometric optics while scattering properties of the medium can be either calculated with wave optics or measured. In either case the distance between individual particles must be much greater than wavelength. At the same time current computer power allows to simulate paint layer with wave optics. We decided to compare paint simulation done by the scalar wave approach and by ray tracing with individual particles. One of the goals of this work is to verify the correctness of ray tracing results for various sizes of metal flakes used often in production of metallic or pearlescent paints. Ray tracing had been done in two variants. One assumes the flakes have perfectly mirror reflection, while in the other variant the reflection is slightly diffuse with the angular distribution taken from the Fraunhofer diffraction on thin disk. For not too large flakes results of these two approaches substantially differ. The second “hybrid” method is considerably closer to the wave optics results. Keywords 1 Dispersed medium, paint simulation, wave optics, ray tracing 1. Introduction Many works on modeling various coloring substances have been published in the literature on computer graphics [1-4]. Most them are based either on the pure geometric ray tracing, or use the various variants of the “continuous medium approximation” [5-7] also known as the “light transport equation” (LTE) or “radiation transfer”. This latter is rather a phenomenological equation which can sometimes work even for sub- wavelength size particles. But it requires that we know scattering properties of the medium which must be calculated elsewhere or measured. Usually they are taken as just the scattering properties of an individual isolated particle scaled proportionally to the volumetric fraction of such particles in the paint medium. Scattering of a single isolated particle can be calculated by wave optics (for example, by the Mie theory or its analogs). This method works quite well when the distance between particles is much greater that wavelength. Exactly the same limitation emerges in the pure ray tracing method. In case of small particles in high concentrations this condition is violated rendering the simulation inaccurate. At the same time, current computer power allows to simulate paint layer by wave optics (solving diffraction problem for the ensemble of particles). The exact equations describing the propagation of light are the Maxwell equations [8] whose unpleasant feature is that normal components of fields can be discontinuous across the boundaries between media. This requires a special grid and finite-difference GraphiCon 2021: 31st International Conference on Computer Graphics and Vision, September 27-30, 2021, Nizhny Novgorod, Russia EMAIL: mephi32@rambler.ru (S. Pozdnyakov); ersh@gin.keldysh.ru (S. Ershov); dek@keldysh.ru (N. Deryabin); voloboy@gin.keldysh.ru (A. Voloboy) ORCID: 0000-0002-5493-1076 (S. Ershov); 0000-0003-1248-6047 (N. Deryabin); 0000-0003-1252-8294 (A. Voloboy) ©️ 2021 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). CEUR Workshop Proceedings (CEUR-WS.org) approximation to treat those discontinuities. Applying these to a piece of medium with randomly distributed particles requires very large computation time and RAM. So instead of the exact electromagnetic calculations for such media various simplifications are often used, e.g. the boundaries between media (i.e. the step function in dielectric constant) are replaced with a continuous transition. The error introduced is unknown in advance (and often after calculations). In case when polarization effects can be neglected, the vector wave from the Maxwell equations can be replaced with the scalar one. Also the stationary diffraction problem (when the fields are monochrome and so the time can be excluded) is enough. This leads to the Helmholtz equation [8]. This approximation is widely used in optics when calculating the propagation of unpolarized radiation. Very roughly, the scalar diffraction corresponds to electromagnetic waves polarized in the plane of separation of different media because for normal incidence solutions of the Maxwell and Helmholtz equations coincide. For a disperse medium with randomly distributed particles in high concentration scattering effectively depolarizes the radiation and so the use of the scalar approximation is quite justified. Let us repeat that the scalar approximation allows to observe at least the qualitative difference between the ray and wave optics and frequently it is quite quantitatively accurate. Some program for scalar calculations in plane geometry had been created at the Keldysh Institute of Applied Mathematics. A fairly large number of calculations were carried out for various types of dispersed medium with its help. A detailed description of the physical and mathematical approach and calculation results can be found in [9]. The model of paint layer used in the above calculations consists of the embedding medium (binder) and particles dispersed in it which can have random size, shape and orientation. The media of the paint layer elements can have an arbitrary refraction index which allows to simulate both dielectrics or metals. For the latter case the imaginary part of refraction dominates. The goal of our study is to compare the result of dispersed medium simulation done by ray optics and wave optics methods. For example, we were interested how the difference between them changes with size and concentration of the metal flakes. Also we would like to find an approach which allows to improve ray simulation results making them closer to wave simulation ones. In this work we consider a paint layer consisting of thin metal disks (flakes) embedded in a dielectric medium (varnish or binder). These flakes are distributed statistically uniformly in the layer. Their orientation is random: the normal to the flake surface slightly deviates from the normal to the whole layer and this deviation is random. Paints of similar structure are called “metallic” or “pearlescent” and are widely used in automotive industry. Simulation with LTE requires the phase function of the whole material which, as explained above, can be approximated by the scaled indicatrix of a single pigment particle, usually calculated from light diffraction on that isolated particle. For rather large planar flakes one can use various approximate diffraction theories, e.g. the Kirchhoff approximation. It gives a quite accurate angular distribution in the far zone. If the characteristic distances between the plates are small then the use of distributions calculated in this way can lead to significant uncontrollable errors so we used Fraunhoffer approximation (most accurate in the near zone) instead. The wave optics was simulated using the scalar diffraction theory. 2. Paint structure In our study the structure of the paint layer is stochastic and statistically uniform. Our scalar diffraction solver assumes periodic boundary conditions on the side boundaries of the domain. This is equivalent to an infinite plane parallel layer made of periodically tiled “pieces”. Such an approach is widely used in diffraction problems [10]. A detailed description of the geometry used in the calculations is contained in [9]. As mentioned above, the layer consists of small metal flakes embedded in the binder. These flakes are spaced within the layer evenly on average. The flakes are not parallel to the surface of the layer but are randomly inclined. The distribution of the angle between the normal to the plate surfaces and the normal to the layer surface is Gaussian. As a rule, real varnishes are good dielectrics with low absorption. They have minimal dispersion and their refractive index is close to 1.5 and we set refraction to exactly 1.5. Silver was chosen as the material of the plates as the best reflector. Refractive indices of silver at various wavelengths were taken from [11]. Real paints use aluminum or zinc. The choice of silver for simulation was dictated by convenience: its refractive index (a small real part and not too large imaginary) allows using grids with not too many vertices without loss of accuracy. In finite-difference wave calculations the size of the grid step is often chosen somewhere between 𝑠 = λ/20 and 𝑠 = λ/10, where λ is the wavelength in the corresponding medium. With the help of reference calculations it was found that for silver it is possible to use horizontal step of 40 nm (notice illumination is directed vertically) and vertical step of 10 nm. For other metals with greater values of the imaginary and real parts of the refractive index the grid steps should be smaller which requires more computer RAM and more time-consuming calculations. So silver is preferable. At each grid vertex the refractive index of the corresponding material is set. All metal flakes were the same in each calculation. They were thin disks of thickness 0.3 μm and diameter ranging from 5 μm to 10 μm. These dimensions are close to that of the really used flakes. To create sample geometry of the paint layer one needs to distribute in it without intersections the flakes whose number is determined by their concentration. Because of periodicity of the boundary conditions at the side walls, the layer consists of periodically tiled elementary fragments (whose horizontal size equals one period). In our calculations it was about 41 μm. The layer thickness was 5μm that is close to the real paint values. Horizontal dimensions of the elementary fragment are mainly determined by the maximum diameter of the metal plates (so that there be enough flakes inside). One of the most important paint characteristic is its hiding power. In our case with opaque flakes the hiding power depends on the total surface area of the flakes per unit area of the layer surface. By area of one flake we mean the area of one basement of its disk. It is easy to see that if the thickness of the flake is the same for all variants then their total area per unit surface area of the layer is proportional to the pigment volume concentration (PVC). For the paint model under consideration the hiding power of the paint becomes quite good even at not too high PVC. For example, for PVC = 30% the ratio of the total surface area of flakes to the corresponding surface area of the layer is about 5. According to our experience in stochastic disperse media the simplest method to generate ensemble of non-intersecting flakes can be used, namely: a simple random selection of the slopes and coordinates of the centers of the disks satisfying the corresponding distributions. In practice, the calculation procedure for the position and orientation of the given disk is quite simple. The polar angle of the flake’s normal is chosen at random according to the formula [13] which gives Gaussian random values: 𝜃 = 𝜃0 |cos(2𝜋𝜉1 ) √−ln 𝜉2 | where 𝜉1 and 𝜉2 are random variables uniformly distributed over the interval (0, 1) and 𝜃0 was 7.09°. In this case the average angle of inclination of the normal is about 4 degrees. The azimuth is uniformly distributed in [0, 2𝜋]. This procedure produces a horizontally isotropic distribution of orientation which well approximates the real paint structure. After that, coordinates of the center of the particle are chosen at random. The horizontal coordinate is chosen uniformly in the interval [0, 41] microns (one period). The vertical coordinate z is also chosen uniformly so that the inclined disk does not cross the layer boundaries. Note that the coordinates in the plane of the layer are chosen along the entire length and width of the elementary fragment in order not to introduce non-uniformity in the distribution of particles near these artificial walls. Thus flakes near the side walls can “penetrate” them and because of periodicity the protruded part of that flake is near the opposite size of the domain, see Figure 1. 1´ 1 2 2´ Elementary fragment (one period) Figure 1: Flakes in a periodically tiled elementary fragment. The binder is shaded violet, flakes are thick green dashes. The center of flake 1 is near the right boundary of the elementary fragment, so repeated periodically as 1´ a part of it is near the left boundary. Similarly for the flake 2 whose center is near the left boundary, a part of it appears near the left boundary as flake 2´. Flakes are added to the elementary fragment one by one until the desired number of them (determined by concentration) is achieved. While adding the next particle, its intersection with the existing ones and with the top and bottom domain boundaries is checked. If there are no intersections, the flake is added. Otherwise we generate its position and orientation anew and repeat the procedure like in the rejection sampling method. Notice that because of periodicity we must check intersection with flakes from the 8 adjacent periodic fragments (see Figure 1). 3. Wave optics simulation As mentioned above, an accurate calculation of light propagation in a complex medium should be based on the Maxwell equations. Although the scalar wave approximation by the Helmholtz equation is quite applicable in case when the incident wave is nearly normal to the boundary between media. For normal or close to it illumination the situation is exactly this: the paint layer boundaries are exactly horizontal while the flakes which are not tilted much are also close to horizontal. Below we shall outline the calculation scheme from [9] and mention the features related to the present case. The Helmholtz equation for the stationary diffraction problem operates a single time-independent scalar field being the amplitude of the monochrome wave. The paint layer is illuminated by a monochrome plane normal to its surface wave and the result of the simulation is the complex wave amplitude in the elementary fragment. The boundary conditions at the side walls of the elementary domain are periodic. At the top and bottom boundary we use the well-known Sommerfield radiation conditions: the wave field and its derivative smoothly conjugate with the field outside the domain which is:  a sum of the incident wave and waves travelling upwards in the upper half-space;  a sum of waves travelling downwards in the lower half-space. which gives linear relations between the field and its normal derivative that the top and the bottom boundaries, see [9]. These boundary conditions are common for most of solutions of the stationary diffraction problems because they allow to operate the wave field only inside the domain of interest and not to keep in memory its values in the uniform space outside it. In case of need the wave field in that uniform space can be analytically calculated from the boundary values using an analog of the Green’s formula. For the normal illumination the wave field is simply periodic in the horizontal coordinates x and y. For a general incidence, the wave field is a periodic function times 𝑒 𝑖𝑘𝑥 𝑥+𝑖𝑘𝑦 𝑦 where k is the wave vector of the incident wave. This is named the Bloch representation. So we can rewrite our Helmholtz equation so that it operates a periodic field. For a periodic function, one of the main methods of analysis is the representation in the form of a Fourier series. Formally the Fourier series contains an infinite number of terms and for a numerical simulation we must somehow limit their number, e.g. as in the Galerkin method. In case the field is operated on a space grid we just use the Discrete Fourier Transform (DFT) which performs a one-to- one transformation from the function to its Fourier image and back again and thus the number of the Fourier modes just equals the number of the grid vertices. If the grid contains a power of 2 vertices per period the DFT is especially efficient and is termed Fast Fourier Transform (FFT) [13]. The 2D FFT is just two one-dimensional transforms applied one after another. Applying the FFT in x, y the original Helmholtz equation becomes a system of the ordinary differential equations in the vertical coordinate z for all the harmonics. In case the medium is inhomogeneous (refraction index is space dependent) these equations are coupled. To solve this system a special procedure similar to the relaxation method had been developed. It introduces an artificial time (as it was mentioned above we solve a stationary problem which does not have genuine time!) changing the Helmholtz equation into an analog of the non-stationary Schrödinger equation. For this latter we begin with some first approximation of the solution and trace its evolution in this artificial time. It can be proved that eventually it converges to the stationary state i.e. solution to the Helmholtz equation in the non-uniform medium [9]. 4. Geometric optics (ray tracing) and its modification Light scattering in a paint layer can be calculated with the Monte-Carlo ray tracing. It uses sample geometry where the individual flakes are positioned so as to reproduce the statistical parameters of their distribution in the paint: concentration and orientation. Since the same geometry is used for wave calculations it is a periodically tiled elementary fragment (see above). When tracing ray in such a geometry one may also use only one elementary fragment but then a ray when hitting its side wall “jumps” to the opposite side (because of periodicity). This must be duly taken into account when searching for the intersection with flakes (see also Figure 1). In pure ray optics simulation the ray is specularly reflected (or absorbed) when it hits a flake. Meanwhile, for the flakes of size about 5 or 10 microns the incident wave would diffract on it and the outgoing light will be not a single parallel wave (an analog of the mirror reflection in ray optics) but a bunch of waves with slightly different directions. Interaction of a plane wave with a thin disk can be calculated analytically, see e.g. [12]. So we can “modify” ray tracing procedure by making the ray reflection by a flake not ideally specular but somehow use the light distribution from diffraction on a disk. Since wave optics does not operate “hit point”, it is extremely difficult to decide how reflection of a ray by the flake surface depends on the hit point position (e.g. if it is near the center, scattering is weaker than when it is near the edge). The simplest solution is to attribute to the flake surface a diffuse BRDF which is the same in all points and provides angular distribution of reflected light as in the diffraction of a plane wave on that disk. We did try this approach. This is something like merging the phase function of particle taken from the wave optics (as in LTE), pure ray tracing and wave optics (to get that phase function). Below this method is termed “ray optics modified”. The angular distribution was taken from the solution to the diffraction on a thin disk in the Fraunhofer approximation [12]. Notice that the angular spread of the reflected beam is about λ⁄D, where λ is the wavelength of the light wave in the binder and D is the flake’s diameter. For D=10μm and wavelength 550 nm in vacuum it is about 2 degrees. This value corresponds to small angles at which a large difference is already observed between the results of the wave and ray calculations. 5. Results The model paint layer has thickness 5 micron. Refraction index of the binder is 1.5 (at all wavelengths). Flakes are thin silver disks of thickness 300 nm and diameter ranging from 5 to 10 microns, distributed uniformly with volumetric concentration PVC=30%. Reflection of flake surface for ray calculations was 90%. For the sake of simplicity we neglected the effect of the Fresnel boundaries at the top and bottom of that layer, i.e. assumed that above and below there is a clear binder without particles. The main goal of this work was to compare the ray and wave approaches. So we did simulations for normal incidence only. Calculations were done for three wavelengths: 450nm, 550nm and 650nm (in vacuum) and three values of flake diameter: 5 μm, 7 μm, and 10 μm. These values were chosen so that the area of the surface of a flake increases by a factor of two. The number of flakes set for the elementary fragment to keep the PVC was 427, 218, and 107 respectively for the above 3 variants of the diameter. We employed 3 methods: wave optics (scalar stationary diffraction), pure ray optics (flake reflect specularly) and “modified ray optics” (see Section 4) when flake surface is sharp diffuse. All the three methods were applied to the same sample geometry which was horizontally periodic. So first we created geometry of the elementary fragment as described in Section 2. Horizontal period was 41 micron which is mainly determined by the maximum diameter of flakes. Geometry sample created this way is random and we must average scattering calculated for several random samples of geometry and average, the more the better. We decided to use exactly the same geometry with all the three methods. Therefore the number of samples had been limited by the speed of the slowest one, i.e. the wave calculations. So we used a limited number of geometry instances, namely just about a dozen. In wave calculations the grid step is often chosen somewhere between 𝑠 = λ/20 and 𝑠 = λ/10, where λ is the wavelength in the corresponding medium. With the help of reference calculations it was found that for silver it is possible to use a horizontal step of 40 nm and a vertical step of 10 nm. For other metals with greater values of the imaginary and real parts of the refractive index the grid steps should be smaller which would require more computer RAM and more time-consuming calculations At each grid vertex the refractive index is determined depending on whether this point belongs to a flake or not; then refraction is that of the flake material or of the binder, respectively. The calculation results are presented in Figures 2, 3 and 4 which show paint’s BRDF in luminance factor units for different wavelength and flake size, calculated with the 3 methods. One can see that while the results of the wave calculations and “modified ray optics” are rather smooth curves, those of the pure geometric optics exhibit strong oscillations in the near-specular region. This is because of the small number of samples of the random elementary fragment (limited by the speed of wave optics). For wave optics and “modified ray optics” this number of samples was sufficient because in these methods a flake diffuses the light incident on it and this is, so to say, additional smoothing. In case of pure geometric optics when flake reflection is pure specular we face a strong “sparkling effect”. The “genuine” BRDF obtained as the average over an infinite number of geometry samples would naturally be smooth too but a finite number of samples happens insufficient to kill the remaining fluctuations. Notice that the BRDF obtained by pure ray tracing has a strong near-specular peak absent in the other approaches. Indeed, flakes are rather well aligned so among them usually one or two are nearly exactly horizontal. Reflection from them makes a sharp and high specular peak of considerable amplitude (in case of diameter 10 microns one flake occupies about 5% of the elementary fragment surface area, thus for 2 horizontal flakes the specular peak contains about 10% of the total reflected energy). 500 600 550nm D=5micron 550nm D=7micron 500 400 Wave optics 400 Wave optics 300 Geom. optics Geom. optics 300 200 Geom. optics modified Geom. optics modified 200 100 100 Theta Theta 0 0 0 1 2 3 4 5 6 7 8 9 10 Figure 2: BRDF (in luminance factor units) at wavelength 550 nm in vacuum for flake diameter 5μm (left) and 7μm (right). θ is the angle between the view direction and mirror reflection of the incidence by a horizontal plane (so in our case between the view direction and the vertical). 800 800 450nm Diameter=10 micron 550nm D=10micron 700 700 600 Wave optics 600 Wave optics 500 Geom. optics 500 Geom. optics 400 Geom. optics modified 400 Geom. optics modified 300 300 200 200 100 100 Theta Theta 0 0 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10 Figure 3: BRDF (in luminance factor units) at wavelength in vacuum 450 nm (left) and 550 nm (right) for flake diameter 10 μm. θ is the angle between the view direction and mirror reflection of the incidence by a horizontal plane (so in our case between the view direction and the vertical). 800 800 550nm D=10micron 650nm D=10micron 700 700 600 600 Wave optics 500 500 Wave optics Geom. optics 400 Geom. optics modified 400 Geom. optics 300 300 Geom. optics modified 200 200 100 100 Theta Theta 0 0 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10 Figure 4: BRDF (in luminance factor units) at wavelength in vacuum 550 nm (left) and 650 nm (right) for flake diameter 10 μm. θ is the angle between the view direction and mirror reflection of the incidence by a horizontal plane (so in our case between the view direction and the vertical). The usually applied criterion for “applicability of ray optics” is that all dimensions (size of particles and distance between them) are large compared to wavelength. In case of paint layer we are “on the eve” for both criterion even for the largest flakes. So the pure ray optics can be formally inaccurate and this is so indeed. But our investigation shows details about that inaccuracy. We see that agreement is quite good in the off-specular zone. For diameter 5 microns we also see that the “modified ray optics” is more or less applicable even in the near-specular region. But for larger flakes of diameter 10 microns it is not so (which is strange). There is practically no dependence on the illumination wavelength in contrast to the dependence on the particle size. For the wave approach this situation is probably due to the rather strong dispersion of silver whose refraction index is 0.13547+2.3808i, 0.14512+3.19i and 0.15943+3.9291i for wavelengths 450nm, 550nm and 650nm (in vacuum) respectively. It is hardly possible to expect the usual weakening of wave effects with decreasing radiation wavelength for such a strong dispersion. For the ray approach this dependence should be very weak. The only effect the wavelength has is only through the same refractive index on which the reflection coefficient from the metal surface depends. The difference in reflectance at the wavelengths used is small and therefore the distributions of reflected light for ray tracing are similar. 6. Discussions and conclusion The purpose of this work was to compare the results of ray and wave calculations of light reflection from a paint layer consisting of small flat metal flakes dispersed in a dielectric binder. Wave calculations were done for the stationary scalar diffraction problem. Applicability of the scalar diffraction theory for small angles on incidence (and observation) is quite obvious. Ray optics calculations were done both for pure ray optics when flakes have specular reflection and for the “modified approach” when they have diffuse BRDF taken from diffraction on thin disk. For our model paint the flake diameter is considerably greater than wavelength in the binder. The (vertical) distance between flakes is smaller but all the same exceeds wavelength. Therefore ray optics can be accurate in our case; but one cannot state that definitely because the criterion “must be substantially greater than the wavelength” is too vague. Our simulations shown that ray optics does give a good approximation in the off-specular region of BRDF while in the near-specular region it strongly deviates from the wave results. The “modified ray optics” results are closer to the wave one but still deviate in the near-specular zone. The angular distribution of reflected light for both wave optics and the “modified ray optics” is quite smooth while for the pure ray optics it has a strong near-specular peak. This is because in the two former approaches the light reflected by a flake is “diffused” because of diffraction effects, accurate in wave optics and approximate in the “modified ray optics”. The difference between the wave and ray calculations means that the diffraction effects are substantial even for flakes of rather large size (10 microns). Therefore the pure ray tracing is of limited applicability in such cases. Even the “modified ray optics” does deviate from the full diffraction calculations. Therefore the difference between the ray and wave results is not only due to the angular “blurring” of the light reflected from a flake but also due to interference between flakes which the “modified ray optics” neglects. 7. Acknowledgements The work is supported by RFBR, grant 20-01-00547. 8. References [1] S. Ershov, K. Kolchin, K. Myszkowski, Rendering pearlescent appearance based on paint- composition modeling, Computer Graphics Forum 20(3) (2001) 227-238. [2] M. Rump, G. Muller, R. Sarlette, D. Koch, R. Klein, Photo-realistic Rendering of Metallic Car Paint from Image-Based Measurements, Comput. Graph. Forum 27(2) (2008) 527-536. doi:10.1111/j.1467-8659.2008.01150.x [3] R. Durikovic, K. Kolchin, S. Ershov, Rendering of Japanese artcraft, in: Proceedings of the EUROGRAPHICS, short presentations, 2002, pp. 131-138. [4] J. Filip, R. Vavra, Image-based appearance acquisition of effect coatings, Computational Visual Media 5(1) (2019) 73-89. doi:10.1007/s41095-019-0134-3 [5] A.G. Voloboy, S.V. Ershov, S.G. Pozdnyakov, Interactive modeling of automotive paints, in: Proceedings of the 22nd International Conference on Computer Graphics and Vision, 2012, pp. 242-247. [6] S. Ergun, S. Onel, A. Ozturk, A general micro-flake model for predicting the appearance of car paint, in: Proc. of the Eurographics Symposium on Rendering: Experimental Ideas & Implementations, 2016, pp. 65-71. doi:10.2312/sre.20161211 [7] G.Y. Kim, K.H. Lee, A Reflectance Model for Metallic Paints Using a Two-Layer Structure Surface with Microfacet Distributions, IEICE Transactions on Information and Systems, E93-D (11), (2010) 2076-3087. [8] N. S. Koshlyakov, E. B. Gliner, M. M. Smirnov, Equations in partial derivatives of mathematical physics, Higher school, Moscow, 1970 [9] A.G. Voloboy, S.V. Ershov, S.G. Pozdnyakov, Ink layer appearance simulation, Matematicheskoe modelirovanie, 26(2) (2014) 81-94 (in Russian). [10] Akira Ishimaru, Wave Propagation and Scattering in Random Media, Academic press, New York, San Francisco, London, 1978. [11] Refractive index database, 2021. URL: https://refractiveindex.info/?shelf=3d&book=metals&page=silverS. [12] Airy disk, 2021. URL: https://en.wikipedia.org/wiki/Airy_disk [13] W. Press, S. Teukolsky, W. Vetterling, B. Flannery. Numerical Recipes in C. The Art of Scientific Computing, Second Edition, Cambridge University Press, 1992.