Physical-geometric optics method for large size faceted particles

Bingqiang Sun; Ping Yang; George W. Kattawar; Xiaodong Zhang

doi:10.1364/OE.25.024044

1. Introduction

The single-scattering properties of dielectric particles are fundamental to the study of radiative transfer and remote sensing [1]. The governing equations are Maxwell’s equations in conjunction with appropriate boundary conditions [2]. For homogeneous spherical particles, solutions can be analytically obtained using the Lorenz-Mie theory (LMT) [3–5]. The generalized Lorenz-Mie theory is given for a case where the incident wave is not a plane-wave, such as a Gaussian beam or other beam forms [6]. The multiple-scattering T-matrix method is given in terms of a synergistic combination of the LMT and addition theorem [7–9] when there are multiple spheres [10–13]. As a semi-analytical solution, T-matrix method stems from the expansion of vector spherical wave functions. Two effective realizations of the T-matrix method are the extended boundary condition method (EBCM) [14–16] and the invariant-imbedding T-matrix (IITM) method [17, 18]. Furthermore, the random orientation condition for a scattering particle can be analytically implemented in conjunction with the T-matrix method [19, 20]. In addition to the aforementioned analytical and semi-analytical methods, a number of methods that are numerically rigorous, such as the finite-difference time-domain (FDTD) [21, 22] and the pseudo-spectral time domain (PSTD) [23] in conjunction with the PML conditions [24], the discrete dipole approximation (DDA) [25–27], and the multilevel fast multi-pole method [28] are widely used to compute the single-scattering properties of particles. For these methods, the random-orientation condition for the single-scattering properties has to be obtained by averaging the corresponding computational results over a finite number of specific orientations [29]. The computational time and the computer memory usage associated with the computation of the single-scattering properties of non-spherical particles increases exponentially with the size of the scattering particle. Correspondingly, the maximum size of applicability of the aforementioned methods is quite limited. Thus, the geometric optics method is inevitably required in the case of light scattering by moderate and large particles with the size parameters beyond the resonance regime where the incident wavelength is comparable to the particle size.

The conventional geometric optics method (CGOM) uses an enormous number of rays to trace the electromagnetic near field using Snell’s law and the Fresnel Formulas but ignores the mapping process from the near field to the far field [30, 31]. The ray-by-ray integral (RBRI) method not only uses the ray-tracing technique to compute the electromagnetic near field but also accurately maps the near field to the far field via the following volume-integral equation [32–34]:

{{\vec{E}}^{s} (\vec{r}) |}_{r \to \infty} = \frac{\exp (i k_{0} r)}{- i k_{0} r} \frac{i k_{0}^{3}}{4 π} \underset{V}{∭} (1 - m^{2}) [\vec{E} (\vec{r}') - (\hat{r} \cdot \vec{E} (\vec{r}'))] \exp (- i k_{0} \hat{r} \cdot \vec{r}') d^{3} \vec{r}',

where the superscript ‘s’ indicates the scattered electric field and

k_{0}

is equal to

2 π / λ

with

λ

being the incident wavelength. Other notations are explained in Table 1. The integral over the volume V includes the contributions from all the rays including all order of their sub-segments. The volume-integral equation includes not only the near field

{\vec{E}}^{s} (\vec{r})

but also the scattering direction.

\vec{r}'

Subsequently, for every ray segment, a mapping computation for all required scattering directions is necessary. Even though the far-field mapping in the RBRI uses the physical optics method, the RBRI is computationally inefficient, particularly in the case of particles with random orientations. However, if particles are faceted, all rays incident on an individual facet can be analytically expressed and be treated collectively as a single beam. Consequently, the near field can be analytically obtained through the beam-tracing technique. The number of beams, dependent only on the particle shape, is significantly reduced as compared to the original number of rays in the RBRI. The new algorithm using the analytical beam-tracing process in the near-filed and the accurate mapping process in the far field is named the physical-geometric optics method to distinguish between this approach and the RBRI [35, 36]. For the physical-geometric optics method, two substantial aspects associated with the analytical near-field achievement are the beam splitting and the process associated with absorptive particles. When a beam impinges on multiple facets, a beam splitting is necessary to ensure that each sub-beam is incident only on one facet so that the analytical formulas can be employed. In addition, if a particle is absorptive, the plane waves become inhomogeneous, that is, the planes of constant phase are not in parallel with those of constant amplitude. The constant phase plane is perpendicular to the propagation direction while the constant amplitude plane for the order p = 0 is perpendicular to the particle facet and has to be traced for the orders p>0 using the boundary conditions [37]. The record of the reflection-refraction events is represented by symbol p, where p = 0 indicates the external reflection and internal refraction occurrence. Regarding a propagating vector, an absorptive vector has to be defined to account for the aforementioned inhomogeneous effect in absorptive particles. An algorithm developed by Bi et al. [35], referred to as the physical-geometric optics hybrid (PGOH) method, uses a simple geometric technique to split the original beam into sub-beams. The geometric beam-splitting technique is only applicable to homogeneous convex particles; moreover, it may substantially increase the number of sub-beams for practical application and, thus, can become inefficient in the near-to-far field mapping process. Though the PGOH defines an absorptive vector to account for the inhomogeneous effect associated with an absorptive particle, the absorptive vector is only defined for a homogeneous particle. Furthermore, the variables with respect to the absorptive vector in the PGOH may be divergent in extremely strong absorptive cases associated with large sizes.

Table 1. Symbols used in the formulas in this paper.

View Table | View all tables in this article

In this study, a general absorption vector is defined and is applicable to arbitrarily faceted particles. Correspondingly, the real parameters, including the refractive index, phase and wave vector, are systematically generalized to effective complex parameters when the particle is absorptive so that the convergence is ensured in a general case. Most importantly, the divide-and-conquer algorithm based on computer graphics [38], including translation, projection, and the polygon-clipping algorithm, are employed to split an original beam into several sub-beams that are incident on individual facets. The beam-splitting technique ensures the most efficient splitting of an original beam. Compared to the PGOH, the method reported in the paper, named the physical-geometric optics method (PGOM) to be differentiated from the PGOH, significantly reduces the computational time by a factor of 3 to 10, depending on particle shape and refractive index. Based on the techniques and the effective parameters, the PGOM can be applied to arbitrarily faceted particles.

This paper is divided into five sections. Sections 2 and 3 present theoretical framework of the amplitude scattering matrix in terms of the geometric optics and then the physical optics methods, respectively. The assessment the numerical performance of the PGOM in comparison with the PGOH and the IITM is given in section 4. The conclusions for the study are given in section 5.

2. Geometric Optics

If the incident wave is not a plane wave, it can be expanded in terms of the plane waves based on the Fourier expansion technique. For this reason, a plane wave incidence is assumed in this study. Geometric optics assumes that an incident wave can be replace by a large number of localized rays with their respective phases, amplitudes, and cross-sections and that the propagation and the reflection/refraction of a localized ray obey the same principle as those of the corresponding plane wave [4, 39]. When an incident wave is intercepted by a particle, the reflected and refracted waves can be computed using Snell’s law and the Fresnel formulas. The same reflection-refraction events are recursively repeated for the following ones. Figure 1 shows the incident planes for the reflection-refraction events with order p = 0 and p>0. The normal to the particle facets in Fig. 1 is such that, ${\hat{n}}_{p} \cdot {\hat{e}}_{p} < 0$ . Media 1 and 2 in the panel (b) of Fig. 1 are to represent the multi-layer cases; otherwise, medium 2 is air. All the symbols used in the paper are summarized in Table 1.

Fig. 1 Diagrams showing the reflection-refraction events at order p = 0 and p>0. For p>0, the superscripts ‘s’, ‘r’, and ‘t’ represent the scattered, reflected and refracted quantities. $\overset{⌢}{β} \times \overset{⌢}{α} = \overset{⌢}{e}$ Here, the subscript ‘i’ represents the incident quantities. All other subscripts correspond to the reflection-refraction order p.

Download Full Size | PDF

2.1 Beam-tracing

For a faceted particle, a bundle of rays incident on the same facet can be treated as one single beam. The beam cross-section on the facet is a polygon. The beam from the p-th order to the (p + 1)-th order is a prism whose two ends lay on the corresponding two facets. Figure 2 illustrates a systematic beam prism assuming that its two ends are pentagons.

Fig. 2 A systematic beam prism obtained from the p-th order to the (p + 1)-th order.

Download Full Size | PDF

The electric field at the first vertex of the p-th order beam cross-section can be written as follows:

[\begin{matrix} \begin{matrix} E_{p, α} & ({\vec{r}}_{p, 1}) \end{matrix} \\ \begin{matrix} E_{p, β} & ({\vec{r}}_{p, 1}) \end{matrix} \end{matrix}] = U_{p} [\begin{matrix} E_{α}^{i} \\ E_{β}^{i} \end{matrix}] \exp (i k_{0} ϕ_{p, 1}),

ϕ_{p, 1} = φ_{p, 1} + i ψ_{p, 1} .

where

U_{p}

indicates the Fresnel coefficients in the previous orders and the corresponding rotations; the real and imaginary components of the effective phase

ϕ_{p, 1}

are calculated by the following recurrence relations:

φ_{0, 1} = {\hat{e}}^{i} \cdot {\vec{r}}_{0, 1}, ψ_{0, 1} = 0,

φ_{p, 1} = φ_{p - 1, 1} + N_{r, p} | {\vec{r}}_{p, 1} - {\vec{r}}_{p - 1, 1} |,

ψ_{p, 1} = ψ_{p - 1, 1} + N_{i, p} | {\vec{r}}_{p, 1} - {\vec{r}}_{p - 1, 1} |,

N_{p} = N_{r, p} + i N_{i, p} .

where

N_{p}

is the effective refractive index. The Fresnel formula and the effective refractive index are given in Chang et al. [37] and partially discussed by Bi et al [40].

The effective wave vector can be defined in terms of a general absorptive vector as follows:

{\vec{k}}_{p} = k_{0} (N_{r, p} {\hat{e}}_{p} + {\vec{A}}_{p}),

{\vec{A}}_{0} = 0, {\vec{A}}_{p} = {\vec{A}}_{p - 1} + (N_{i, p} - {\hat{e}}_{p} \cdot {\vec{A}}_{p - 1}) \frac{{\hat{n}}_{p - 1}}{{\hat{e}}_{p} \cdot {\hat{n}}_{p - 1}} .

The absorption vector

{\vec{A}}_{p}

depicting the inhomogeneity of absorptive particles is given in terms of a recurrence relation. The electric field in any position of the p-th order beam cross-section is

[\begin{matrix} \begin{matrix} E_{p, α} \\ E_{p, β} \end{matrix} & \begin{matrix} ({\vec{r}}_{p, 1} + {\vec{ω}}_{p}) \\ ({\vec{r}}_{p, 1} + {\vec{ω}}_{p}) \end{matrix} \end{matrix}] = U_{p} [\begin{matrix} E_{α}^{i} \\ E_{β}^{i} \end{matrix}] \exp [i (k_{0} ϕ_{p, 1} + {\vec{k}}_{p} \cdot {\vec{ω}}_{p})],

{\vec{ω}}_{p}

is a position vector on the p-th order beam cross-section, originating at the first vertex. Correspondingly, any position in the beam prism is

[\begin{matrix} \begin{matrix} E_{p, α} \\ E_{p, β} \end{matrix} & \begin{matrix} (\vec{r}) \\ (\vec{r}) \end{matrix} \end{matrix}] = U_{p} [\begin{matrix} E_{α}^{i} \\ E_{β}^{i} \end{matrix}] \exp [i (k_{0} ϕ_{p, 1} + {\vec{k}}_{p} \cdot {\vec{ω}}_{p} + {\vec{k}}_{p + 1} \cdot {\hat{e}}_{p + 1} l)],

\vec{r} = {\vec{r}}_{p, 1} + {\vec{ω}}_{p} + l {\hat{e}}_{p + 1} .

The energy in the p-th order beam cross-section can be obtained

F_{p} = \frac{N_{r, p}}{2} ({| U_{p}^{11} |}^{2} + {| U_{p}^{12} |}^{2} + {| U_{p}^{21} |}^{2} + {| U_{p}^{22} |}^{2}) | {\hat{e}}_{p + 1} \cdot {\hat{n}}_{p} | {\tilde{D}}_{p},

\begin{array}{l} {\tilde{D}}_{p} = \exp (- 2 k_{0} φ_{p, 1}) \iint_{Λ_{p}} \exp (- 2 k_{0} {\vec{A}}_{p} \cdot {\vec{ω}}_{p}) d^{2} {\vec{ω}}_{p} \\ = \frac{1}{2 k_{0}^{2}} \sum_{j = 1}^{N_{ν}} [\frac{k_{0} ({\vec{r}}_{p, j + 1} - {\vec{r}}_{p, j}) \cdot ({\vec{A}}_{p} \times {\hat{n}}_{p})}{{\vec{A}}_{p} \cdot {\vec{A}}_{p} - {({\vec{A}}_{p} \cdot {\hat{n}}_{p})}^{2}}] [\frac{\exp (- 2 k_{0} φ_{p, j + 1}) - \exp (- 2 k_{0} φ_{p, j})}{- 2 k_{0} (φ_{p, j + 1} - φ_{p, j})}] . \end{array}

In the above two equations,

N_{r, p}

denotes the light speed ratio in vacuum and the scattering medium, the factor ½ is from the average of the parallel and perpendicular components,

U_{p}^{i j}, i, j = 1, 2

are the matrix elements of,

U_{p}

and

{\tilde{D}}_{p}

stands for the absorption.

2.1 Beam-splitting

When a beam in the p-th order is incident on multiple facets, a beam splitting is necessary to ensure the split sub-beam is incident only on a single facet. A systematic computer graphics method is used to implement the beam splitting procedure [38]. This approach for the beam-splitting procedure is the so-called divide-and-conquer algorithm [38]. The beam-splitting procedure is a complicated three-dimensional problem. Thus, it substantially simplifies the computation if we transform the three-dimensional problem into a two-dimensional problem in terms of some useful transformations. Specifically, the transformations include translation, rotation, and projection. After using an efficient algorithm to solve the two-dimensional problem, we use the transformations to transform the two-dimensional results to the three-dimensional solution.

We use a convex hexahedron shown in Fig. 3 as an example to show the beam-splitting technique. We suppose that a plane wave (not illustrated in Fig. 3(a)) is incident on the bottom surface. The external reflected wave goes back to the air and no further beam-splitting is necessary. However, the refracted beam (the green arrows shown in Fig. 3(a)) is incident partially on the top facet and partially on the side facets. Here the incidence on the top surface is used as an example to illustrate the standard beam splitting procedure. The incidence on the side facets can be accomplished using the following standard process. The red arrows shown in Fig. 3(a) are the local frame of reference attached to the bottom surface. Figure 3(b) shows the local frame of reference (in red arrows) and the global frame of reference (in blue arrows). The first step implements translational and rotational operations to transfer all the quantities from the global frame of reference (the blue arrows) to the local frame of reference (the red arrows) and to define the bottom surface as the x-o-y plane. The corresponding translational and rotational matrices in the affine space are:

T (x_{L}, y_{L}, z_{L}) = [\begin{matrix} 1 & 0 & 0 & - x_{L} \\ 0 & 1 & 0 & - y_{L} \\ 0 & 0 & 1 & - z_{L} \\ 0 & 0 & 0 & 1 \end{matrix}],

R = [\begin{matrix} R_{3} & 0 \\ 0 & 1 \end{matrix}],

where

x_{L}, y_{L}, z_{L}

is the origin of the local frame of reference relative to the global frame of reference and R₃ is the 3-dimensional rotational matrix in vector space. All quantities are expressed in the local frame of reference. Fig. 3(c) shows the projection process that all vertices on the top facet are projected along the opposite beam propagation direction (the green arrows) onto the bottom surface (the x-o-y plane of the local frame of reference) and the projection matrix in affine space can be written as:

M_{proj} = [\begin{matrix} 1 & 0 & - e_{x} / e_{z} & 0 \\ 0 & 1 & - e_{y} / e_{z} & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{matrix}],

where

(e_{x}, e_{y}, e_{z})

is the propagation direction (the green arrows in Fig. 3(a)). After these transformations, the points from the projection of the top vertices and the points from the bottom vertices locate on the x-o-y plane of the local frame of reference. Correspondingly, in Fig. 3(c), the green polygon is from the projection of the top facet. The bottom facet is actually the cross-section of the original beam. The green polygon is used to clip the bottom polygon and the beam with the clipped polygon is definitely incident on the top facet. The beam-splitting problem now turns into the polygon-clipping problem.

Fig. 3 Translation, rotation, and projection process.

Download Full Size | PDF

For clarity, the polygon of the original beam cross-section (the bottom facet here) and the polygon from the projection of the potentially incident facet (the top facet here) are referred to as the subject polygon and the clipping polygon, respectively., The polygon clipped off the subject polygon by the clipping polygon is referred to as the clipped polygon. An effective method, called the Weiler-Atherton algorithm [38], is used to illustrate the polygon-clipping procedure. The algorithm is robust and cannot only deal with the case where the subject and clipping polygons are both convex, such as polygons shown in the bottom facet of Fig. 3(c), but can also deal with the cases where the polygons can even be concave or even have holes in them. Figure 4 shows the schematic process of the Weiler-Atherton algorithm, where S, C, and I represent subject, clipping vertices, and intersection point, respectively. The new subject and clipping polygons are used to better illustrate the Weiler-Atherton algorithm in Fig. 4.

Fig. 4 The schematic Weiler-Atherton algorithm.

Download Full Size | PDF

The upper-left panel in Fig. 4 shows the subject and the clipping polygons and the corresponding vertices and intersection points; the upper-right panel is the subject list and the lower-right panel is the clipping list; the lower-left lists the vertices of the clipped polygon. We first need to find all intersection points between the subject polygon and the clipping polygon and divide these points into the entering points (I₁, I₃, and I₅) and leaving points (I₂, I₄, and I₆) clockwise along the subject polygon. Then, the vertices from the subject polygon and from the clipping polygon and the intersection points can be formed to create the subject and clipping cyclic linked lists, respectively. The vertices of the clipped polygon are started from one of the entering points (here I₁ in the upper-right panel), named as the starting entering point, in the subject linked list and followed in a clockwise manner to points after the entering point until a leaving point (here I₂ in the upper-right panel) is encountered. Then, the vertices of the clipped polygon have to switch to the corresponding leaving point (here I₂ in the lower-right panel) in the clipping list. After the leaving point is selected, the preceding procedure needs to be continued with the remaining points (here only I₃ in the lower-right panel) in a clockwise manner until an entering point is encountered. The process of finding the vertices of the clipped polygon is repeated by going back to the entering point in the subject list and to the leaving point in the clipping list (here I3 and I4 in the upper-right panel, then I4 and I5 in the lower-right panel, then I5, S4, and I6 in the upper-right panel, the I6 and I1 in the lower-right panel, and then I1 in the upper-right panel). The process can be stopped unless the starting entering point (here I1 in the upper-right panel) is encountered. The clipped polygon is obtained by the aforementioned algorithm. If there are still no encountered entering points, the same processes are needed to find the next clipped polygon until all the entering points are exhausted. Even though there are three entering points in Fig. 4, staring at any one of the entering points one can encounter all the entering points so that only one clipped polygon is obtained as shown in the lower-left panel. If no intersection is obtained, a ray-casting algorithm [38] is employed to determine if the clipping polygon is inside of the subject polygon or not so that the clipping polygon is the clipped polygon or not.

The clipped polygons are transformed using the reverse rotational and translational processes back to the three-dimensional polygons in the global frame of reference. The transformed clipped polygons are the required sub-beams, which are definitely incident on individual facets.

3. Physical Optics: Mapping from Near Field to Far Field

3.1 Amplitude Scattering Matrix

The near electric field in Eq. (1) can be approximately obtained in terms of the beam tracing process and the near and far fields can be decomposed into parallel and perpendicular components. The Eq. (1) can be represented as a summation over the scattering order p:

[\begin{matrix} E_{α}^{s} \\ E_{β}^{s} \end{matrix}] = \frac{\exp (i k_{0} r)}{- k_{0} r} \frac{i k^{3}}{4 π} \sum_{p = 0}^{\infty} \underset{V_{p}}{∭} (1 - m_{p + 1}^{2}) K_{p} [\begin{matrix} \begin{matrix} E_{p, α} & (\vec{r}') \end{matrix} \\ \begin{matrix} E_{p, β} & (\vec{r}') \end{matrix} \end{matrix}] \exp (- i k_{0} \hat{r} \cdot \vec{r}') d^{3} \vec{r}',

K_{p} = [\begin{matrix} {\hat{α}}^{s} \cdot {\hat{α}}_{p} & {\hat{α}}^{s} \cdot {\hat{β}}_{p} \\ {\hat{β}}^{s} \cdot {\hat{α}}_{p} & {\hat{β}}^{s} \cdot {\hat{β}}_{p} \end{matrix}] .

The matrix

K_{p}

illustrates the transformation to different frames of reference and the near electric field is given in Eqs. (2), (10), and (11). Substituting Eq. (11) into Eq. (18), a compact form can be given by:

[\begin{matrix} E_{α}^{s} \\ E_{β}^{s} \end{matrix}] = \frac{\exp (i k_{0} r)}{- k_{0} r} \sum_{p = 0}^{\infty} (1 - m_{p + 1}^{2}) q_{p} (\hat{r}) K_{p} U_{p} [\begin{matrix} E_{α}^{i} \\ E_{β}^{i} \end{matrix}],

where the integral in Eq. (18) is included in the integral factor

q_{p}

:

\begin{array}{l} q_{p} (\hat{r}) = \frac{i k_{0}^{3}}{4 π} \exp (i k_{0} ϕ_{p, 1}) \underset{V_{p}}{∭} \exp [i ({\vec{k}}_{p} \cdot {\vec{ω}}_{p} + {\vec{k}}_{p + 1} \cdot {\hat{e}}_{p + 1} l)] \exp (- i k_{0} \hat{r} \cdot \vec{r}') d^{3} \vec{r}' \\ = \frac{i k_{0}^{3}}{4 π} \exp (i k_{0} Φ_{p, 1}) \iint_{Λ_{p}} d^{2} {\vec{ω}}_{p} {| {\hat{e}}_{p + 1} \cdot {\hat{n}}_{p} | \exp (i {\vec{κ}}_{p} \cdot {\vec{ω}}_{p}) \int_{0}^{[{\vec{r}}_{p + 1, 1} - {\vec{r}}_{p, 1}] + \frac{{\vec{ω}}_{p + 1} \cdot {\hat{n}}_{p}}{{\hat{e}}_{p + 1} \cdot {\hat{n}}_{p}}} d l \exp [i ({\vec{κ}}_{p + 1} \cdot {\hat{e}}_{p + 1}) l]}, \end{array}

Φ_{p, j} = ϕ_{p, j} - \hat{r} \cdot {\vec{r}}_{p, j},

{\vec{κ}}_{p} = {\vec{k}}_{p} - k_{0} \hat{r},

k_{0} Φ_{p + 1, j} = k_{0} Φ_{p, j} + {\vec{κ}}_{p + 1} \cdot ({\vec{r}}_{p + 1, j} - {\vec{r}}_{p, j}),

k_{0} Φ_{p, j + 1} = k_{0} Φ_{p, j} + {\vec{κ}}_{p} \cdot ({\vec{r}}_{p, j + 1} - {\vec{r}}_{p, j}) .

Using the following identities

{\vec{ω}}_{p + 1} = {\vec{ω}}_{p} + \frac{{\vec{ω}}_{p + 1} \cdot {\hat{n}}_{p}}{{\hat{e}}_{p + 1} \cdot {\hat{n}}_{p}} {\hat{e}}_{p + 1},

{\vec{A}}_{p + 1} \cdot {\vec{ω}}_{p + 1} = {\vec{A}}_{p} \cdot {\vec{ω}}_{p} + N_{i, p + 1} \frac{{\vec{ω}}_{p + 1} \cdot {\hat{n}}_{p}}{{\hat{e}}_{p + 1} \cdot {\hat{n}}_{p}} .

We can write the integral factor

q_{p}

in the form of:

q_{p} (\hat{r}) = \frac{k_{0} (| {\hat{e}}_{p + 1} \cdot {\hat{n}}_{p + 1} | D_{p + 1} - | {\hat{e}}_{p + 1} \cdot {\hat{n}}_{p} | D_{p})}{{\vec{κ}}_{p + 1} \cdot {\hat{e}}_{p + 1}},

\begin{array}{l} D_{p} = \frac{k_{0}^{2}}{4 π} \exp (- i k_{0} Φ_{p, 1}) \iint_{Λ_{p}} \exp (i {\vec{κ}}_{p} \cdot {\vec{ω}}_{p}) d^{2} {\vec{ω}}_{p} \\ = \frac{k_{0}^{2}}{4 π} \sum_{j = 1}^{N_{ν}} [\frac{({\vec{r}}_{p, j + 1} - {\vec{r}}_{p, j}) \cdot (i {\vec{κ}}_{p} \times {\hat{n}}_{p})}{{\vec{κ}}_{p} \cdot {\vec{κ}}_{p} - {({\vec{κ}}_{p} \cdot {\hat{n}}_{p})}^{2}}] [\frac{\exp (i k_{0} Φ_{p, j + 1}) - \exp (i k_{0} Φ_{p, j})}{k_{0} (Φ_{p, j + 1} - Φ_{p, j})}] . \end{array}

where

\iint_{Λ} \exp (\vec{c} \cdot \vec{ω}) d^{2} \vec{ω} = \sum_{j = 1}^{N_{ν}} [\frac{\exp (\vec{c} \cdot {\vec{ν}}_{j + 1}) - \exp (\vec{c} \cdot {\vec{ν}}_{j})}{\vec{c} \cdot \vec{c} - {(\vec{c} \cdot \hat{n})}^{2}}] [\frac{(\hat{n} \times \vec{c}) \cdot ({\vec{ν}}_{j + 1} - {\vec{ν}}_{j})}{\vec{c} \cdot ({\vec{ν}}_{j + 1} - {\vec{ν}}_{j})}] .

Consequently, the far-field in Eq. (20) can be analytically expressed as:

[\begin{array}{l} E_{α}^{s} \\ E_{β}^{s} \end{array}] = \frac{\exp (i k_{0} r)}{- i k_{0} r} \sum_{p = 0}^{\infty} (1 - m_{p + 1}^{2}) [| {\hat{e}}_{p} \cdot {\hat{n}}_{p + 1} | D_{p + 1} - | {\hat{e}}_{p} \cdot {\hat{n}}_{p} | D_{p}] \frac{K_{p} U_{p}}{{\vec{κ}}_{p} \cdot {\hat{e}}_{p}} [\begin{array}{l} E_{α}^{i} \\ E_{β}^{i} \end{array}] .

The amplitude scattering matrix related to the incident and the scattering electric field can be given as:

\begin{array}{l} S = \sum_{p = 0}^{\infty} (1 - m_{p + 1}^{2}) [| {\hat{e}}_{p + 1} \cdot {\hat{n}}_{p + 1} | D_{p + 1} - | {\hat{e}}_{p + 1} \cdot {\hat{n}}_{p} | D_{p}] \frac{K_{p} U_{p} Γ}{{\vec{κ}}_{p + 1} \cdot {\hat{e}}_{p + 1}} \\ = | {\hat{e}}_{1} \cdot {\hat{n}}_{0} | D_{1} (m_{1}^{2} - 1) \frac{K_{1} U_{1} Γ}{{\vec{κ}}_{1} \cdot {\hat{e}}_{1}} + \\ \sum_{p = 1}^{\infty} D_{p} [(1 - m_{p}^{2}) \frac{K_{p - 1} U_{p - 1}}{{\vec{κ}}_{p} \cdot {\hat{e}}_{p}} | {\hat{e}}_{p} \cdot {\hat{n}}_{p} | - 〈 (1 - m_{p + 1}^{2}) \frac{K_{p} U_{p}}{{\vec{κ}}_{p + 1} \cdot {\hat{e}}_{p + 1}} | {\hat{e}}_{p + 1} \cdot {\hat{n}}_{p} | 〉] Γ . \end{array}

Γ = [\begin{matrix} {\hat{β}}^{i} \cdot {\hat{β}}^{s} & α^{i} \cdot {\hat{β}}^{s} \\ - α^{i} \cdot {\hat{β}}^{s} & {\hat{β}}^{i} \cdot {\hat{β}}^{s} \end{matrix}] .

Γ

is related to the rotation of the frame of reference. The sign “

〈 〉

” means the summation of reflected and refracted beams if the refracted beam is not the air.

3.2 Integrated Scattering Properties

The extinction and absorption cross sessions can be expressed using the same processes described in section 3.1 as follows:

\begin{array}{l} σ_{ext} = Im [\frac{k_{0}}{{| {\vec{E}}^{i} |}^{2}} \underset{V}{∭} (m^{2} - 1) (\vec{E} (\vec{r}') \cdot {\vec{E}}^{i, *} (\vec{r}')) d^{3} \vec{r}'] \\ = \frac{2 π}{k^{2}} Re [S_{11} ({\hat{e}}^{i}) + S_{22} ({\hat{e}}^{i})], \end{array}

\begin{array}{l} σ_{abs} = \frac{k_{0}}{{| {\vec{E}}^{i} |}^{2}} \underset{V}{∭} ε_{i} (\vec{E} (\vec{r}') \cdot {\vec{E}}^{*} (\vec{r}')) d^{3} \vec{r}' \\ = \frac{1}{2} \sum_{p = 0}^{\infty} N_{r, p} {| U_{p} |}^{2} [| {\hat{e}}_{p} \cdot {\hat{n}}_{p} | {\tilde{D}}_{p} - | {\hat{e}}_{p} \cdot {\hat{n}}_{p + 1} | {\tilde{D}}_{p + 1}] . \end{array}

In Eq. (34),

S_{11}

and

S_{22}

are the matrix elements of the amplitude scattering matrix given by Eq. (32). The extinction cross section is a manifestation of the optical theorem. The absorption cross section in Eq. (34) heuristically shows the difference of the incoming and outgoing energy for each order p in terms of Eq. (13). The factor ½ stands for the average of the perpendicular and parallel components.

4. Validations

With no loss in generality, two randomly generated irregular hexahedra are employed to compare the phase matrix elements in both homogeneous and inhomogeneous cases. Table 2 lists the vertices of the two hexahedra, where the radii of its smallest circumscribed sphere are unity. Figure 5 shows the two hexahedra used to do the validations. The size given in use is scaled according to the shape. All codes used here are parallelized using the Message Passing Interface (MPI).

Table 2. Vertices of two hexahedra. The two shapes are normalized so that the radii of their circumscribed spheres are unity.

View Table | View all tables in this article

Fig. 5 Three hexahedra used in the comparisons. The panel (a) is marked as hexahedron 1 and the panel (b) is hexahedron 2. The normalized vertices are given in Table 2. The panel (c) is an inhomogeneous hexahedron with hexahedron 1 included with hexahedron 2.

Download Full Size | PDF

4.1 Compared to PGOH

The physical-geometric optics hybrid (PGOH) method is developed by Bi et al. in [41] and is extensively used in the light scattering by faceted particles (e.g [42].). The PGOH is employed to verify the current method on the basis of the physical-geometric optics. The computational time of the physical-geometric optics method is highly dependent on the beam number. Compared to the PGOH, the current beam-splitting technique using computer graphics ensures the most effective splitting of the beams so that the computational time is significantly reduced.

Figure 6 shows the comparisons of the phase matrix elements of a single hexahedron, calculated by the PGOH and the PGOM under the random orientation condition. The hexahedron is shown in Fig. 5(b). The two methods give the same results. However, the PGOH is three times faster than the PGOM because of the introduction of the new beam-splitting technique that leads to the most efficient beam-splitting implementation.

Fig. 6 Comparisons of the phase matrix elements of the hexahedron 2 in random orientation calculated by the PGOH and the present PGOM. The incident wavelength is 0.658µm and the equivalent-volume radius is 8.0µm. The refractive index is 1.12 + i0.0005.

Download Full Size | PDF

Figure 7 shows the comparisons of the phase matrix elements of the inhomogeneous hexahedra, calculated by the PGOH and the PGOM under random orientation condition. The PGOH code is only for simulations for the homogeneous cases. Two cases are set up to verify the PGOM code and are effectively homogeneous cases. Consequently, the results of Fig. 7(b) are the same as the ones of Fig. 6. Even if the two cases are equivalent to a homogeneous particle, the beam-splitting process in the PGOM is still triggered because this process is only dependent on the shape and not the refractive index. The PGOM agrees perfectly with the PGOH. The two cases verify the PGOM in the beam-splitting process and the inhomogeneous scenarios.

Fig. 7 Comparisons of the phase matrix elements of inhomogeneous hexahedra, calculated by the PGOH and the PGOM under the random orientation condition. The incident wavelength is 0.658 µm. The inner layer is hexahedron 1 and the equivalent-volume radius is 4.68 µm while the outer layer is hexahedron 2 and the corresponding radius is 8.0 µm. The refractive indices of the hexahedra 1 and 2 are 1.12 + i0.0005 and 1.0 for panel (a) and 1.12 + i0.0005 and 1.12 + i0.0005 for panel (b).

Download Full Size | PDF

4.2 Comparison to highly accurate IITM method

The invariant-imbedding T-matrix (IITM) method is an accurate semi-analytical method. It uses the vector spherical wave functions to expand the incident and the scattered fields. The T-matrix that transforms the incident field to the scattered field can be written as a matrix Riccati equation or an iterative equation. The equation can be solved in terms of the invariant imbedding technique. It was pioneered by Johnson [17] and the state-of-the-art numerical implementation of it has been accomplished by Bi et al. [18]. The IITM and the EBCM are the two most efficient realizations for obtaining the T-matrix of a scattering particle. The EBCM is basically a surface-integral algorithm so that the T-matrix method using the EBCM is not stable for large particles or for a particle with a large aspect ratio. The IITM, however, is a volume-integral algorithm and the imbedding process for each iteration is within the same radius so that the IITM is extremely stable for large particles. For example, the phase matrix of a spheroid with size parameters of its semi-major and semi-minor axes around 445 and 390 can be obtained using the IITM [43]. Furthermore, the IITM can be applied to particles without any symmetry and inhomogeneous as well [20]. Due to its stability, the process involving T-matrix iteration is parallelized using MPI to terminate the memory limitation for T-matrix computation [44]. A significant advantage of the T-matrix method is that the random orientation results can be automatically obtained once the T-matrix is given. The formulas for particles with rotational symmetry are given by Mishchenko [19] and have been generalized to particles without any symmetry by Bi and Yang [20]. In order to obtain the single-scattering properties of large particles under random orientation, the random orientation process in the IITM is parallelized using MPI by the current authors so that the memory limitation toward the IITM is ultimately terminated. The IITM code is used as a benchmark in this study to validate the accuracy of the PGOM.

Figure 8 shows the comparisons of phase matrix elements of hexahedron 2, shown in Fig. 5(b), calculated by the PGOM and the IITM under the random orientation condition. The inset of the P₁₁ panel shows the same forward scattering between the PGOM and the IITM while the CGOM shows some discrepancies because the phase interference between the diffraction and the transmitted rays is automatically taken into consideration by the PGOM. The CGOM and the PGOM show the similar accuracy except the forward scattering because the particle size is approaching the geometric regime. Even though the near field for the PGOM is obtained using the ray-tracing method, the results calculated by the PGOM still agree well with the ones calculated by the IITM.

Fig. 8 Comparisons of the phase matrix elements of hexahedron 2, calculated by the PGOM, the CGOM and the IITM under random orientation condition. The equivalent-volume radius is 8.0µm and the incident wavelength is 0.658µm. The refractive index is 1.12 + i0.0005.

Download Full Size | PDF

Figure 9 shows the comparisons of the phase matrix elements of a hexagonal column with a unit aspect ratio, calculated by the PGOM and the IITM under the random orientation condition. The refractive index for the column is 1.33. The results obtained by the PGOM and the IITM are in agreement even for the oscillations of elements P₁₂/P₁₁ and P₂₂/P₁₁.

Fig. 9 Comparisons of the phase matrix elements of a hexagonal column with a unit aspect ratio (defined as the height over the diameter of the circumscribed circle of the bottom hexagon), calculated by the PGOM and the IITM under the random orientation condition. The height of the column is 20.94µm and the incident wavelength is 0.658µm. The refractive index is 1.33.

Download Full Size | PDF

Figure 10 shows comparisons of the phase matrix elements of the inhomogeneous hexahedron, shown in the Fig. 5(c), calculated by the PGOM and the IITM under the random orientation condition. The incident wavelength is 0.658 µm and the radius of the circumscribed sphere of the hexahedron is roughly 14.8µm, that is, its size parameter is around 141. The PGOM agrees quite well with the IITM in spite of the fact that the PGOM is an approximate method. There are some small discrepancies in the backward directions. Since the refractive indices of the inhomogeneous hexahedron are close to unity, the right four panels show features similar to those of soft particles. These results validate the accuracy of the present physical-geometric optics method.

Fig. 10 Comparisons of the phase matrix elements of the inhomogeneous hexahedron, calculated by the PGOM and the IITM under the random orientation condition. The incident wavelength is 0.658µm. The equivalent-volume radii of hexahedron 1 and 2 are 4.6 µm and 8.0 µm. The refractive indices of the inner layer and outer layer are 1.12 + i0.0005 and 1.03 + i0.005, respectively.

Download Full Size | PDF

The total beam number determines the computational time of the PGOM. Consequently, the simulation time is increased with the increase of the number of facets of a particle. Moreover, the simulation of an inhomogeneous particle takes longer time than that of a homogeneous particle. However, the advantage of the PGOM is that its computational time depends only on the particle shape rather than the particle size. As for a particle with non-facet geometry, the light scattering simulation can be accomplished by approximately faceting the geometry. Therefore, the PGOM can be employed for particles with arbitrary geometries.

5. Conclusion

The present physical-geometric optics method has focused exclusively on faceted particles. It uses a ray-tracing technique to analytically obtain the electromagnetic near field and subsequently maps the near field to the far field. All rays incident on a specific facet can be simultaneously traced and treated as a single beam. If a beam propagates toward multiple facets, a beam-splitting technique is necessary to ensure that each split beam is incident on an individual facet. The beam-splitting process is efficiently accomplished by the use of computer graphics technique. If the medium is absorptive, the plane wave during propagation is inhomogeneous that the planes of constant phase are not in parallel with those of constant amplitude. A general absorption vector is defined to account for the absorption of the beams for any homogeneous and inhomogeneous particle. A systematic derivation from the geometric optics to the physical optics method is given in this study. Parameters, such as phase and the refractive index, are generalized to a series of effective counterparts in the general formulas, which are then reduced to the corresponding quantities in a non-absorptive medium. The application regime of the present method is generalized from a convex homogeneous particle to an inhomogeneous convex or concave particle. The assessment of the PGOM in comparison with the PGOH is performed in a homogeneous case and two inhomogeneous cases. The comparisons between the results obtained with the PGOM and the highly accurate method IITM show close agreement.

Acknowledgments

The authors thank two anonymous reviewers for constructive comments. This study was partly supported by a NSF grant (OCE-1459180) and by the endowment funds related to the David Bullock Harris Chair in Geosciences at the College of Geosciences, Texas A&M University. A major portion of the computations involved in this study was performed with the advanced computing resources provided by Texas A&M High Performance Research Computing.

References and links

1. K. N. Liou, An Introduction to Atmospheric Radiation (Academic Press, 2002), Vol. 84.

2. J. D. Jackson, Classical Electrodynamics (John Wiley & Sons, 2007).

3. G. Mie, “Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen,” Ann. Phys. 330(3), 377–445 (1908). [CrossRef]

4. H. C. van de Hulst, Light Scattering by Small Particles (Dover Publications, 1981).

5. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley-Vch, 2008).

6. G. Gouesbet and G. Gréhan, Generalized Lorenz-Mie Theories (Springer, 2011), Vol. 31.

7. B. Friedman and J. Russek, “Addition theorems for spherical waves,” Q. Appl. Math. 12(1), 13–23 (1954). [CrossRef]

8. S. Stein, “Addition theorems for spherical wave functions,” Q. Appl. Math. 19(1), 15–24 (1961). [CrossRef]

9. O. R. Cruzan, “Translational addition theorems for spherical vector wave functions,” Q. Appl. Math. 20(1), 33–40 (1962). [CrossRef]

10. J. H. Bruning and Y. Lo, “Multiple scattering by spheres,” (University of Illinois at Urbana-Champaign, 1969).

11. K. A. Fuller and G. W. Kattawar, “Consummate solution to the problem of classical electromagnetic scattering by an ensemble of spheres. I: Linear chains,” Opt. Lett. 13(2), 90–92 (1988). [CrossRef] [PubMed]

12. D. W. Mackowski and M. I. Mishchenko, “Calculation of the T matrix and the scattering matrix for ensembles of spheres,” J. Opt. Soc. Am. A 13(11), 2266–2278 (1996). [CrossRef]

13. Y. L. Xu, “Electromagnetic scattering by an aggregate of spheres,” Appl. Opt. 34(21), 4573–4588 (1995). [CrossRef] [PubMed]

14. P. Waterman, “Matrix formulation of electromagnetic scattering,” Proc. IEEE 53(8), 805–812 (1965). [CrossRef]

15. P. Waterman, “Symmetry, unitarity, and geometry in electromagnetic scattering,” Phys. Rev. D Part. Fields 3(4), 825–839 (1971). [CrossRef]

16. M. I. Mishchenko and L. D. Travis, “Capabilities and limitations of a current FORTRAN implementation of the T-matrix method for randomly oriented, rotationally symmetric scatterers,” J. Quant. Spectrosc. Radiat. Transf. 60(3), 309–324 (1998). [CrossRef]

17. B. R. Johnson, “Invariant imbedding T matrix approach to electromagnetic scattering,” Appl. Opt. 27(23), 4861–4873 (1988). [CrossRef] [PubMed]

18. L. Bi, P. Yang, G. W. Kattawar, and M. I. Mishchenko, “Efficient implementation of the invariant imbedding T-matrix method and the separation of variables method applied to large nonspherical inhomogeneous particles,” J. Quant. Spectrosc. Radiat. Transf. 116, 169–183 (2013). [CrossRef]

19. M. Mishchenko, “Light scattering by randomly oriented axially symmetric particles,” J. Opt. Soc. Am. A 8(6), 871–882 (1991). [CrossRef]

20. L. Bi and P. Yang, “Accurate simulation of the optical properties of atmospheric ice crystals with the invariant imbedding T-matrix method,” J. Quant. Spectrosc. Radiat. Transf. 138, 17–35 (2014). [CrossRef]

21. K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propagation 14(3), 302–307 (1966). [CrossRef]

22. P. Yang and K. Liou, “Finite-difference time domain method for light scattering by small ice crystals in three-dimensional space,” J. Opt. Soc. Am. A 13(10), 2072–2085 (1996). [CrossRef]

23. Q. H. Liu, “The pseudospectral time-domain (PSTD) algorithm for acoustic waves in absorptive media,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control 45(4), 1044–1055 (1998). [CrossRef] [PubMed]

24. J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114(2), 185–200 (1994). [CrossRef]

25. E. M. Purcell and C. R. Pennypacker, “Scattering and absorption of light by nonspherical dielectric grains,” Astrophys. J. 186, 705–714 (1973). [CrossRef]

26. B. T. Draine, “The discrete-dipole approximation and its application to interstellar graphite grains,” Astrophys. J. 333, 848–872 (1988). [CrossRef]

27. M. A. Yurkin and A. G. Hoekstra, “The discrete dipole approximation: an overview and recent developments,” J. Quant. Spectrosc. Radiat. Transf. 106(1-3), 558–589 (2007). [CrossRef]

28. J. Song and W. C. Chew, “Multilevel fast‐multipole algorithm for solving combined field integral equations of electromagnetic scattering,” Microw. Opt. Technol. Lett. 10(1), 14–19 (1995). [CrossRef]

29. M. I. Mishchenko and M. A. Yurkin, “On the concept of random orientation in far-field electromagnetic scattering by nonspherical particles,” Opt. Lett. 42(3), 494–497 (2017). [CrossRef] [PubMed]

30. Q. Cai and K. N. Liou, “Polarized light scattering by hexagonal ice crystals: theory,” Appl. Opt. 21(19), 3569–3580 (1982). [CrossRef] [PubMed]

31. A. Macke, J. Mueller, and E. Raschke, “Single scattering properties of atmospheric ice crystals,” J Atmo. Sci. 53, 2813–2825 (1996).

32. P. Yang and K. Liou, “Light scattering by hexagonal ice crystals: solutions by a ray-by-ray integration algorithm,” J. Opt. Soc. Am. A 14(9), 2278–2289 (1997). [CrossRef]

33. B. Sun, G. W. Kattawar, P. Yang, M. S. Twardowski, and J. M. Sullivan, “Simulation of the scattering properties of a chain-forming triangular prism oceanic diatom,” J. Quant. Spectrosc. Radiat. Transf. 178, 390–399 (2016). [CrossRef]

34. G. Tang, P. Yang, B. Sun, R. L. Panetta, and G. W. Kattawar, “Enhancement of the computational efficiency of the near-to-far field mapping in the finite-difference method and ray-by-ray method with the fast multi-pole plane wave expansion approach,” J. Quant. Spectrosc. Radiat. Transf. 176, 70–81 (2016). [CrossRef]

35. L. Bi, P. Yang, G. W. Kattawar, Y. Hu, and B. A. Baum, “Scattering and absorption of light by ice particles: Solution by a new physical-geometric optics hybrid method,” J. Quant. Spectrosc. Radiat. Transf. 112(9), 1492–1508 (2011). [CrossRef]

36. A. Borovoi, A. Konoshonkin, and N. Kustova, “The physical-optics approximation and its application to light backscattering by hexagonal ice crystals,” J. Quant. Spectrosc. Radiat. Transf. 146, 181–189 (2014). [CrossRef]

37. P. C. Chang, J. Walker, and K. Hopcraft, “Ray tracing in absorbing media,” J. Quant. Spectrosc. Radiat. Transf. 96(3-4), 327–341 (2005). [CrossRef]

38. J. F. Hughes and J. D. Foley, Computer Graphics: Principles and Practice (Pearson Education, 2014).

39. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (Elsevier, 2013).

40. L. Bi and P. Yang, “Improved ice particle optical property simulations in the ultraviolet to far-infrared regime,” J. Quant. Spectrosc. Radiat. Transf. 189, 228–237 (2017). [CrossRef]

41. L. Bi and P. Yang, “Physical-geometric optics hybrid methods for computing the scattering and absorption properties of ice crystals and dust aerosols,” in Light Scattering Reviews 8 (Springer, 2013), pp. 69–114.

42. G. Xu, B. Sun, S. D. Brooks, P. Yang, G. W. Kattawar, and X. Zhang, “Modeling the inherent optical properties of aquatic particles using an irregular hexahedral ensemble,” J. Quant. Spectrosc. Radiat. Transf. 191, 30–39 (2017). [CrossRef]

43. B. Sun, G. W. Kattawar, P. Yang, and K. F. Ren, “Rigorous 3-D vectorial complex ray model applied to light scattering by an arbitrary spheroid,” J. Quant. Spectrosc. Radiat. Transf. 179, 1–10 (2016). [CrossRef]

44. J. Liu, P. Yang, and K. Muinonen, “Dust-aerosol optical modeling with Gaussian spheres: Combined invariant-imbedding T-matrix and geometric-optics approach,” J. Quant. Spectrosc. Radiat. Transf. 161, 136–144 (2015). [CrossRef]

m	Complex refractive index (the imaginary part is positive)
N	Effective complex refractive index (the imaginary part is positive)
φ	Real phase
ψ	Imaginary phase (amplitude attenuation factor)
ϕ	Complex phase
Φ	Complex phase difference
$\vec{k}$	Complex wave vector
$\vec{κ}$	Complex wave vector difference
$\vec{r}$	Position vector in any position, originating from the Laboratory frame origin
$\vec{ω}, \vec{υ}$	Position vector in a beam cross-section, originating from the first vertex
$\overset{⌢}{r}$	Unit vector in scattering direction
$\overset{⌢}{n}$	Outward normal to a plane
$\overset{⌢}{e}$	Unit vector normal to the plane of constant phase
$\overset{⌢}{α}, \overset{⌢}{β}$	Unit vectors parallel and perpendicular to the incident plane, respectively
$\vec{A}$	Vector normal to the plane of constant amplitude
p	Reflection-refraction order
N_v	Total vertex number of a beam cross-section
j	Vertex enumeration of a beam cross-section in subscript
i	Imaginary unit in expression; incident quantities in superscript or subscript

	Hexahedron 1	Hexahedron 2
V1	(−0.145435, −0.108740, 0.255752)	(−0.147164, −0.096923, 0.268810)
V2	(0.036896, −0.197196, 0.200786)	(0.075715, −0.123441, 0.335384)
V3	(0.032048, 0.211640, 0.442543)	(0.178159, 0.212427, 0.456195)
V4	(−0.106788, 0.206650, 0.441630)	(−0.255321, 0.207136, 0.311977)
V5	(−0.265424, −0.077406, −0.447325)	(−0.126923, −0.265593, −0.322892)
V6	(−0.245286, 0.159181, −0.359473)	(−0.300638, 0.190611, −0.421459)
V7	(0.303837, 0.190634, −0.135239)	(0.350472, 0.196127, −0.314243)
V8	(0.256585, −0.320537, −0.326506)	(0.145396, −0.306316, −0.271473)

m	Complex refractive index (the imaginary part is positive)
N	Effective complex refractive index (the imaginary part is positive)
φ	Real phase
ψ	Imaginary phase (amplitude attenuation factor)
ϕ	Complex phase
Φ	Complex phase difference
$\vec{k}$	Complex wave vector
$\vec{κ}$	Complex wave vector difference
$\vec{r}$	Position vector in any position, originating from the Laboratory frame origin
$\vec{ω}, \vec{υ}$	Position vector in a beam cross-section, originating from the first vertex
$\overset{⌢}{r}$	Unit vector in scattering direction
$\overset{⌢}{n}$	Outward normal to a plane
$\overset{⌢}{e}$	Unit vector normal to the plane of constant phase
$\overset{⌢}{α}, \overset{⌢}{β}$	Unit vectors parallel and perpendicular to the incident plane, respectively
$\vec{A}$	Vector normal to the plane of constant amplitude
p	Reflection-refraction order
N_v	Total vertex number of a beam cross-section
j	Vertex enumeration of a beam cross-section in subscript
i	Imaginary unit in expression; incident quantities in superscript or subscript

	Hexahedron 1	Hexahedron 2
V1	(−0.145435, −0.108740, 0.255752)	(−0.147164, −0.096923, 0.268810)
V2	(0.036896, −0.197196, 0.200786)	(0.075715, −0.123441, 0.335384)
V3	(0.032048, 0.211640, 0.442543)	(0.178159, 0.212427, 0.456195)
V4	(−0.106788, 0.206650, 0.441630)	(−0.255321, 0.207136, 0.311977)
V5	(−0.265424, −0.077406, −0.447325)	(−0.126923, −0.265593, −0.322892)
V6	(−0.245286, 0.159181, −0.359473)	(−0.300638, 0.190611, −0.421459)
V7	(0.303837, 0.190634, −0.135239)	(0.350472, 0.196127, −0.314243)
V8	(0.256585, −0.320537, −0.326506)	(0.145396, −0.306316, −0.271473)

Physical-geometric optics method for large size faceted particles

Abstract

1. Introduction

2. Geometric Optics

2.1 Beam-tracing

2.1 Beam-splitting

3. Physical Optics: Mapping from Near Field to Far Field

3.1 Amplitude Scattering Matrix

3.2 Integrated Scattering Properties

4. Validations

4.1 Compared to PGOH

4.2 Comparison to highly accurate IITM method

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (10)

Tables (2)

Equations (35)

Optics Express