Beam-splitting code for light scattering by ice crystal particles within geometric-optics approximation

Journal of Quantitative Spectroscopy & Radiative Transfer 164 (2015) 175–183

Contents lists available at ScienceDirect

Journal of Quantitative Spectroscopy & Radiative Transfer journal homepage: www.elsevier.com/locate/jqsrt

Beam-splitting code for light scattering by ice crystal particles within geometric-optics approximation Alexander V. Konoshonkin a,b,n, Natalia V. Kustova b, Anatoli G. Borovoi a,b a b

National Research Tomsk State University, Lenina Street 36, 634050 Tomsk, Russia V. E. Zuev Institute of Atmospheric Optics SB RAS, Academician Zuev Square 1, 634021 Tomsk, Russia

a r t i c l e i n f o

abstract

Article history: Received 6 April 2015 Accepted 8 June 2015 Available online 17 June 2015

The open-source beam-splitting code is described which implements the geometricoptics approximation to light scattering by convex faceted particles. This code is written in C þ þ as a library which can be easy applied to a particular light scattering problem. The code uses only standard components, that makes it to be a cross-platform solution and provides its compatibility to popular Integrated Development Environments (IDE's). The included example of solving the light scattering by a randomly oriented ice crystal is written using Qt 5.1, consequently it is a cross-platform solution, too. Both physical and computational aspects of the beam-splitting algorithm are discussed. Computational speed of the beam-splitting code is obviously higher compared to the conventional raytracing codes. A comparison of the phase matrix as computed by our code with the raytracing code by A. Macke shows excellent agreement. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Beam-splitting code Light scattering Geometric-optics approximation Ice crystals

1. Introduction The problem of light scattering by nonspherical and large, as compared with incident wavelengths, particles is a challenging problem for the light-scattering community. The conventional methods solving directly the Maxwell equations like the T-matrix method [1], the finitedifference time-domain (FDTD) method [2], and the dipole-discrete approximation (DDA) [3] are capable to solve this light scattering problem up to the sizes of, say, 20 μm for the visible. For the larger particles, computer resources restrict further simulations. In this case, the geometric-optics approach looks as an obvious and reliable alternative to the abovementioned methods. Usually geometric-optics solutions to the scattering problem are obtained by use of ray-tracing techniques. For example, this is the ray-tracing code developed by Macke [4] that is

n

Corresponding author. Tel.: þ 7 962 777 26 02. E-mail address: [email protected] (A.V. Konoshonkin).

http://dx.doi.org/10.1016/j.jqsrt.2015.06.008 0022-4073/& 2015 Elsevier Ltd. All rights reserved.

freely available to the scientific community and therefore it is widespread. In the problem of light scattering by large nonspherical particles, there is a specific case of ice crystal particles of cirrus clouds where the crystal sizes range from a few micrometers to millimeters. For such faceted particles, the conventional ray-tracing procedures become mathematically excessive. Indeed, there are a lot of photon trajectories inside the crystals that are parallel to each other. All these trajectories bear the same information except of their perpendicular shifts. A set of these similar trajectories create a plane-parallel light beam of some crosssection and shape. Therefore the propagation of light inside crystal particles with the reflection/refraction events on faceted surfaces can be compactly simulated as the propagation and reflection/refraction of the whole plane-parallel beams which are initiated from illuminated crystal facets. Polarization characteristics of any ray inside the beam are the same. Such an alternative to ray-tracing procedures were discussed and independently explored by a number of

176

A.V. Konoshonkin et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 164 (2015) 175–183

authors. In particular, Popov was the first who discussed such an approach [5]. Then del Guasta [6] developed a code named the facet tracing and applied it to calculate the backscatter by hexagonal ice crystals. Also the analogous codes were developed and used by Romashov [7], by Borovoi and Grishin [8] and by Borovoi et al. [9]. Recently, the same beam-splitting code was considered by Bi et al. [10] where the case of absorbing crystals was included as well. However, none of these codes are open-source and freely available. The aim of this paper is to present our open-source beam-splitting code [11] that has some advantages as compared to the conventional ray-tracing techniques in the case of ice crystals. The strengths and weaknesses of this code are discussed.

2. Theoretical basis 2.1. The scattered field in the near zone It is obvious that, within the framework of geometric optics, the scattering of an incident plane electromagnetic wave E0 is reduced to a sequence of reflection/refraction events produced by crystal facets. The incident planeparallel wave is split by the crystal into a set of planeparallel beams (see Fig. 1). As a result, the scattered field Es on the crystal surface and in the near zone becomes a superposition of the beams Ej leaving the crystal surface in different directions. The incident plane electromagnetic wave with the propagation direction i has the following view: 0 E0 ðrÞ ¼ E expðikir þ iψ 0 Þ ¼ @ 0

E0jj E0?

1 Aexpðikir þiψ Þ; 0

ð1Þ

where E0 is a constant polarization vector [12] and ψ 0 is a constant phase shift. Then the scattered field Es outside the crystal is the superposition of the plane-parallel beams Ej Es ðrÞ ¼

X

Ej ðrÞ:

ð2Þ

Every beam is described as

0

Ej ðrÞ ¼ E ηj ðrÞexpðiknj r þiψ j Þ ¼ @ j

Ejjj Ej?

1 Aη ðrÞexpðiknj r þ iψ Þ; j j ð3Þ

where the additional factor ηj ðrÞ determining the beam shape appears ( 1 inside the j  th beam; ηj ðrÞ ¼ ð4Þ 0 otherwise; and the constant polarization vector Ej perpendicular to the propagation direction nj is the product Ej ¼ Lm þ 1 Fm Lm :::F2 L2 F1 L1 E0 :

ð5Þ

In Eq. (5), every refraction/reflection event changes the polarization vector by multiplication on the Fresnel matrix F. This matrix is diagonal in the appropriate coordinate system (see Fig. 5), i.e. ! F jj 0 ; ð6Þ F¼ 0 F? where F jj and F ? are the Fresnel reflection/refraction coefficients for a plane interface [13]. The rotation matrix   sin α cos α L¼ ð7Þ  sin α cos α adjusts these coordinate systems, where α is an angle between previous and current coordinate systems. In Eq. (5), we have m refraction/reflection events. The last rotation matrix Lm þ 1 in Eq. (3) transforms the outgoing beam to a final coordinate system. The basis vectors determining the parallel and perpendicular components of the constant polarization vector Ej are defined in Section 2.3. Instead of the polarization vector Ej , it is more convenient to use the ð2  2Þ Jones matrix Jj connecting the polarization vectors Ej and E0 as follows: Ej ¼ Jj E0 :

ð8Þ

Thus, the Jones matrix of the j-th beam is equal to Jj ¼ Lm þ 1 Fm :::F2 L2 F1 L1 :

ð9Þ

Fig. 1. Example of the plane-parallel beams split by a crystal. (a) Light is incident normally from a reader. One of the illuminated facets is selected (colored). (b) and (c) illustrates the beam-splitting initiated by the colored facet. (b) shows only three beams, for simplicity. The shapes of two transmitted beams are shown in c. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

A.V. Konoshonkin et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 164 (2015) 175–183

177

the Jones matrix of Eq. (9) according to the general rule M ¼ ΓðJ  Jn ÞΓ

1

;

ð12Þ

where 0

1

0

0

1 B1 2@ 0

0

0

1

1

0

i

i

Γ ¼ pffiffiffiB B

1

1

1 C C C: 0A

ð13Þ

0

Any scattering to the far zone is usually described by the Mueller (scattering) matrix Mðn; iÞ connecting the Stokes parameters of the far-zone-scattered field Is ðnÞ with the Stokes parameters of the incident wave I0 as Is ðnÞ ¼ Mðn; iÞI0 :

Fig. 2. Example of the outgoing beams mapped into the scattering direction hemisphere.

2.2. The scattered field in the far zone The scattered field is always observed at large distance from the crystal particles R-1 where the scattered field has a view of a diverging spherical wave. Consequently, any beam of Eq. (6) being transformed into a diverging spherical wave in the far zone creates a diffraction pattern which is distributed over the scattering direction sphere n. If a transverse size of a beam in the near zone is much larger than the wavelength b 4 4 λ, the diffraction pattern in the far zone is narrow, it occupies a cone of the angle Δθ  λ=bo o 1 about the scattering direction n ¼ nj . The scattered field is often detected with rough angular resolution that is larger than Δθ  λ=b. In this case, the diffraction pattern can be effectively replaced by the Dirac δ-function δðn  nj Þ (Fig. 2). In the geometric-optics approximation, scattering of light is consistently formulated in terms of the Stokes parameters as we discussed earlier in details [14,15]. We define the Stokes parameters as 0 1 0 1 E2jj þ E2? I B C B Q C B E2  E 2 C C B C B jj ? I¼B C¼B ð10Þ C: @ U A B  2ReEjj En C ? A @ V 2ImEjj En?

ð14Þ

This matrix is dependent on the coordinate systems assumed. Majority of authors refer the coordinate systems to the scattering plane, i.e. the plane containing the incident i and scattering n directions (Fig. 3). The matrix in these coordinates is called the phase matrix P(n, i) [17]. Here ~ ϕ ÞMðn; iÞLð ~ ϕ Þ; Pðn; iÞ ¼ Lð 2 1 where   L~ ϕ ¼ ΓðL  LÞΓ

1

0

1

B0 B ¼B @0 0

ð15Þ

0

0

cos 2ϕ

sin 2ϕ

 sin 2ϕ

cos 2ϕ

0

0

0

1

0C C C: 0A 1

It is worthwhile to mention that all symmetry properties of the matrixes discussed in [16–19] are inherent only to the phase matrix. Eq. (11) results in the following view of the total Mueller matrix where all beams are added incoherently X ð17Þ Mðn; iÞ ¼ sj δðn nj Þmj : This matrix obeys the energy conservation law since every reflection/refraction event inside a crystal satisfies this law. For the incident wave of the unit intensity and arbitrary polarization, we obtain the obvious equations for the differential cross section σ ði; nÞ and total cross section

This definition is the same as in [16,17] but the signs of U and V are opposite in [18,19]. This difference is explained by opposite directions of the basis perpendicular vectors e ? ¼ e0? assumed. The j-th beam of Eq. (3) in the far zone is described by means of the Stokes parameters as follows: Ij ðnÞ ¼ sj δðn nj Þmj I0 ;

ð11Þ

where the column-vectors of the scattered Ij and incident I0 fields are formed from the polarization vectors Ej and E0 , respectively, according to Eq. (10). The transverse area of the beam sj is the integral of the factor ηj ðrÞ of Eq. (4) over any plane perpendicular to the propagation direction nj . The matrix mj is the ð4  4Þ Mueller matrix obtained from

ð16Þ

Fig. 3. Relations between the Mueller and phase matrices.

178

A.V. Konoshonkin et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 164 (2015) 175–183

σ ðiÞ as σ ði; nÞ ¼ σ ðiÞ ¼

Z

X

sj δðn nj ÞI j ;

σ ði; nÞdn ¼

X

ð18Þ

and two perpendicular unit vectors are directed parallel esjj and perpendicular es? to the meridional plane. We define 8 θs ¼ 0; > < ð0; 1; 0Þ; es? ¼

sj I j ;

ð19Þ

ð  sin φs ; cos φs ; 0Þ; > : ð0; 1; 0Þ;

0 o θs o π ;

ð23Þ

θs ¼ π :

Then the vector esjj is found as where I j is the intensity of the j-th beam defined by Eqs. (3)–(5). For the non-absorptive particles, the scattering crosssection is equal to area of the particle shadow S for the given incident direction and a fixed particle orientation

σ ðiÞ ¼ S;

ð20Þ

remembering that the intensity of the incident field is I 0 ¼ 1. Also we notice a useful analytical equation for the averaged shadow area 〈S〉 that is valid for any convex particles in the case of their random orientation 〈S〉 ¼ B=4;

ð21Þ

esjj ¼ es?  n;

ð24Þ (esjj ,es? ,

that corresponds to the right-hand triple n) of the basis vectors. The fourth coordinate system is similarly constructed for the incident wave where the incident propagation direction i is defined as i ¼ ð  sin θi cos φi ;  sin θi sin φi ;  cos θi Þ: Here the right-hand basis vector triple used where 8 θi ¼ 0; > < ð0; 1; 0Þ; ei? ¼ ð  sin φi ; cos φi ; 0Þ; 0 o θi o π ; > : ð0; 1; 0Þ; θ ¼ π;

ð25Þ i) is also

ð26Þ

i

and

where B is the total area of the particle surface.

eijj ¼ ei?  i: 2.3. The coordinate systems It is important to specify clearly the coordinate systems used because of a lack of uniformity in the literature yet. We use four coordinate systems shown in Fig. 4. The first is the right-hand Cartesian system (x, y, z) referred to as the laboratory coordinate system. Its origin usually coincides with the center of a scattering particle. Then the azimuth θ and zenith φ angles of the conventional spherical coordinate system determine the coordinates (θ, φ). The third coordinate system defines a righthand triple of unit vectors on the scattering direction sphere that describes polarization of the scattered field Es : The propagation direction of the scattered field is defined as n ¼ ð sin θs cos φs ; sin θs sin φs ; cos θs Þ

(eijj ,ei? ,

ð22Þ

ð27Þ

It is worthwhile to note that in the books [18] and [19] the incident direction is always taken along the z-axis. It is important to note that their basis vector ei? has the opposite direction by definition. Therefore the offdiagonal elements of their Jones matrices have the opposite sign as compared to our quantities. Our basis vectors are the same as in [16,17]. Also, there are the coordinate systems used for any reflection/refraction event on a facet with the normal N (Fig. 5). Here we use the right-hand triples of the unit vectors (Ijj ; I ? ; i), where i is the propagation direction before the event and Ijj and I ? are the parallel and perpendicular vectors relative to the plane containing the vectors n and N. We note that these basis vectors are the same as in [13] but are different from [19]. 2.4. Shapes and orientation of ice crystals The current version of the code considers only convex crystals with the real-valued refractive index, i.e. absorption of

Fig. 4. The coordinated systems.

Fig. 5. The coordinate systems for the reflection/refraction event.

A.V. Konoshonkin et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 164 (2015) 175–183

light is not included. For convenience, some basic shapes of ice crystals are predefined, they are hexagonal prism, bullet, droxtal, etc. A list of them is presented in the user manual [11]. Any crystal particle is determined by coordinates of its vertexes in the laboratory coordinate system. These coordinates are calculated in response to specification of the necessary sizes in microns. For example, the sizes for a hexagonal prism are its height h and the diameter d of the circle circumscribing the hexagonal facet. The initial particle orientation is indicated in the user manual (see Fig. 6a). An arbitrary particle orientation is specified by three Euler angles α, β, γ (see Fig. 6b). There are a number of their definitions. In our code, they are defined as follows: a) rotation about the zp axis through an angle α A [0, 2π]; b) rotation about the y0 axis through an angle β A[0, π]; c) rotation about the z″ axis through an angle γ A [0, 2π]. An angle of rotation is positive if the rotation is in the clockwise direction, when looking in the positive direction of the rotation axis.

3. Interface There are a lot of different applications of the light scattering problems for ice crystals. For example, scattering by either fixed or randomly oriented crystals; scattering either at some fixed direction or over all directions; and so on. It is not expedient to develop one universal program for all the applications. However, the same beamsplitting technique can be used. This situation led us to the idea of implementing beam-splitting algorithm as a library or as a software development kit (SDK), which can be easy applied to a particular light scattering problem. The SDK includes an example project of solving the light scattering by randomly oriented ice crystals. The Reference Manual is also included in the SDK. The source code of beam-splitting algorithm is written in C þ þ as free and open source software. The code uses only standard components this makes it to be a crossplatform solution and provides its compatibility to popular Integrated Development Environments (IDE's). The source code is distributed under the GNU General Public License (GPL).

Fig. 6. The initial particle orientation (a) and the Euler angles (b). Table 1 Input parameters. Number of line 1 2 3 4 5 6 7 8 9 10 11 12 13 14

179

Meanings

Particle kind. There are 5 particle types: hexagonal prism, hexagonal bullet, pyramid, tapered prism, cup. Not used for hexagonal prism. If it is chosen as 1 the lines 5 and 6 are ignored. Otherwise the tip angle is set to 56 degrees. Radius of the circle circumscribed around a hexagonal facet (μm). Half-height of the particle (μm). Height of the tip (μm) Radius of the circle circumscribed around the tip (μm) Real-valued refractive index. Number of steps around the axis of symmetry. Number of steps along the zenith angle. Number of bins in the output file «M.dat» relative to the zenith scattering angle with equal intervals. The output file supplies the nonzero elements of the Mueller matrix. Number of reflection/refraction events where 0 means only the external specular reflection 1 Means the forward and backward peaks are output to separate files, otherwise they are presented in the output file «M.dat». 1 Means discrete steps for particle orientations. 0 means that the particle orientation is chosen by use of “random()” function. 0 Means the program takes all outgoing beams into account. Any other integer number specifies the beams taken into account that are listed below in the file.

180

A.V. Konoshonkin et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 164 (2015) 175–183

The included example of solving the light scattering by a randomly oriented ice crystal is written using Qt 5.1, consequently it is a cross-platform solution, too. This solution is a single-threaded. At present, this version of the code does not take into account light absorption and nonconvex crystals. There are no fundamental restrictions for including them. 3.1. Description of the included example The example consists of the main file «main.cpp», the settings file «params.dat», the Qt Creator project file «Random.pro», the compiled binary file «random.exe», the library of beam-splitting algorithm «Lib\…» and Reference Manual «refman.pdf», which contains the class reference. All parameters to light scattering problem are set in the text file «params.dat». Each line contains one parameter which is separated from the description by double slash symbol «//». All these parameters are shown in Table 1. The solution to the problem has several steps. First, all the parameters have to be set in the settings file «params. dat». Then, the project has to be compiled and run. The settings file must be located in the same folder as the binary file («random.exe» for Windows). Once the program runs, the console window appears with the solution progress shown there. As a result of the calculation, two files will be created in the directory: «out.dat» and «M.dat». These files contain the solution to the light scattering problem as a look up table of the Mueller matrices. To solve another light scattering problem, the developer has to write down a suitable main file «main.cpp» according to the example. 3.2. Description of the main file Let us discuss the peculiarities of the main file «main.cpp». Special header file must be included in the main file. This file «\Lib\particle.hpp» contains all necessary references to other header files, function prototypes, type definitions, etc. It is useful to include other files: «\Lib\Mueller.hpp» contains the function of evaluating the Mueller matrix from the Johns matrix by Eq. (12), «\Lib\PhysMtr.hpp» contains the definitions of two-dimensional real-valued and complex-valued matrices, «\Lib\trajectory.hpp» contains the definitions of a trajectory chain.

Fig. 7. Contributions of different kinds of beams to the halo phenomena. 1 means externally reflected beams; 2 and 3 are the beams passed through the dihedral angles of 601 and 901, respectively; the heavy curve is the total phase function.

User must write down the handler function prototype as follows void Handler(Beam& bm); and add the realization of this function according to the task. The start point of using the algorithm is to create the particle. Example: Crystal* Body ¼ NULL; Body ¼new Prism(RefIndex, Radius, Halh_Height, rec_depth, i, Et, E_min, D, S_min). Several parameters are specified here: refractive index, particle size, recursion depth, incident direction i, the basis vector ei? , threshold value of the energy for omitting a beam, distance D for calculation the optical path, threshold area for omitting a beam. Next, user must specify whether the algorithm calculates the optical path or not: Body-4Phase() ¼ false; The particle is rotated at the Euler angles with the function: Body-4ChangePosition(betta, gamma, alpha). The beam-splitting process is run by the function Body-4FTforConvexCrystal(Handler); where the handler function is passed as a parameter. The function returns a value of the cross section for this particle orientation. As a result of these steps, the beam-splitting algorithm has been run. It finds out all outgoing beams, evaluates the cross section, John matrices, and optical path for each beam in the framework of the geometrical optics approximation. Then a user has to handle these beams in the handler function according to his task. The handler function takes every outgoing beam (bm) void Handler(Beam& bm) which contains all the information. The cross section, for example, CrossSection(bm); Jones matrix bm(); outgoing direction bm.r; optical path bm.lng; and so on. See Reference manual for details.

In the included example, the program rotates the particle according to the settings file, runs the beamsplitting function and summarizes the Mueller matrices

Fig. 8. The phase function for only external reflection. The solid curve is the analytical equation; the dashed and dotted curves correspond to the corrected and initial Macke's codes, respectively.

A.V. Konoshonkin et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 164 (2015) 175–183

for all outgoing beams in array. Then the array is output to the file «M.dat». For convenience, the file «main.cpp» has comments. Most of the geometric primitives, which are used in the library, are taken from the book Laszlo [20] and are significantly modified.

4. Beam-splitting versus ray-tracing The main difference between the ray-tracing and beamsplitting techniques is that the first method deals with stochastic rays while the second method operates with geometrically determined beams. In the ray-tracing

181

method, the ray stochastics leads to some additional errors and a necessity for their statistical estimates. So an advantage of a beam-splitting algorithm is that any obtained solution does not include random input and output data and the solution becomes exact in this meaning. Operating with the whole beams, we get a number of advantages. The beam-splitting algorithm allows anybody to take such wave phenomena as diffraction and interference into account providing the exact boundaries and phases of beams. The information concerning beam crosssections makes it easy to separate and reject the beams with small transverse sizes whose contributions are negligible. There is an ability to consider separately the beams that are of interest. For example, as known, the halo 221 and

Fig. 9. Non-zero elements of the phase matrix calculated with the ray-tracing and beam-splitting codes for the depth of recursion of 10. The solid curve corresponds to the beam-splitting code; the dashed and dotted curves correspond to the corrected and initial Macke's codes, respectively.

182

A.V. Konoshonkin et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 164 (2015) 175–183

461 are created by the beams passed through the dihedral angles of 601 and 901, respectively. Contributions of these beams to the total phase function are illustrated in Fig.7 by the curves 2 and 3. The contributions shown in Fig. 7 prove that while the halo 221 is formed predominately by only the beams that pass the dihedral angle 601, the halo 461 is created more complicatedly. 5. Comparison of the beam-splitting and ray-tracing simulations Generally speaking, the beam-splitting codes are considerably faster since only vertexes of the beams are traced. However, our present version of the beamsplitting code operates recursively, i.e. every reflected/ refracted beam is used in the next iteration as an initial one. In the ray-tracing techniques, one ray can be divided twice while a beam can be split to several beams. As a result of these two circumstances the speed of the beamsplitting simulation is decreased faster with increasing of the recursion depth. Therefore there is a threshold when the ray-tracing algorithm becomes faster. In this section, we illustrate this conclusion comparing the numerical data obtained by our beam-splitting code [11] and by the ray-tracing code developed by Macke [21] By the way, we noticed a small mistake in the code written by Macke. Namely, in the code [4], incident rays are emitted from a rectangular bounded the particle projection. With a change of particle orientation, the area of the rectangular is variable. This area variation must be taken into account. However, this code includes this variation only to the cross-section evaluation (line 376) but it is omitted by mistake to the Mueller-matrix evaluation (line 267). To prove this statement, we have considered a simple case of only externally reflected rays. Here there is an analytical solution to the phase matrix [14]. In particular, the first element is equal to     2 2   B F ? θ=2 þ F J θ=2 P 11 θ ¼ ; ð28Þ 2 16π where B is the total area of a crystal surface, F jj and F ? are the reflection Fresnel coefficients, θ is the scattering angle. The phase function is defined in [21] as     2P 11 θ   : ð29Þ P θ ¼R P 11 θ sin θdθ In Fig. 8, the solid curve corresponds to Eq. (29), the dashed and dotted curves are simulated by the corrected and original Macke's code, respectively, that are simulated for only specular reflection of randomly oriented hexagonal ice crystal particle with size of 80  200 μm and refractive index of 1.31. We see that the corrected code agrees with the analytical solution very well. Now we can compare the beam-splitting algorithm with ray-tracing. Fig. 9 shows all non-zero elements of the phase matrix calculated with the ray-tracing and beam-splitting codes for the depth of recursion of 10. The number of rays in the ray-tracing code was assumed 1000. The number of crystal orientations was 1,000,000.

Fig. 10. The execution time versus the recursion depth. The solid curve corresponds to the beam-splitting code; the dashed and dotted curves correspond to the Macke's code with 1000 and 10,000 initial rays, respectively.

We see the excellent agreement between the beamsplitting and corrected ray-tracing codes. Let us compare the execution times of both algorithms. In the ray-tracing code, the execution time predominantly depends on two parameters: the recursion depth, i.e. the number of the reflection/refraction events, and the number of the incident rays. In the beam-splitting code, only the recursion depth remains. Fig. 10 presents the execution time depending on the recursion depth obtained for the randomly oriented test crystal. We see that the beamsplitting code essentially faster at small recursion depths. We also see that the execution time increases almost linearly for the ray-tracing and exponentially for the beam-splitting code. As a result, there is a threshold where the ray-tracing code becomes faster. In the case of 1000 incident rays, the threshold is 6. Usually, anyone takes more than 10,000 rays with the recursion depth less than 10. In this case, the beam-splitting code is preferable. 6. Conclusions We have presented the free and open-source beamsplitting code solving the problem of light scattering by faceted particles in the framework of geometric optics. Comparison with the conventional ray-tracing code by Macke showed an excellent agreement. The light scattered by faceted particles reveals a number of specific features that are fully taken into account in our code. Therefore this code is faster than conventional ray-tracing techniques. This code has the following advantages. It excludes stochastics in input and output data inherent to ray-tracing. Specific beams responsible for some physical phenomena like halos can be easy separated and studied. Basing on shapes and energies of the whole beams it is easy to apply any physical limitations and to expand the code with the wave phenomena like diffraction and interference. The code is written in C þ þ as a separate library that can be easy included in the third party programs. The SDK

A.V. Konoshonkin et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 164 (2015) 175–183

also includes Reference Manual and an example, for convenience. This code is cross-platform and singlethreated. Acknowledgments The authors are grateful to Dr. I.A. Grishin for his contribution to development of the first version of the code. This work is supported by the Russian Foundation for Basic Research under Grants nos. 15-05-06100 and 15-5553081, and partly by the Russian Science Foundation (Agreement no. 14-27-00022). A. Konoshonkin acknowledges the support of the RF President Grant no. MK6680.2015.5 and the RF Ministry of Education and Science under the program rising competitiveness of the TSU. References [1] Mishchenko MI, Travis LD, Mackowski DW. T-matrix computations of light scattering by nonspherical particles: a review. J Quant Spectrosc Radiat Transf 1996;55:535–75. [2] Yang P, Liou KN. Finite-difference time domain method for light scattering by small ice crystals in three-dimensional space. J Opt Soc Am A 1996;13:2072–85. [3] Yurkin MA, Maltsev VP, Hoekstra AG. The discrete dipole approximation for simulation of light scattering by particles much larger than the wavelength. J Quant Spectrosc Radiat Transf 2007;106: 546–57. [4] Ray-tracing code for light scattering at non-spherical particles containing scattering and absorbing inclusions. Available from: 〈http://tools.tropos.de/〉. [5] Popov AA. New method for calculating the characteristics of light scattering by spatially oriented atmospheric crystals. Proc SPIE 1996;2822:186–94.

183

[6] Del Guasta M. Simulation of LIDAR returns from pristine and deformed hexagonal ice prisms in cold cirrus by means of face tracing. J Geophys Res 2001;106(12):589–602. [7] Romashov DN. Light scattering by hexagonal water ice crystals. Atmos Ocean Opt 2001;14:102–10. [8] Borovoi AG, Grishin IA. Scattering matrices for large ice crystal particles. J Opt Soc Am 2003;A20:2071–80. [9] Borovoi AG, Kustova NV, Oppel UG. Light backscattering by hexagonal ice crystal particles in the geometrical optics approximation. Opt Eng 2005;44:071208. [10] Bi L, Yang P, Kattawar GW, Hu Y, Baum BA. Scattering and absorption of light by ice particles: solution by a new physical-geometric optics hybrid method. J Quant Spectrosc Radiat Transf 2011;112:1492–508. [11] Beam Splitting algorithm. Branch: geometrical optics. Available from: 〈https://github.com/sasha-tvo/Beam-Splitting〉. [12] Jackson JD. Classical electrodynamics. 3rd ed.New York: John Wiley & Sons; 1999. [13] Born M, Wolf E. Principles of Optics.Oxford: Pergamon Press; 1959. [14] Borovoi AG, Kustova NV, Oppel UG. Light backscattering by hexagonal ice crystal particles in the geometrical optics approximation. Opt Eng 2005;44:071208. [15] Borovoi A, Konoshonkin A, Kustova N. The physical-optics approximation and its application to light backscattering by hexagonal ice crystals. J Quant Spectrosc Radiat Transf 2014;146:181–9. [16] Light scattering by nonspherical particles: theory, measurements, and applications. In: Mishchenko MI, Hovenier JW, Travis LD, editors. San Diego: Academic Press; 2000. [17] Mishchenko MI, Travis LD, Lacis AA. Scattering, absorption, and emission of light by small particles.Cambridge: Cambridge University Press; 2002. [18] van de Hulst HC. Light scattering by small particles.New York: Dover; 1981. [19] Bohren CF, Huffman DR. Absorption and scattering of light by small particles.New York: Wiley; 1983. [20] Laszlo MJ. Computational geometry and computer graphics in C þ þ. Lebanon: Prentice-Hall; 1995. [21] Macke A, Mueller J, Raschke E. Single scattering properties of atmospheric ice crystals. J Atmos Sci 1996;53:2813–25.