=Paper= {{Paper |id=Vol-3027/paper12 |storemode=property |title=Comparison of the Results of Modeling a Dispersed Medium by Wave and Ray Methods |pdfUrl=https://ceur-ws.org/Vol-3027/paper12.pdf |volume=Vol-3027 |authors=Sergey Pozdnyakov,Sergey Ershov,Nikolay Deryabin,Alexey Voloboy }} ==Comparison of the Results of Modeling a Dispersed Medium by Wave and Ray Methods== https://ceur-ws.org/Vol-3027/paper12.pdf
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.